1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/volume.cc	Mon Apr 09 23:42:57 2012 +0300
     1.3 @@ -0,0 +1,212 @@
     1.4 +#include <stdio.h>
     1.5 +#include <stdlib.h>
     1.6 +#include <string.h>
     1.7 +#include <ctype.h>
     1.8 +#include <list>
     1.9 +#include <string>
    1.10 +
    1.11 +#ifndef __APPLE__
    1.12 +#include <GL/gl.h>
    1.13 +#else
    1.14 +#include <OpenGL/gl.h>
    1.15 +#endif
    1.16 +
    1.17 +#include <imago2.h>
    1.18 +#include "volume.h"
    1.19 +
    1.20 +static void calc_gradients(float *voxels, int xsz, int ysz, int zsz);
    1.21 +static struct slice *read_voldesc(const char *fname, struct volume *vol);
    1.22 +static char *trim(char *str);
    1.23 +
    1.24 +Volume::Volume()
    1.25 +{
    1.26 +	sz[0] = sz[1] = sz[2] = 0;
    1.27 +	zaspect = 1.0;
    1.28 +	tex = 0;
    1.29 +}
    1.30 +
    1.31 +Volume::~Volume()
    1.32 +{
    1.33 +	if(tex) {
    1.34 +		glDeleteTextures(1, &tex);
    1.35 +	}
    1.36 +}
    1.37 +
    1.38 +bool Volume::load(const char *fname)
    1.39 +{
    1.40 +	std::list<std::string> slist;
    1.41 +	float *voxels = 0, *vptr;
    1.42 +	struct img_pixmap img;
    1.43 +
    1.44 +	if(!(read_voldesc(fname, &slist))) {
    1.45 +		return false;
    1.46 +	}
    1.47 +
    1.48 +	/* load the first image to determine slice dimensions */
    1.49 +	img_init(&img);
    1.50 +	if(img_load(&img, slist.begin().name.c_str()) == -1) {
    1.51 +		fprintf(stderr, "failed to load volume slice: %s\n", slist->name);
    1.52 +		return false;
    1.53 +	}
    1.54 +
    1.55 +	sz[0] = img.width;
    1.56 +	sz[1] = img.height;
    1.57 +	printf("volume dimensions: %dx%dx%d\n", img.width, img.height, vol->sz[2]);
    1.58 +
    1.59 +	/* allocate the whole volume at once */
    1.60 +	voxels = new float[sz[0] * sz[1] * sz[2] * 4];
    1.61 +	vptr = voxels;
    1.62 +	img_destroy(&img);
    1.63 +
    1.64 +	/* put the volume data into the alpha component */
    1.65 +	for(auto slice : slist) {
    1.66 +		int x, y, xsz, ysz;
    1.67 +		float *pixels, *src;
    1.68 +
    1.69 +		if(!(pixels = img_load_pixels(slice.name.c_str(), &xsz, &ysz, IMG_FMT_GREYF))) {
    1.70 +			fprintf(stderr, "failed to load volume slice: %s\n", slice.name.c_str());
    1.71 +			delete [] voxels;
    1.72 +			return false;
    1.73 +		}
    1.74 +
    1.75 +		if(xsz != sz[0] || ysz != sz[1]) {
    1.76 +			fprintf(stderr, "inconsistent dimensions for slice %s: %dx%d (%dx%d expected)\n",
    1.77 +					slice.name.c_str(), xsz, ysz, sz[0], sz[1]);
    1.78 +			delete [] voxels;
    1.79 +			return false;
    1.80 +		}
    1.81 +
    1.82 +		src = pixels;
    1.83 +		for(y=0; y<ysz; y++) {
    1.84 +			for(x=0; x<xsz; x++) {
    1.85 +				vptr[0] = vptr[1] = vptr[2] = 0.0f;
    1.86 +				vptr[3] = *src++;
    1.87 +				vptr += 4;
    1.88 +			}
    1.89 +		}
    1.90 +		img_free_pixels(pixels);
    1.91 +	}
    1.92 +
    1.93 +	/* calculate gradients */
    1.94 +	calc_gradients(voxels, sz[0], sz[1], sz[2]);
    1.95 +
    1.96 +	/* create the volume texture */
    1.97 +	glGenTextures(1, &tex);
    1.98 +	glBindTexture(GL_TEXTURE_3D, tex);
    1.99 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
   1.100 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
   1.101 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
   1.102 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
   1.103 +	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
   1.104 +	glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA32F_ARB, sz[0], sz[1], sz[2], 0, GL_RGBA, GL_FLOAT, voxels);
   1.105 +
   1.106 +	/* ... and destroy our copy */
   1.107 +	delete [] voxels;
   1.108 +	return true;
   1.109 +}
   1.110 +
   1.111 +unsigned int Volume::get_texture() const
   1.112 +{
   1.113 +	return tex;
   1.114 +}
   1.115 +
   1.116 +void Volume::bind(int tex_unit) const
   1.117 +{
   1.118 +	glActiveTexture(GL_TEXTURE0 + tex_unit);
   1.119 +	glBindTexture(GL_TEXTURE_3D, tex);
   1.120 +}
   1.121 +
   1.122 +struct vec4 { float x, y, z, w; };
   1.123 +
   1.124 +static void calc_gradients(float *voxels, int xsz, int ysz, int zsz)
   1.125 +{
   1.126 +	int x, y, z, slice_pixels = xsz * ysz;
   1.127 +	vec4 *vptr = (vec4*)voxels;
   1.128 +
   1.129 +	for(z=0; z<zsz; z++) {
   1.130 +		for(y=0; y<ysz; y++) {
   1.131 +			for(x=0; x<xsz; x++) {
   1.132 +				float x0, x1, y0, y1, z0, z1;
   1.133 +
   1.134 +				x0 = x > 0 ? (vptr - 1)->w : vptr->w;
   1.135 +				y0 = y > 0 ? (vptr - xsz)->w : vptr->w;
   1.136 +				z0 = z > 0 ? (vptr - slice_pixels)->w : vptr->w;
   1.137 +				x1 = x < xsz - 1 ? (vptr + 1)->w : vptr->w;
   1.138 +				y1 = y < ysz - 1 ? (vptr + xsz)->w : vptr->w;
   1.139 +				z1 = z < zsz - 1 ? (vptr + slice_pixels)->w : vptr->w;
   1.140 +
   1.141 +				vptr->x = x1 - x0;
   1.142 +				vptr->y = y1 - y0;
   1.143 +				vptr->z = z1 - z0;
   1.144 +
   1.145 +				vptr++;
   1.146 +			}
   1.147 +		}
   1.148 +	}
   1.149 +}
   1.150 +
   1.151 +bool read_voldesc(const char *fname, std::list<std::string> *slist) const
   1.152 +{
   1.153 +	FILE *fp = 0;
   1.154 +	char buf[512], *slash, *prefix = "";
   1.155 +	int mode_slices = 0;
   1.156 +
   1.157 +	slist->clear();
   1.158 +
   1.159 +	if(!(fp = fopen(fname, "r"))) {
   1.160 +		fprintf(stderr, "failed to open volume descriptor: %s\n", fname);
   1.161 +		return false;
   1.162 +	}
   1.163 +	fgets(buf, sizeof buf, fp);
   1.164 +	if(strstr(buf, "VOLDESC") != buf) {
   1.165 +		fprintf(stderr, "invalid file format while trying to read volume descriptor: %s\n", fname);
   1.166 +		fclose(fp);
   1.167 +		return false;
   1.168 +	}
   1.169 +
   1.170 +	if((slash = strrchr(fname, '/'))) {
   1.171 +		int len = slash - fname + 1;
   1.172 +		prefix = alloca(len + 1);
   1.173 +		memcpy(prefix, fname, len);
   1.174 +		prefix[len] = 0;
   1.175 +	}
   1.176 +
   1.177 +	while(fgets(buf, sizeof buf, fp)) {
   1.178 +		char *line = trim(buf);
   1.179 +
   1.180 +		if(!*line || *line == '#')
   1.181 +			continue;
   1.182 +
   1.183 +		if(mode_slices) {
   1.184 +			/* we're in the part that contains filenames... append to the list */
   1.185 +			slist->push_back(std::string(prefix) + std::string(line));
   1.186 +			sz[2]++;
   1.187 +		} else {
   1.188 +			/* TODO we're in the options part... parse */
   1.189 +			if(strcmp(line, "SLICES") == 0) {
   1.190 +				mode_slices = 1;
   1.191 +			}
   1.192 +		}
   1.193 +	}
   1.194 +	fclose(fp);
   1.195 +
   1.196 +	if(slist->empty()) {
   1.197 +		fprintf(stderr, "file: %s contains no slices\n", fname);
   1.198 +		return false;
   1.199 +	}
   1.200 +	return true;
   1.201 +}
   1.202 +
   1.203 +static char *trim(char *str)
   1.204 +{
   1.205 +	char *tmp = str + strlen(str) - 1;
   1.206 +
   1.207 +	while(isspace(*tmp))
   1.208 +		tmp--;
   1.209 +	tmp[1] = 0;
   1.210 +
   1.211 +	tmp = str;
   1.212 +	while(isspace(*tmp))
   1.213 +		tmp++;
   1.214 +	return tmp;
   1.215 +}