glviewvol

view src/volume.cc @ 1:cc9e0d8590e2

find changesets by author, revision, files, or words in the commit message
foo
author John Tsiombikas <nuclear@member.fsf.org>
date Sat, 27 Dec 2014 06:32:28 +0200
parents
children 701507c8238f
line source
1 #include <stdio.h>
2 #include <string.h>
3 #include <ctype.h>
4 #include <errno.h>
5 #include <math.h>
6 #include "volume.h"
7
8 static char *strip_space(char *s);
9
10
11 Volume::~Volume()
12 {
13 }
14
15 int Volume::num_samples(int dim) const
16 {
17 return 0;
18 }
19
20 void Volume::normalf(float *norm, float x, float y, float z, float delta)
21 {
22 float dx, dy, dz;
23 dx = dy = dz = delta;
24
25 if(num_samples(0) > 0) {
26 // discrete volume
27 dx /= (float)num_samples(0);
28 dy /= (float)num_samples(1);
29 dz /= (float)num_samples(2);
30 }
31
32 norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z);
33 norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z);
34 norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz);
35
36 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
37 if(len != 0.0) {
38 norm[0] /= len;
39 norm[1] /= len;
40 norm[2] /= len;
41 }
42 }
43
44 void Volume::normali(float *norm, int x, int y, int z)
45 {
46 int sz[3];
47
48 if((sz[0] = num_samples(0)) <= 0) {
49 // if it's a continuous volume, just call normalf
50 normalf(norm, x, y, z);
51 return;
52 }
53 sz[1] = num_samples(1);
54 sz[2] = num_samples(2);
55
56 int prevx = x <= 0 ? 0 : x;
57 int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x;
58 int prevy = y <= 0 ? 0 : y;
59 int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y;
60 int prevz = z <= 0 ? 0 : z;
61 int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z;
62
63 norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z);
64 norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z);
65 norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz);
66
67 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
68 if(len != 0.0) {
69 norm[0] /= len;
70 norm[1] /= len;
71 norm[2] /= len;
72 }
73 }
74
75 // ---- VoxelVolume (discrete) implementation ----
76 VoxelVolume::VoxelVolume()
77 {
78 size[0] = size[1] = size[2] = 0;
79 }
80
81 bool VoxelVolume::load(const char *fname)
82 {
83 FILE *fp = fopen(fname, "r");
84 if(!fp) {
85 fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno));
86 return false;
87 }
88
89 char buf[512];
90 while(fgets(buf, sizeof buf, fp)) {
91 char *line = strip_space(buf);
92
93 Image img;
94 if(!img.load(line)) {
95 slices.clear();
96 return false;
97 }
98
99 if(slices.empty()) {
100 size[0] = img.width;
101 size[1] = img.height;
102 } else {
103 if(img.width != size[0] || img.height != size[1]) {
104 fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n",
105 (int)slices.size(), line, img.width, img.height, size[0], size[1]);
106 slices.clear();
107 return false;
108 }
109 }
110
111 slices.push_back(img);
112 img.pixels = 0; // otherwise the destructor will free it
113 }
114
115 size[2] = slices.size();
116 fclose(fp);
117 return true;
118 }
119
120 int VoxelVolume::num_samples(int dim) const
121 {
122 return size[dim];
123 }
124
125 static inline int clamp(int x, int low, int high)
126 {
127 return x < low ? low : (x > high ? high : x);
128 }
129
130 static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty,
131 float *pixels, int width, int height)
132 {
133 float v00, v01, v10, v11;
134
135 v00 = pixels[y0 * width + x0];
136 v01 = pixels[y1 * width + x0];
137 v10 = pixels[y0 * width + x1];
138 v11 = pixels[y1 * width + x1];
139
140 float v0 = v00 + (v01 - v00) * ty;
141 float v1 = v10 + (v11 - v10) * ty;
142
143 return v0 + (v1 - v0) * tx;
144 }
145
146 float VoxelVolume::valuef(float x, float y, float z) const
147 {
148 if(slices.empty()) {
149 return 0.0f;
150 }
151
152 float floor_x = floor(x);
153 float ceil_x = ceil(x);
154 float tx = x - floor_x;
155
156 float floor_y = floor(y);
157 float ceil_y = ceil(y);
158 float ty = y - floor_y;
159
160 float floor_z = floor(z);
161 float ceil_z = ceil(z);
162 float tz = z - floor_z;
163
164 int x0 = clamp(floor_x, 0, size[0] - 1);
165 int x1 = clamp(ceil_x, 0, size[0] - 1);
166
167 int y0 = clamp(floor_y, 0, size[1] - 1);
168 int y1 = clamp(ceil_y, 0, size[1] - 1);
169
170 int s0 = clamp(floor_z, 0, size[2] - 1);
171 int s1 = clamp(ceil_z, 0, size[2] - 1);
172
173 float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]);
174 float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]);
175
176 return val_s0 + (val_s1 - val_s0) * tz;
177 }
178
179 float VoxelVolume::valuei(int x, int y, int z) const
180 {
181 x = clamp(x, 0, size[0] - 1);
182 y = clamp(y, 0, size[1] - 1);
183 z = clamp(z, 0, size[2] - 1);
184
185 return slices[z].pixels[y * size[0] + x];
186 }
187
188 static char *strip_space(char *s)
189 {
190 while(*s && isspace(*s)) s++;
191 if(!*s) return 0;
192
193 char *end = s + strlen(s) - 1;
194 while(end > s && isspace(*end)) end--;
195 end[1] = 0;
196
197 return s;
198 }