volmetrics

view src/volume.cc @ 16:add30e2d5253

find changesets by author, revision, files, or words in the commit message
foo
author Eleni Maria Stea <elene.mst@gmail.com>
date Mon, 03 Feb 2014 02:52:46 +0200
parents 77a376198d38
children 4e120dcd55ec
line source
1 #include <GL/glew.h>
2
3 #include <assert.h>
4
5 #include <ctype.h>
6 #include <errno.h>
7 #include <stdio.h>
8 #include <string.h>
9
10 #include <string>
11 #include <vector>
12
13 #include <metasurf.h>
14
15 #include "mesh.h"
16 #include "volume.h"
17
18 static char *strip_whitespaces(char *buf);
19
20 Volume::Volume()
21 {
22 width = height = 0;
23 zaspect = 1;
24 vol_tex_valid = false;
25
26 glGenTextures(1, &vol_tex);
27 }
28
29 Volume::~Volume()
30 {
31 glDeleteTextures(1, &vol_tex);
32 }
33
34 void Volume::create_vol_tex() const
35 {
36 if(slices.empty())
37 return;
38
39 glBindTexture(GL_TEXTURE_3D, vol_tex);
40 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
41 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
42 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
43 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
44 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
45
46 assert(glGetError() == GL_NO_ERROR);
47 glTexImage3D(GL_TEXTURE_3D, 0, GL_LUMINANCE, width, height,
48 slices.size(), 0, GL_LUMINANCE, GL_FLOAT, 0);
49 assert(glGetError() == GL_NO_ERROR);
50
51 for(size_t i=0; i<slices.size(); i++) {
52 glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, i, width, height,
53 1, GL_LUMINANCE, GL_FLOAT, slices[i].get_pixels());
54 }
55 }
56
57 bool Volume::load(const char *fname)
58 {
59 FILE *fp = fopen(fname, "r");
60 if(!fp) {
61 fprintf(stderr, "Failed to open file: %s: %s\n", fname, strerror(errno));
62 return false;
63 }
64
65 char buf[512];
66 if(!fgets(buf, sizeof buf, fp)) {
67 fprintf(stderr, "Empty file: %s.\n", fname);
68 return false;
69 }
70 if(strstr(buf, "VOLUME") != buf) {
71 fprintf(stderr, "Invalid volume file format: %s\n", fname);
72 return false;
73 }
74
75 bool reading_slices = false;
76 std::vector<std::string> img_fnames;
77
78 while(fgets(buf, sizeof buf, fp)) {
79 char *line = strip_whitespaces(buf);
80 if(!line || *line == 0 || *line == '#')
81 continue;
82
83 if(reading_slices == false) {
84 if(strcmp(line, "SLICES") == 0)
85 reading_slices = true;
86
87 sscanf(line, "zaspect = %f", &zaspect);
88 }
89 else {
90 img_fnames.push_back(line);
91 }
92 }
93
94 fclose(fp);
95
96 for(size_t i=0; i<img_fnames.size(); i++) {
97 Image img;
98 if(img.load(img_fnames[i].c_str())) {
99 push_slice(std::move(img));
100 }
101 else {
102 fprintf(stderr, "Failed to load slice: %s\n", img_fnames[i].c_str());
103 }
104 }
105
106 return true;
107 }
108
109 bool Volume::push_slice(Image &&slice)
110 {
111 if(slices.empty()) {
112 width = slice.get_width();
113 height = slice.get_height();
114 }
115 else {
116 if(width != slice.get_width() || height != slice.get_height()) {
117 fprintf(stderr, "failed to load slice no: %d\n", (int)slices.size() + 1);
118 return false;
119 }
120 }
121
122 slices.push_back(slice);
123 return true;
124 }
125
126 const Image *Volume::get_slice(int num_slice) const
127 {
128 if(num_slice >= (int)slices.size() || num_slice < 0)
129 return 0;
130
131 return &slices[num_slice];
132 }
133
134 const Image *Volume::get_slice_by_z(float z) const
135 {
136 int idx = get_slice_idx_by_z(z);
137 return get_slice(idx);
138 }
139
140 int Volume::get_slice_count() const
141 {
142 return (int)slices.size();
143 }
144
145 int Volume::get_slice_idx_by_z(float z) const
146 {
147 int last_slice_idx = (int)slices.size() - 1;
148 int idx = z * last_slice_idx;
149 if(idx < 0) {
150 idx = 0;
151 }
152 else if(idx > last_slice_idx) {
153 idx = last_slice_idx;
154 }
155
156 return idx;
157 }
158
159 unsigned int Volume::get_texture() const
160 {
161 if(!vol_tex_valid) {
162 create_vol_tex();
163 vol_tex_valid = true;
164 }
165
166 return vol_tex;
167 }
168
169 static Mesh *cur_mesh;
170 static Volume *cur_vol;
171 static float low_thres, high_thres;
172
173 static float smoothstep(float x, float low, float high)
174 {
175 if(x < low) return 0.0;
176 if(x >= high) return 1.0;
177
178 x = (x - low) / (high - low);
179 return x * x * (3.0 - 2.0 * x);
180 }
181
182 float transfer_function(float x, float low_thres, float high_thres)
183 {
184 float dt = 0.25 * (high_thres - low_thres);
185 return smoothstep(x, low_thres - dt, low_thres + dt) *
186 (1 - smoothstep(x, high_thres - dt, high_thres + dt));
187 }
188
189 static float cb_eval(float x, float y, float z)
190 {
191 const Image *img = cur_vol->get_slice_by_z((z + 1) / 2.0);
192 const float *pixels = img->get_pixels();
193
194 int px = (x + 1) / 2.0 * img->get_width();
195 int py = (y + 1) / 2.0 * img->get_height();
196
197 float val = pixels[px + img->get_width() * py];
198 return transfer_function(val, low_thres, high_thres);
199 }
200
201 static void cb_vertex(float x, float y, float z)
202 {
203 float dx = 1.0 / cur_vol->get_slice(0)->get_width();
204 float dy = 1.0 / cur_vol->get_slice(0)->get_height();
205 float dz = 1.0 / cur_vol->get_slice_count();
206 float dfdx = cb_eval(x - dx, y, z) - cb_eval(x + dx, y, z);
207 float dfdy = cb_eval(x, y - dy, z) - cb_eval(x, y + dy, z);
208 float dfdz = cb_eval(x, y, z - dz) - cb_eval(x, y, z + dz);
209
210 cur_mesh->add_normal(Vector3(dfdx, dfdy, dfdz));
211 cur_mesh->add_vertex(Vector3(x, y, z));
212 }
213
214 void Volume::create_mesh(Mesh *mesh, float tmin, float tmax, int xres, int yres, int zres)
215 {
216 if (tmin > tmax) {
217 float t = tmax;
218 tmax = tmin;
219 tmin = t;
220 }
221 low_thres = tmin; high_thres = tmax;
222
223 if(xres <= 0)
224 xres = width;
225 if(yres <= 0)
226 yres = height;
227 if(zres <= 0)
228 zres = slices.size();
229
230 cur_mesh = mesh;
231 cur_vol = this;
232
233 metasurface *ms = msurf_create();
234 msurf_threshold(ms, 0.5);
235 msurf_resolution(ms, xres, yres, zres);
236 msurf_vertex_func(ms, cb_vertex);
237 msurf_eval_func(ms, cb_eval);
238 msurf_polygonize(ms);
239 msurf_free(ms);
240 }
241
242 void Volume::draw() const
243 {
244 }
245
246 static char *strip_whitespaces(char *buf)
247 {
248 while(*buf && isspace(*buf))
249 buf++;
250
251 if(!*buf)
252 return 0;
253
254 char *end = buf + strlen(buf) - 1;
255 while(end > buf && isspace(*end))
256 end--;
257 end[1] = 0;
258 return buf;
259 }