nuclear@12: /*
nuclear@12: glviewvol is an OpenGL 3D volume data viewer
nuclear@12: Copyright (C) 2014  John Tsiombikas <nuclear@member.fsf.org>
nuclear@12: 
nuclear@12: This program is free software: you can redistribute it and/or modify
nuclear@12: it under the terms of the GNU General Public License as published by
nuclear@12: the Free Software Foundation, either version 3 of the License, or
nuclear@12: (at your option) any later version.
nuclear@12: 
nuclear@12: This program is distributed in the hope that it will be useful,
nuclear@12: but WITHOUT ANY WARRANTY; without even the implied warranty of
nuclear@12: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
nuclear@12: GNU General Public License for more details.
nuclear@12: 
nuclear@12: You should have received a copy of the GNU General Public License
nuclear@12: along with this program.  If not, see <http://www.gnu.org/licenses/>.
nuclear@12: */
nuclear@1: #include <stdio.h>
nuclear@1: #include <string.h>
nuclear@1: #include <ctype.h>
nuclear@2: #include <alloca.h>
nuclear@1: #include <errno.h>
nuclear@1: #include <math.h>
nuclear@1: #include "volume.h"
nuclear@1: 
nuclear@1: static char *strip_space(char *s);
nuclear@1: 
nuclear@1: 
nuclear@1: Volume::~Volume()
nuclear@1: {
nuclear@1: }
nuclear@1: 
nuclear@1: int Volume::num_samples(int dim) const
nuclear@1: {
nuclear@1: 	return 0;
nuclear@1: }
nuclear@1: 
nuclear@1: void Volume::normalf(float *norm, float x, float y, float z, float delta)
nuclear@1: {
nuclear@1: 	float dx, dy, dz;
nuclear@1: 	dx = dy = dz = delta;
nuclear@1: 
nuclear@1: 	if(num_samples(0) > 0) {
nuclear@1: 		// discrete volume
nuclear@1: 		dx /= (float)num_samples(0);
nuclear@1: 		dy /= (float)num_samples(1);
nuclear@1: 		dz /= (float)num_samples(2);
nuclear@1: 	}
nuclear@1: 
nuclear@1: 	norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z);
nuclear@1: 	norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z);
nuclear@1: 	norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz);
nuclear@1: 
nuclear@1: 	float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
nuclear@1: 	if(len != 0.0) {
nuclear@1: 		norm[0] /= len;
nuclear@1: 		norm[1] /= len;
nuclear@1: 		norm[2] /= len;
nuclear@1: 	}
nuclear@1: }
nuclear@1: 
nuclear@1: void Volume::normali(float *norm, int x, int y, int z)
nuclear@1: {
nuclear@1: 	int sz[3];
nuclear@1: 
nuclear@1: 	if((sz[0] = num_samples(0)) <= 0) {
nuclear@1: 		// if it's a continuous volume, just call normalf
nuclear@1: 		normalf(norm, x, y, z);
nuclear@1: 		return;
nuclear@1: 	}
nuclear@1: 	sz[1] = num_samples(1);
nuclear@1: 	sz[2] = num_samples(2);
nuclear@1: 
nuclear@1: 	int prevx = x <= 0 ? 0 : x;
nuclear@1: 	int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x;
nuclear@1: 	int prevy = y <= 0 ? 0 : y;
nuclear@1: 	int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y;
nuclear@1: 	int prevz = z <= 0 ? 0 : z;
nuclear@1: 	int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z;
nuclear@1: 
nuclear@1: 	norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z);
nuclear@1: 	norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z);
nuclear@1: 	norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz);
nuclear@1: 
nuclear@1: 	float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
nuclear@1: 	if(len != 0.0) {
nuclear@1: 		norm[0] /= len;
nuclear@1: 		norm[1] /= len;
nuclear@1: 		norm[2] /= len;
nuclear@1: 	}
nuclear@1: }
nuclear@1: 
nuclear@1: // ---- VoxelVolume (discrete) implementation ----
nuclear@1: VoxelVolume::VoxelVolume()
nuclear@1: {
nuclear@1: 	size[0] = size[1] = size[2] = 0;
nuclear@14: 
nuclear@14: 	hist_valid = false;
nuclear@14: 	hist = 0;
nuclear@14: 	num_hist_samples = 0;
nuclear@1: }
nuclear@1: 
nuclear@3: VoxelVolume::~VoxelVolume()
nuclear@3: {
nuclear@3: 	for(size_t i=0; i<slices.size(); i++) {
nuclear@3: 		slices[i].destroy();
nuclear@3: 	}
nuclear@3: }
nuclear@3: 
nuclear@1: bool VoxelVolume::load(const char *fname)
nuclear@1: {
nuclear@3: 	if(!fname) return false;
nuclear@3: 
nuclear@2: 	char *prefix = (char*)alloca(strlen(fname) + 1);
nuclear@2: 	strcpy(prefix, fname);
nuclear@2: 	char *slash = strrchr(prefix, '/');
nuclear@2: 	if(slash) {
nuclear@2: 		*slash = 0;
nuclear@2: 	} else {
nuclear@2: 		prefix = 0;
nuclear@2: 	}
nuclear@2: 
nuclear@4: 	printf("loading volume dataset: %s\n", fname);
nuclear@1: 	FILE *fp = fopen(fname, "r");
nuclear@1: 	if(!fp) {
nuclear@1: 		fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno));
nuclear@1: 		return false;
nuclear@1: 	}
nuclear@1: 
nuclear@2: 	char buf[256], path[300];
nuclear@1: 	while(fgets(buf, sizeof buf, fp)) {
nuclear@1: 		char *line = strip_space(buf);
nuclear@2: 		sprintf(path, "%s/%s", prefix, line);
nuclear@1: 
nuclear@4: 		printf("  loading slice %d: %s\n", (int)slices.size(), path);
nuclear@1: 		Image img;
nuclear@2: 		if(!img.load(path)) {
nuclear@1: 			slices.clear();
nuclear@1: 			return false;
nuclear@1: 		}
nuclear@1: 
nuclear@1: 		if(slices.empty()) {
nuclear@1: 			size[0] = img.width;
nuclear@1: 			size[1] = img.height;
nuclear@1: 		} else {
nuclear@1: 			if(img.width != size[0] || img.height != size[1]) {
nuclear@1: 				fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n",
nuclear@1: 						(int)slices.size(), line, img.width, img.height, size[0], size[1]);
nuclear@1: 				slices.clear();
nuclear@1: 				return false;
nuclear@1: 			}
nuclear@1: 		}
nuclear@1: 
nuclear@1: 		slices.push_back(img);
nuclear@1: 	}
nuclear@1: 
nuclear@1: 	size[2] = slices.size();
nuclear@14: 	hist_valid = false;
nuclear@1: 	fclose(fp);
nuclear@1: 	return true;
nuclear@1: }
nuclear@1: 
nuclear@1: int VoxelVolume::num_samples(int dim) const
nuclear@1: {
nuclear@1: 	return size[dim];
nuclear@1: }
nuclear@1: 
nuclear@1: static inline int clamp(int x, int low, int high)
nuclear@1: {
nuclear@1: 	return x < low ? low : (x > high ? high : x);
nuclear@1: }
nuclear@1: 
nuclear@1: static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty,
nuclear@1: 		float *pixels, int width, int height)
nuclear@1: {
nuclear@1: 	float v00, v01, v10, v11;
nuclear@1: 
nuclear@1: 	v00 = pixels[y0 * width + x0];
nuclear@1: 	v01 = pixels[y1 * width + x0];
nuclear@1: 	v10 = pixels[y0 * width + x1];
nuclear@1: 	v11 = pixels[y1 * width + x1];
nuclear@1: 
nuclear@1: 	float v0 = v00 + (v01 - v00) * ty;
nuclear@1: 	float v1 = v10 + (v11 - v10) * ty;
nuclear@1: 
nuclear@1: 	return v0 + (v1 - v0) * tx;
nuclear@1: }
nuclear@1: 
nuclear@1: float VoxelVolume::valuef(float x, float y, float z) const
nuclear@1: {
nuclear@1: 	if(slices.empty()) {
nuclear@1: 		return 0.0f;
nuclear@1: 	}
nuclear@1: 
nuclear@1: 	float floor_x = floor(x);
nuclear@1: 	float ceil_x = ceil(x);
nuclear@1: 	float tx = x - floor_x;
nuclear@1: 
nuclear@1: 	float floor_y = floor(y);
nuclear@1: 	float ceil_y = ceil(y);
nuclear@1: 	float ty = y - floor_y;
nuclear@1: 
nuclear@1: 	float floor_z = floor(z);
nuclear@1: 	float ceil_z = ceil(z);
nuclear@1: 	float tz = z - floor_z;
nuclear@1: 
nuclear@1: 	int x0 = clamp(floor_x, 0, size[0] - 1);
nuclear@1: 	int x1 = clamp(ceil_x, 0, size[0] - 1);
nuclear@1: 
nuclear@1: 	int y0 = clamp(floor_y, 0, size[1] - 1);
nuclear@1: 	int y1 = clamp(ceil_y, 0, size[1] - 1);
nuclear@1: 
nuclear@1: 	int s0 = clamp(floor_z, 0, size[2] - 1);
nuclear@1: 	int s1 = clamp(ceil_z, 0, size[2] - 1);
nuclear@1: 
nuclear@1: 	float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]);
nuclear@1: 	float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]);
nuclear@1: 
nuclear@1: 	return val_s0 + (val_s1 - val_s0) * tz;
nuclear@1: }
nuclear@1: 
nuclear@1: float VoxelVolume::valuei(int x, int y, int z) const
nuclear@1: {
nuclear@1: 	x = clamp(x, 0, size[0] - 1);
nuclear@1: 	y = clamp(y, 0, size[1] - 1);
nuclear@1: 	z = clamp(z, 0, size[2] - 1);
nuclear@1: 
nuclear@1: 	return slices[z].pixels[y * size[0] + x];
nuclear@1: }
nuclear@1: 
nuclear@14: float *VoxelVolume::calc_histogram(int hist_samples)
nuclear@14: {
nuclear@14: 	if(hist && hist_samples == num_hist_samples && hist_valid) {
nuclear@14: 		return hist;
nuclear@14: 	}
nuclear@14: 
nuclear@14: 	int num_pixels = size[0] * size[1] * size[2];
nuclear@14: 	if(!num_pixels) {
nuclear@14: 		return 0;
nuclear@14: 	}
nuclear@14: 
nuclear@14: 	delete [] hist;
nuclear@14: 	hist = new float[hist_samples];
nuclear@14: 	memset(hist, 0, hist_samples * sizeof *hist);
nuclear@14: 
nuclear@14: 	for(int i=0; i<size[2]; i++) {
nuclear@14: 		float *pptr = slices[i].pixels;
nuclear@14: 		for(int j=0; j<size[0] * size[1]; j++) {
nuclear@14: 			int idx = (int)(*pptr++ * (float)hist_samples);
nuclear@14: 
nuclear@14: 			hist[idx] += 1.0;
nuclear@14: 		}
nuclear@14: 	}
nuclear@14: 
nuclear@14: 	for(int i=0; i<hist_samples; i++) {
nuclear@14: 		hist[i] /= (float)num_pixels;
nuclear@14: 	}
nuclear@14: 
nuclear@14: 	hist_valid = true;
nuclear@14: 	num_hist_samples = hist_samples;
nuclear@14: 	return hist;
nuclear@14: }
nuclear@14: 
nuclear@1: static char *strip_space(char *s)
nuclear@1: {
nuclear@1: 	while(*s && isspace(*s)) s++;
nuclear@1: 	if(!*s) return 0;
nuclear@1: 
nuclear@1: 	char *end = s + strlen(s) - 1;
nuclear@1: 	while(end > s && isspace(*end)) end--;
nuclear@1: 	end[1] = 0;
nuclear@1: 
nuclear@1: 	return s;
nuclear@1: }