volmetrics

annotate src/volume.cc @ 34:eb87d4a12bd3

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