volmetrics

annotate src/volume.cc @ 31:cfeecd093730

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