volmetrics

annotate src/volume.cc @ 37:1d935677d3ab

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