1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/volume.cc	Sat Dec 27 06:32:28 2014 +0200
     1.3 @@ -0,0 +1,198 @@
     1.4 +#include <stdio.h>
     1.5 +#include <string.h>
     1.6 +#include <ctype.h>
     1.7 +#include <errno.h>
     1.8 +#include <math.h>
     1.9 +#include "volume.h"
    1.10 +
    1.11 +static char *strip_space(char *s);
    1.12 +
    1.13 +
    1.14 +Volume::~Volume()
    1.15 +{
    1.16 +}
    1.17 +
    1.18 +int Volume::num_samples(int dim) const
    1.19 +{
    1.20 +	return 0;
    1.21 +}
    1.22 +
    1.23 +void Volume::normalf(float *norm, float x, float y, float z, float delta)
    1.24 +{
    1.25 +	float dx, dy, dz;
    1.26 +	dx = dy = dz = delta;
    1.27 +
    1.28 +	if(num_samples(0) > 0) {
    1.29 +		// discrete volume
    1.30 +		dx /= (float)num_samples(0);
    1.31 +		dy /= (float)num_samples(1);
    1.32 +		dz /= (float)num_samples(2);
    1.33 +	}
    1.34 +
    1.35 +	norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z);
    1.36 +	norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z);
    1.37 +	norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz);
    1.38 +
    1.39 +	float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
    1.40 +	if(len != 0.0) {
    1.41 +		norm[0] /= len;
    1.42 +		norm[1] /= len;
    1.43 +		norm[2] /= len;
    1.44 +	}
    1.45 +}
    1.46 +
    1.47 +void Volume::normali(float *norm, int x, int y, int z)
    1.48 +{
    1.49 +	int sz[3];
    1.50 +
    1.51 +	if((sz[0] = num_samples(0)) <= 0) {
    1.52 +		// if it's a continuous volume, just call normalf
    1.53 +		normalf(norm, x, y, z);
    1.54 +		return;
    1.55 +	}
    1.56 +	sz[1] = num_samples(1);
    1.57 +	sz[2] = num_samples(2);
    1.58 +
    1.59 +	int prevx = x <= 0 ? 0 : x;
    1.60 +	int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x;
    1.61 +	int prevy = y <= 0 ? 0 : y;
    1.62 +	int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y;
    1.63 +	int prevz = z <= 0 ? 0 : z;
    1.64 +	int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z;
    1.65 +
    1.66 +	norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z);
    1.67 +	norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z);
    1.68 +	norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz);
    1.69 +
    1.70 +	float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
    1.71 +	if(len != 0.0) {
    1.72 +		norm[0] /= len;
    1.73 +		norm[1] /= len;
    1.74 +		norm[2] /= len;
    1.75 +	}
    1.76 +}
    1.77 +
    1.78 +// ---- VoxelVolume (discrete) implementation ----
    1.79 +VoxelVolume::VoxelVolume()
    1.80 +{
    1.81 +	size[0] = size[1] = size[2] = 0;
    1.82 +}
    1.83 +
    1.84 +bool VoxelVolume::load(const char *fname)
    1.85 +{
    1.86 +	FILE *fp = fopen(fname, "r");
    1.87 +	if(!fp) {
    1.88 +		fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno));
    1.89 +		return false;
    1.90 +	}
    1.91 +
    1.92 +	char buf[512];
    1.93 +	while(fgets(buf, sizeof buf, fp)) {
    1.94 +		char *line = strip_space(buf);
    1.95 +
    1.96 +		Image img;
    1.97 +		if(!img.load(line)) {
    1.98 +			slices.clear();
    1.99 +			return false;
   1.100 +		}
   1.101 +
   1.102 +		if(slices.empty()) {
   1.103 +			size[0] = img.width;
   1.104 +			size[1] = img.height;
   1.105 +		} else {
   1.106 +			if(img.width != size[0] || img.height != size[1]) {
   1.107 +				fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n",
   1.108 +						(int)slices.size(), line, img.width, img.height, size[0], size[1]);
   1.109 +				slices.clear();
   1.110 +				return false;
   1.111 +			}
   1.112 +		}
   1.113 +
   1.114 +		slices.push_back(img);
   1.115 +		img.pixels = 0;	// otherwise the destructor will free it
   1.116 +	}
   1.117 +
   1.118 +	size[2] = slices.size();
   1.119 +	fclose(fp);
   1.120 +	return true;
   1.121 +}
   1.122 +
   1.123 +int VoxelVolume::num_samples(int dim) const
   1.124 +{
   1.125 +	return size[dim];
   1.126 +}
   1.127 +
   1.128 +static inline int clamp(int x, int low, int high)
   1.129 +{
   1.130 +	return x < low ? low : (x > high ? high : x);
   1.131 +}
   1.132 +
   1.133 +static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty,
   1.134 +		float *pixels, int width, int height)
   1.135 +{
   1.136 +	float v00, v01, v10, v11;
   1.137 +
   1.138 +	v00 = pixels[y0 * width + x0];
   1.139 +	v01 = pixels[y1 * width + x0];
   1.140 +	v10 = pixels[y0 * width + x1];
   1.141 +	v11 = pixels[y1 * width + x1];
   1.142 +
   1.143 +	float v0 = v00 + (v01 - v00) * ty;
   1.144 +	float v1 = v10 + (v11 - v10) * ty;
   1.145 +
   1.146 +	return v0 + (v1 - v0) * tx;
   1.147 +}
   1.148 +
   1.149 +float VoxelVolume::valuef(float x, float y, float z) const
   1.150 +{
   1.151 +	if(slices.empty()) {
   1.152 +		return 0.0f;
   1.153 +	}
   1.154 +
   1.155 +	float floor_x = floor(x);
   1.156 +	float ceil_x = ceil(x);
   1.157 +	float tx = x - floor_x;
   1.158 +
   1.159 +	float floor_y = floor(y);
   1.160 +	float ceil_y = ceil(y);
   1.161 +	float ty = y - floor_y;
   1.162 +
   1.163 +	float floor_z = floor(z);
   1.164 +	float ceil_z = ceil(z);
   1.165 +	float tz = z - floor_z;
   1.166 +
   1.167 +	int x0 = clamp(floor_x, 0, size[0] - 1);
   1.168 +	int x1 = clamp(ceil_x, 0, size[0] - 1);
   1.169 +
   1.170 +	int y0 = clamp(floor_y, 0, size[1] - 1);
   1.171 +	int y1 = clamp(ceil_y, 0, size[1] - 1);
   1.172 +
   1.173 +	int s0 = clamp(floor_z, 0, size[2] - 1);
   1.174 +	int s1 = clamp(ceil_z, 0, size[2] - 1);
   1.175 +
   1.176 +	float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]);
   1.177 +	float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]);
   1.178 +
   1.179 +	return val_s0 + (val_s1 - val_s0) * tz;
   1.180 +}
   1.181 +
   1.182 +float VoxelVolume::valuei(int x, int y, int z) const
   1.183 +{
   1.184 +	x = clamp(x, 0, size[0] - 1);
   1.185 +	y = clamp(y, 0, size[1] - 1);
   1.186 +	z = clamp(z, 0, size[2] - 1);
   1.187 +
   1.188 +	return slices[z].pixels[y * size[0] + x];
   1.189 +}
   1.190 +
   1.191 +static char *strip_space(char *s)
   1.192 +{
   1.193 +	while(*s && isspace(*s)) s++;
   1.194 +	if(!*s) return 0;
   1.195 +
   1.196 +	char *end = s + strlen(s) - 1;
   1.197 +	while(end > s && isspace(*end)) end--;
   1.198 +	end[1] = 0;
   1.199 +
   1.200 +	return s;
   1.201 +}