glviewvol

annotate src/volume.cc @ 15:2a67ea257ac0

find changesets by author, revision, files, or words in the commit message
makefile fixes for macosx
author John Tsiombikas <nuclear@member.fsf.org>
date Fri, 06 Feb 2015 23:47:28 +0200
parents 773f89037a35
children
rev   line source
nuclear@12 1 /*
nuclear@12 2 glviewvol is an OpenGL 3D volume data viewer
nuclear@12 3 Copyright (C) 2014 John Tsiombikas <nuclear@member.fsf.org>
nuclear@12 4
nuclear@12 5 This program is free software: you can redistribute it and/or modify
nuclear@12 6 it under the terms of the GNU General Public License as published by
nuclear@12 7 the Free Software Foundation, either version 3 of the License, or
nuclear@12 8 (at your option) any later version.
nuclear@12 9
nuclear@12 10 This program is distributed in the hope that it will be useful,
nuclear@12 11 but WITHOUT ANY WARRANTY; without even the implied warranty of
nuclear@12 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
nuclear@12 13 GNU General Public License for more details.
nuclear@12 14
nuclear@12 15 You should have received a copy of the GNU General Public License
nuclear@12 16 along with this program. If not, see <http://www.gnu.org/licenses/>.
nuclear@12 17 */
nuclear@1 18 #include <stdio.h>
nuclear@1 19 #include <string.h>
nuclear@1 20 #include <ctype.h>
nuclear@2 21 #include <alloca.h>
nuclear@1 22 #include <errno.h>
nuclear@1 23 #include <math.h>
nuclear@1 24 #include "volume.h"
nuclear@1 25
nuclear@1 26 static char *strip_space(char *s);
nuclear@1 27
nuclear@1 28
nuclear@1 29 Volume::~Volume()
nuclear@1 30 {
nuclear@1 31 }
nuclear@1 32
nuclear@1 33 int Volume::num_samples(int dim) const
nuclear@1 34 {
nuclear@1 35 return 0;
nuclear@1 36 }
nuclear@1 37
nuclear@1 38 void Volume::normalf(float *norm, float x, float y, float z, float delta)
nuclear@1 39 {
nuclear@1 40 float dx, dy, dz;
nuclear@1 41 dx = dy = dz = delta;
nuclear@1 42
nuclear@1 43 if(num_samples(0) > 0) {
nuclear@1 44 // discrete volume
nuclear@1 45 dx /= (float)num_samples(0);
nuclear@1 46 dy /= (float)num_samples(1);
nuclear@1 47 dz /= (float)num_samples(2);
nuclear@1 48 }
nuclear@1 49
nuclear@1 50 norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z);
nuclear@1 51 norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z);
nuclear@1 52 norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz);
nuclear@1 53
nuclear@1 54 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
nuclear@1 55 if(len != 0.0) {
nuclear@1 56 norm[0] /= len;
nuclear@1 57 norm[1] /= len;
nuclear@1 58 norm[2] /= len;
nuclear@1 59 }
nuclear@1 60 }
nuclear@1 61
nuclear@1 62 void Volume::normali(float *norm, int x, int y, int z)
nuclear@1 63 {
nuclear@1 64 int sz[3];
nuclear@1 65
nuclear@1 66 if((sz[0] = num_samples(0)) <= 0) {
nuclear@1 67 // if it's a continuous volume, just call normalf
nuclear@1 68 normalf(norm, x, y, z);
nuclear@1 69 return;
nuclear@1 70 }
nuclear@1 71 sz[1] = num_samples(1);
nuclear@1 72 sz[2] = num_samples(2);
nuclear@1 73
nuclear@1 74 int prevx = x <= 0 ? 0 : x;
nuclear@1 75 int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x;
nuclear@1 76 int prevy = y <= 0 ? 0 : y;
nuclear@1 77 int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y;
nuclear@1 78 int prevz = z <= 0 ? 0 : z;
nuclear@1 79 int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z;
nuclear@1 80
nuclear@1 81 norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z);
nuclear@1 82 norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z);
nuclear@1 83 norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz);
nuclear@1 84
nuclear@1 85 float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
nuclear@1 86 if(len != 0.0) {
nuclear@1 87 norm[0] /= len;
nuclear@1 88 norm[1] /= len;
nuclear@1 89 norm[2] /= len;
nuclear@1 90 }
nuclear@1 91 }
nuclear@1 92
nuclear@1 93 // ---- VoxelVolume (discrete) implementation ----
nuclear@1 94 VoxelVolume::VoxelVolume()
nuclear@1 95 {
nuclear@1 96 size[0] = size[1] = size[2] = 0;
nuclear@14 97
nuclear@14 98 hist_valid = false;
nuclear@14 99 hist = 0;
nuclear@14 100 num_hist_samples = 0;
nuclear@1 101 }
nuclear@1 102
nuclear@3 103 VoxelVolume::~VoxelVolume()
nuclear@3 104 {
nuclear@3 105 for(size_t i=0; i<slices.size(); i++) {
nuclear@3 106 slices[i].destroy();
nuclear@3 107 }
nuclear@3 108 }
nuclear@3 109
nuclear@1 110 bool VoxelVolume::load(const char *fname)
nuclear@1 111 {
nuclear@3 112 if(!fname) return false;
nuclear@3 113
nuclear@2 114 char *prefix = (char*)alloca(strlen(fname) + 1);
nuclear@2 115 strcpy(prefix, fname);
nuclear@2 116 char *slash = strrchr(prefix, '/');
nuclear@2 117 if(slash) {
nuclear@2 118 *slash = 0;
nuclear@2 119 } else {
nuclear@2 120 prefix = 0;
nuclear@2 121 }
nuclear@2 122
nuclear@4 123 printf("loading volume dataset: %s\n", fname);
nuclear@1 124 FILE *fp = fopen(fname, "r");
nuclear@1 125 if(!fp) {
nuclear@1 126 fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno));
nuclear@1 127 return false;
nuclear@1 128 }
nuclear@1 129
nuclear@2 130 char buf[256], path[300];
nuclear@1 131 while(fgets(buf, sizeof buf, fp)) {
nuclear@1 132 char *line = strip_space(buf);
nuclear@2 133 sprintf(path, "%s/%s", prefix, line);
nuclear@1 134
nuclear@4 135 printf(" loading slice %d: %s\n", (int)slices.size(), path);
nuclear@1 136 Image img;
nuclear@2 137 if(!img.load(path)) {
nuclear@1 138 slices.clear();
nuclear@1 139 return false;
nuclear@1 140 }
nuclear@1 141
nuclear@1 142 if(slices.empty()) {
nuclear@1 143 size[0] = img.width;
nuclear@1 144 size[1] = img.height;
nuclear@1 145 } else {
nuclear@1 146 if(img.width != size[0] || img.height != size[1]) {
nuclear@1 147 fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n",
nuclear@1 148 (int)slices.size(), line, img.width, img.height, size[0], size[1]);
nuclear@1 149 slices.clear();
nuclear@1 150 return false;
nuclear@1 151 }
nuclear@1 152 }
nuclear@1 153
nuclear@1 154 slices.push_back(img);
nuclear@1 155 }
nuclear@1 156
nuclear@1 157 size[2] = slices.size();
nuclear@14 158 hist_valid = false;
nuclear@1 159 fclose(fp);
nuclear@1 160 return true;
nuclear@1 161 }
nuclear@1 162
nuclear@1 163 int VoxelVolume::num_samples(int dim) const
nuclear@1 164 {
nuclear@1 165 return size[dim];
nuclear@1 166 }
nuclear@1 167
nuclear@1 168 static inline int clamp(int x, int low, int high)
nuclear@1 169 {
nuclear@1 170 return x < low ? low : (x > high ? high : x);
nuclear@1 171 }
nuclear@1 172
nuclear@1 173 static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty,
nuclear@1 174 float *pixels, int width, int height)
nuclear@1 175 {
nuclear@1 176 float v00, v01, v10, v11;
nuclear@1 177
nuclear@1 178 v00 = pixels[y0 * width + x0];
nuclear@1 179 v01 = pixels[y1 * width + x0];
nuclear@1 180 v10 = pixels[y0 * width + x1];
nuclear@1 181 v11 = pixels[y1 * width + x1];
nuclear@1 182
nuclear@1 183 float v0 = v00 + (v01 - v00) * ty;
nuclear@1 184 float v1 = v10 + (v11 - v10) * ty;
nuclear@1 185
nuclear@1 186 return v0 + (v1 - v0) * tx;
nuclear@1 187 }
nuclear@1 188
nuclear@1 189 float VoxelVolume::valuef(float x, float y, float z) const
nuclear@1 190 {
nuclear@1 191 if(slices.empty()) {
nuclear@1 192 return 0.0f;
nuclear@1 193 }
nuclear@1 194
nuclear@1 195 float floor_x = floor(x);
nuclear@1 196 float ceil_x = ceil(x);
nuclear@1 197 float tx = x - floor_x;
nuclear@1 198
nuclear@1 199 float floor_y = floor(y);
nuclear@1 200 float ceil_y = ceil(y);
nuclear@1 201 float ty = y - floor_y;
nuclear@1 202
nuclear@1 203 float floor_z = floor(z);
nuclear@1 204 float ceil_z = ceil(z);
nuclear@1 205 float tz = z - floor_z;
nuclear@1 206
nuclear@1 207 int x0 = clamp(floor_x, 0, size[0] - 1);
nuclear@1 208 int x1 = clamp(ceil_x, 0, size[0] - 1);
nuclear@1 209
nuclear@1 210 int y0 = clamp(floor_y, 0, size[1] - 1);
nuclear@1 211 int y1 = clamp(ceil_y, 0, size[1] - 1);
nuclear@1 212
nuclear@1 213 int s0 = clamp(floor_z, 0, size[2] - 1);
nuclear@1 214 int s1 = clamp(ceil_z, 0, size[2] - 1);
nuclear@1 215
nuclear@1 216 float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]);
nuclear@1 217 float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]);
nuclear@1 218
nuclear@1 219 return val_s0 + (val_s1 - val_s0) * tz;
nuclear@1 220 }
nuclear@1 221
nuclear@1 222 float VoxelVolume::valuei(int x, int y, int z) const
nuclear@1 223 {
nuclear@1 224 x = clamp(x, 0, size[0] - 1);
nuclear@1 225 y = clamp(y, 0, size[1] - 1);
nuclear@1 226 z = clamp(z, 0, size[2] - 1);
nuclear@1 227
nuclear@1 228 return slices[z].pixels[y * size[0] + x];
nuclear@1 229 }
nuclear@1 230
nuclear@14 231 float *VoxelVolume::calc_histogram(int hist_samples)
nuclear@14 232 {
nuclear@14 233 if(hist && hist_samples == num_hist_samples && hist_valid) {
nuclear@14 234 return hist;
nuclear@14 235 }
nuclear@14 236
nuclear@14 237 int num_pixels = size[0] * size[1] * size[2];
nuclear@14 238 if(!num_pixels) {
nuclear@14 239 return 0;
nuclear@14 240 }
nuclear@14 241
nuclear@14 242 delete [] hist;
nuclear@14 243 hist = new float[hist_samples];
nuclear@14 244 memset(hist, 0, hist_samples * sizeof *hist);
nuclear@14 245
nuclear@14 246 for(int i=0; i<size[2]; i++) {
nuclear@14 247 float *pptr = slices[i].pixels;
nuclear@14 248 for(int j=0; j<size[0] * size[1]; j++) {
nuclear@14 249 int idx = (int)(*pptr++ * (float)hist_samples);
nuclear@14 250
nuclear@14 251 hist[idx] += 1.0;
nuclear@14 252 }
nuclear@14 253 }
nuclear@14 254
nuclear@14 255 for(int i=0; i<hist_samples; i++) {
nuclear@14 256 hist[i] /= (float)num_pixels;
nuclear@14 257 }
nuclear@14 258
nuclear@14 259 hist_valid = true;
nuclear@14 260 num_hist_samples = hist_samples;
nuclear@14 261 return hist;
nuclear@14 262 }
nuclear@14 263
nuclear@1 264 static char *strip_space(char *s)
nuclear@1 265 {
nuclear@1 266 while(*s && isspace(*s)) s++;
nuclear@1 267 if(!*s) return 0;
nuclear@1 268
nuclear@1 269 char *end = s + strlen(s) - 1;
nuclear@1 270 while(end > s && isspace(*end)) end--;
nuclear@1 271 end[1] = 0;
nuclear@1 272
nuclear@1 273 return s;
nuclear@1 274 }