glviewvol
view src/volume.cc @ 14:0d2447b9c512
did the histogram... doesn't seem to work right
| author | John Tsiombikas <nuclear@member.fsf.org> |
|---|---|
| date | Thu, 01 Jan 2015 06:36:45 +0200 |
| parents | 773f89037a35 |
| children |
line source
1 /*
2 glviewvol is an OpenGL 3D volume data viewer
3 Copyright (C) 2014 John Tsiombikas <nuclear@member.fsf.org>
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18 #include <stdio.h>
19 #include <string.h>
20 #include <ctype.h>
21 #include <alloca.h>
22 #include <errno.h>
23 #include <math.h>
24 #include "volume.h"
26 static char *strip_space(char *s);
29 Volume::~Volume()
30 {
31 }
33 int Volume::num_samples(int dim) const
34 {
35 return 0;
36 }
38 void Volume::normalf(float *norm, float x, float y, float z, float delta)
39 {
40 float dx, dy, dz;
41 dx = dy = dz = delta;
43 if(num_samples(0) > 0) {
44 // discrete volume
45 dx /= (float)num_samples(0);
46 dy /= (float)num_samples(1);
47 dz /= (float)num_samples(2);
48 }
50 norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z);
51 norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z);
52 norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz);
54 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
55 if(len != 0.0) {
56 norm[0] /= len;
57 norm[1] /= len;
58 norm[2] /= len;
59 }
60 }
62 void Volume::normali(float *norm, int x, int y, int z)
63 {
64 int sz[3];
66 if((sz[0] = num_samples(0)) <= 0) {
67 // if it's a continuous volume, just call normalf
68 normalf(norm, x, y, z);
69 return;
70 }
71 sz[1] = num_samples(1);
72 sz[2] = num_samples(2);
74 int prevx = x <= 0 ? 0 : x;
75 int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x;
76 int prevy = y <= 0 ? 0 : y;
77 int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y;
78 int prevz = z <= 0 ? 0 : z;
79 int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z;
81 norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z);
82 norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z);
83 norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz);
85 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
86 if(len != 0.0) {
87 norm[0] /= len;
88 norm[1] /= len;
89 norm[2] /= len;
90 }
91 }
93 // ---- VoxelVolume (discrete) implementation ----
94 VoxelVolume::VoxelVolume()
95 {
96 size[0] = size[1] = size[2] = 0;
98 hist_valid = false;
99 hist = 0;
100 num_hist_samples = 0;
101 }
103 VoxelVolume::~VoxelVolume()
104 {
105 for(size_t i=0; i<slices.size(); i++) {
106 slices[i].destroy();
107 }
108 }
110 bool VoxelVolume::load(const char *fname)
111 {
112 if(!fname) return false;
114 char *prefix = (char*)alloca(strlen(fname) + 1);
115 strcpy(prefix, fname);
116 char *slash = strrchr(prefix, '/');
117 if(slash) {
118 *slash = 0;
119 } else {
120 prefix = 0;
121 }
123 printf("loading volume dataset: %s\n", fname);
124 FILE *fp = fopen(fname, "r");
125 if(!fp) {
126 fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno));
127 return false;
128 }
130 char buf[256], path[300];
131 while(fgets(buf, sizeof buf, fp)) {
132 char *line = strip_space(buf);
133 sprintf(path, "%s/%s", prefix, line);
135 printf(" loading slice %d: %s\n", (int)slices.size(), path);
136 Image img;
137 if(!img.load(path)) {
138 slices.clear();
139 return false;
140 }
142 if(slices.empty()) {
143 size[0] = img.width;
144 size[1] = img.height;
145 } else {
146 if(img.width != size[0] || img.height != size[1]) {
147 fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n",
148 (int)slices.size(), line, img.width, img.height, size[0], size[1]);
149 slices.clear();
150 return false;
151 }
152 }
154 slices.push_back(img);
155 }
157 size[2] = slices.size();
158 hist_valid = false;
159 fclose(fp);
160 return true;
161 }
163 int VoxelVolume::num_samples(int dim) const
164 {
165 return size[dim];
166 }
168 static inline int clamp(int x, int low, int high)
169 {
170 return x < low ? low : (x > high ? high : x);
171 }
173 static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty,
174 float *pixels, int width, int height)
175 {
176 float v00, v01, v10, v11;
178 v00 = pixels[y0 * width + x0];
179 v01 = pixels[y1 * width + x0];
180 v10 = pixels[y0 * width + x1];
181 v11 = pixels[y1 * width + x1];
183 float v0 = v00 + (v01 - v00) * ty;
184 float v1 = v10 + (v11 - v10) * ty;
186 return v0 + (v1 - v0) * tx;
187 }
189 float VoxelVolume::valuef(float x, float y, float z) const
190 {
191 if(slices.empty()) {
192 return 0.0f;
193 }
195 float floor_x = floor(x);
196 float ceil_x = ceil(x);
197 float tx = x - floor_x;
199 float floor_y = floor(y);
200 float ceil_y = ceil(y);
201 float ty = y - floor_y;
203 float floor_z = floor(z);
204 float ceil_z = ceil(z);
205 float tz = z - floor_z;
207 int x0 = clamp(floor_x, 0, size[0] - 1);
208 int x1 = clamp(ceil_x, 0, size[0] - 1);
210 int y0 = clamp(floor_y, 0, size[1] - 1);
211 int y1 = clamp(ceil_y, 0, size[1] - 1);
213 int s0 = clamp(floor_z, 0, size[2] - 1);
214 int s1 = clamp(ceil_z, 0, size[2] - 1);
216 float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]);
217 float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]);
219 return val_s0 + (val_s1 - val_s0) * tz;
220 }
222 float VoxelVolume::valuei(int x, int y, int z) const
223 {
224 x = clamp(x, 0, size[0] - 1);
225 y = clamp(y, 0, size[1] - 1);
226 z = clamp(z, 0, size[2] - 1);
228 return slices[z].pixels[y * size[0] + x];
229 }
231 float *VoxelVolume::calc_histogram(int hist_samples)
232 {
233 if(hist && hist_samples == num_hist_samples && hist_valid) {
234 return hist;
235 }
237 int num_pixels = size[0] * size[1] * size[2];
238 if(!num_pixels) {
239 return 0;
240 }
242 delete [] hist;
243 hist = new float[hist_samples];
244 memset(hist, 0, hist_samples * sizeof *hist);
246 for(int i=0; i<size[2]; i++) {
247 float *pptr = slices[i].pixels;
248 for(int j=0; j<size[0] * size[1]; j++) {
249 int idx = (int)(*pptr++ * (float)hist_samples);
251 hist[idx] += 1.0;
252 }
253 }
255 for(int i=0; i<hist_samples; i++) {
256 hist[i] /= (float)num_pixels;
257 }
259 hist_valid = true;
260 num_hist_samples = hist_samples;
261 return hist;
262 }
264 static char *strip_space(char *s)
265 {
266 while(*s && isspace(*s)) s++;
267 if(!*s) return 0;
269 char *end = s + strlen(s) - 1;
270 while(end > s && isspace(*end)) end--;
271 end[1] = 0;
273 return s;
274 }