glviewvol
view src/volume.cc @ 7:71b479ffb9f7
curve manipulation works
| author | John Tsiombikas <nuclear@member.fsf.org> |
|---|---|
| date | Tue, 30 Dec 2014 17:28:38 +0200 |
| parents | 32c4a7160350 |
| children | 773f89037a35 |
line source
1 #include <stdio.h>
2 #include <string.h>
3 #include <ctype.h>
4 #include <alloca.h>
5 #include <errno.h>
6 #include <math.h>
7 #include "volume.h"
9 static char *strip_space(char *s);
12 Volume::~Volume()
13 {
14 }
16 int Volume::num_samples(int dim) const
17 {
18 return 0;
19 }
21 void Volume::normalf(float *norm, float x, float y, float z, float delta)
22 {
23 float dx, dy, dz;
24 dx = dy = dz = delta;
26 if(num_samples(0) > 0) {
27 // discrete volume
28 dx /= (float)num_samples(0);
29 dy /= (float)num_samples(1);
30 dz /= (float)num_samples(2);
31 }
33 norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z);
34 norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z);
35 norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz);
37 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
38 if(len != 0.0) {
39 norm[0] /= len;
40 norm[1] /= len;
41 norm[2] /= len;
42 }
43 }
45 void Volume::normali(float *norm, int x, int y, int z)
46 {
47 int sz[3];
49 if((sz[0] = num_samples(0)) <= 0) {
50 // if it's a continuous volume, just call normalf
51 normalf(norm, x, y, z);
52 return;
53 }
54 sz[1] = num_samples(1);
55 sz[2] = num_samples(2);
57 int prevx = x <= 0 ? 0 : x;
58 int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x;
59 int prevy = y <= 0 ? 0 : y;
60 int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y;
61 int prevz = z <= 0 ? 0 : z;
62 int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z;
64 norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z);
65 norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z);
66 norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz);
68 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
69 if(len != 0.0) {
70 norm[0] /= len;
71 norm[1] /= len;
72 norm[2] /= len;
73 }
74 }
76 // ---- VoxelVolume (discrete) implementation ----
77 VoxelVolume::VoxelVolume()
78 {
79 size[0] = size[1] = size[2] = 0;
80 }
82 VoxelVolume::~VoxelVolume()
83 {
84 for(size_t i=0; i<slices.size(); i++) {
85 slices[i].destroy();
86 }
87 }
89 bool VoxelVolume::load(const char *fname)
90 {
91 if(!fname) return false;
93 char *prefix = (char*)alloca(strlen(fname) + 1);
94 strcpy(prefix, fname);
95 char *slash = strrchr(prefix, '/');
96 if(slash) {
97 *slash = 0;
98 } else {
99 prefix = 0;
100 }
102 printf("loading volume dataset: %s\n", fname);
103 FILE *fp = fopen(fname, "r");
104 if(!fp) {
105 fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno));
106 return false;
107 }
109 char buf[256], path[300];
110 while(fgets(buf, sizeof buf, fp)) {
111 char *line = strip_space(buf);
112 sprintf(path, "%s/%s", prefix, line);
114 printf(" loading slice %d: %s\n", (int)slices.size(), path);
115 Image img;
116 if(!img.load(path)) {
117 slices.clear();
118 return false;
119 }
121 if(slices.empty()) {
122 size[0] = img.width;
123 size[1] = img.height;
124 } else {
125 if(img.width != size[0] || img.height != size[1]) {
126 fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n",
127 (int)slices.size(), line, img.width, img.height, size[0], size[1]);
128 slices.clear();
129 return false;
130 }
131 }
133 slices.push_back(img);
134 }
136 size[2] = slices.size();
137 fclose(fp);
138 return true;
139 }
141 int VoxelVolume::num_samples(int dim) const
142 {
143 return size[dim];
144 }
146 static inline int clamp(int x, int low, int high)
147 {
148 return x < low ? low : (x > high ? high : x);
149 }
151 static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty,
152 float *pixels, int width, int height)
153 {
154 float v00, v01, v10, v11;
156 v00 = pixels[y0 * width + x0];
157 v01 = pixels[y1 * width + x0];
158 v10 = pixels[y0 * width + x1];
159 v11 = pixels[y1 * width + x1];
161 float v0 = v00 + (v01 - v00) * ty;
162 float v1 = v10 + (v11 - v10) * ty;
164 return v0 + (v1 - v0) * tx;
165 }
167 float VoxelVolume::valuef(float x, float y, float z) const
168 {
169 if(slices.empty()) {
170 return 0.0f;
171 }
173 float floor_x = floor(x);
174 float ceil_x = ceil(x);
175 float tx = x - floor_x;
177 float floor_y = floor(y);
178 float ceil_y = ceil(y);
179 float ty = y - floor_y;
181 float floor_z = floor(z);
182 float ceil_z = ceil(z);
183 float tz = z - floor_z;
185 int x0 = clamp(floor_x, 0, size[0] - 1);
186 int x1 = clamp(ceil_x, 0, size[0] - 1);
188 int y0 = clamp(floor_y, 0, size[1] - 1);
189 int y1 = clamp(ceil_y, 0, size[1] - 1);
191 int s0 = clamp(floor_z, 0, size[2] - 1);
192 int s1 = clamp(ceil_z, 0, size[2] - 1);
194 float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]);
195 float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]);
197 return val_s0 + (val_s1 - val_s0) * tz;
198 }
200 float VoxelVolume::valuei(int x, int y, int z) const
201 {
202 x = clamp(x, 0, size[0] - 1);
203 y = clamp(y, 0, size[1] - 1);
204 z = clamp(z, 0, size[2] - 1);
206 return slices[z].pixels[y * size[0] + x];
207 }
209 static char *strip_space(char *s)
210 {
211 while(*s && isspace(*s)) s++;
212 if(!*s) return 0;
214 char *end = s + strlen(s) - 1;
215 while(end > s && isspace(*end)) end--;
216 end[1] = 0;
218 return s;
219 }