volmetrics

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