nuclear@8: #include <stdio.h>
nuclear@8: #include <stdlib.h>
nuclear@8: #include <string.h>
nuclear@8: #include <ctype.h>
nuclear@11: #include <list>
nuclear@11: #include <string>
nuclear@17: #include <GL/glew.h>
nuclear@8: #include <imago2.h>
nuclear@9: #include "volume.h"
nuclear@8: 
nuclear@9: static void calc_gradients(float *voxels, int xsz, int ysz, int zsz);
nuclear@8: static struct slice *read_voldesc(const char *fname, struct volume *vol);
nuclear@8: static char *trim(char *str);
nuclear@8: 
nuclear@11: Volume::Volume()
nuclear@8: {
nuclear@11: 	sz[0] = sz[1] = sz[2] = 0;
nuclear@11: 	zaspect = 1.0;
nuclear@11: 	tex = 0;
nuclear@11: }
nuclear@11: 
nuclear@11: Volume::~Volume()
nuclear@11: {
nuclear@11: 	if(tex) {
nuclear@11: 		glDeleteTextures(1, &tex);
nuclear@11: 	}
nuclear@11: }
nuclear@11: 
nuclear@18: bool Volume::create(int xsz, int ysz, int zsz, float *data)
nuclear@18: {
nuclear@18: 	if(!tex)
nuclear@18: 		glGenTextures(1, &tex);
nuclear@18: 	glBindTexture(GL_TEXTURE_3D, tex);
nuclear@18: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
nuclear@18: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
nuclear@18: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
nuclear@18: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
nuclear@18: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
nuclear@18: 	glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA32F_ARB, xsz, ysz, zsz, 0, GL_RGBA, GL_FLOAT, data);
nuclear@18: }
nuclear@18: 
nuclear@11: bool Volume::load(const char *fname)
nuclear@11: {
nuclear@11: 	std::list<std::string> slist;
nuclear@9: 	float *voxels = 0, *vptr;
nuclear@9: 	struct img_pixmap img;
nuclear@8: 
nuclear@11: 	if(!(read_voldesc(fname, &slist))) {
nuclear@11: 		return false;
nuclear@8: 	}
nuclear@8: 
nuclear@9: 	/* load the first image to determine slice dimensions */
nuclear@9: 	img_init(&img);
nuclear@13: 	if(img_load(&img, slist.begin()->c_str()) == -1) {
nuclear@13: 		fprintf(stderr, "failed to load volume slice: %s\n", slist.begin()->c_str());
nuclear@11: 		return false;
nuclear@9: 	}
nuclear@9: 
nuclear@11: 	sz[0] = img.width;
nuclear@11: 	sz[1] = img.height;
nuclear@13: 	printf("volume dimensions: %dx%dx%d\n", img.width, img.height, sz[2]);
nuclear@9: 
nuclear@9: 	/* allocate the whole volume at once */
nuclear@11: 	voxels = new float[sz[0] * sz[1] * sz[2] * 4];
nuclear@9: 	vptr = voxels;
nuclear@9: 	img_destroy(&img);
nuclear@9: 
nuclear@9: 	/* put the volume data into the alpha component */
nuclear@17: 	//for(auto slice : slist) { // fucking vs2010 doesn't support range-for
nuclear@17: 	for(auto iter = slist.begin(); iter != slist.end(); iter++) {
nuclear@17: 		auto slice = *iter;
nuclear@17: 
nuclear@9: 		int x, y, xsz, ysz;
nuclear@9: 		float *pixels, *src;
nuclear@9: 
nuclear@13: 		if(!(pixels = (float*)img_load_pixels(slice.c_str(), &xsz, &ysz, IMG_FMT_GREYF))) {
nuclear@13: 			fprintf(stderr, "failed to load volume slice: %s\n", slice.c_str());
nuclear@11: 			delete [] voxels;
nuclear@11: 			return false;
nuclear@9: 		}
nuclear@9: 
nuclear@11: 		if(xsz != sz[0] || ysz != sz[1]) {
nuclear@9: 			fprintf(stderr, "inconsistent dimensions for slice %s: %dx%d (%dx%d expected)\n",
nuclear@13: 					slice.c_str(), xsz, ysz, sz[0], sz[1]);
nuclear@11: 			delete [] voxels;
nuclear@11: 			return false;
nuclear@9: 		}
nuclear@9: 
nuclear@9: 		src = pixels;
nuclear@9: 		for(y=0; y<ysz; y++) {
nuclear@9: 			for(x=0; x<xsz; x++) {
nuclear@9: 				vptr[0] = vptr[1] = vptr[2] = 0.0f;
nuclear@9: 				vptr[3] = *src++;
nuclear@9: 				vptr += 4;
nuclear@9: 			}
nuclear@9: 		}
nuclear@9: 		img_free_pixels(pixels);
nuclear@9: 	}
nuclear@9: 
nuclear@9: 	/* calculate gradients */
nuclear@11: 	calc_gradients(voxels, sz[0], sz[1], sz[2]);
nuclear@9: 
nuclear@9: 	/* create the volume texture */
nuclear@18: 	create(sz[0], sz[1], sz[2], voxels);
nuclear@8: 
nuclear@9: 	/* ... and destroy our copy */
nuclear@11: 	delete [] voxels;
nuclear@11: 	return true;
nuclear@9: }
nuclear@9: 
nuclear@11: unsigned int Volume::get_texture() const
nuclear@9: {
nuclear@11: 	return tex;
nuclear@11: }
nuclear@11: 
nuclear@11: void Volume::bind(int tex_unit) const
nuclear@11: {
nuclear@11: 	glActiveTexture(GL_TEXTURE0 + tex_unit);
nuclear@11: 	glBindTexture(GL_TEXTURE_3D, tex);
nuclear@9: }
nuclear@9: 
nuclear@9: struct vec4 { float x, y, z, w; };
nuclear@9: 
nuclear@9: static void calc_gradients(float *voxels, int xsz, int ysz, int zsz)
nuclear@9: {
nuclear@9: 	int x, y, z, slice_pixels = xsz * ysz;
nuclear@11: 	vec4 *vptr = (vec4*)voxels;
nuclear@9: 
nuclear@9: 	for(z=0; z<zsz; z++) {
nuclear@9: 		for(y=0; y<ysz; y++) {
nuclear@9: 			for(x=0; x<xsz; x++) {
nuclear@9: 				float x0, x1, y0, y1, z0, z1;
nuclear@9: 
nuclear@9: 				x0 = x > 0 ? (vptr - 1)->w : vptr->w;
nuclear@9: 				y0 = y > 0 ? (vptr - xsz)->w : vptr->w;
nuclear@9: 				z0 = z > 0 ? (vptr - slice_pixels)->w : vptr->w;
nuclear@9: 				x1 = x < xsz - 1 ? (vptr + 1)->w : vptr->w;
nuclear@9: 				y1 = y < ysz - 1 ? (vptr + xsz)->w : vptr->w;
nuclear@9: 				z1 = z < zsz - 1 ? (vptr + slice_pixels)->w : vptr->w;
nuclear@9: 
nuclear@9: 				vptr->x = x1 - x0;
nuclear@9: 				vptr->y = y1 - y0;
nuclear@9: 				vptr->z = z1 - z0;
nuclear@9: 
nuclear@9: 				vptr++;
nuclear@9: 			}
nuclear@9: 		}
nuclear@9: 	}
nuclear@8: }
nuclear@8: 
nuclear@13: bool Volume::read_voldesc(const char *fname, std::list<std::string> *slist)
nuclear@8: {
nuclear@8: 	FILE *fp = 0;
nuclear@9: 	char buf[512], *slash, *prefix = "";
nuclear@8: 	int mode_slices = 0;
nuclear@11: 
nuclear@11: 	slist->clear();
nuclear@8: 
nuclear@8: 	if(!(fp = fopen(fname, "r"))) {
nuclear@8: 		fprintf(stderr, "failed to open volume descriptor: %s\n", fname);
nuclear@11: 		return false;
nuclear@8: 	}
nuclear@8: 	fgets(buf, sizeof buf, fp);
nuclear@8: 	if(strstr(buf, "VOLDESC") != buf) {
nuclear@8: 		fprintf(stderr, "invalid file format while trying to read volume descriptor: %s\n", fname);
nuclear@11: 		fclose(fp);
nuclear@11: 		return false;
nuclear@8: 	}
nuclear@8: 
nuclear@13: 	if((slash = strrchr((char*)fname, '/'))) {
nuclear@9: 		int len = slash - fname + 1;
nuclear@13: 		prefix = (char*)alloca(len + 1);
nuclear@9: 		memcpy(prefix, fname, len);
nuclear@9: 		prefix[len] = 0;
nuclear@9: 	}
nuclear@9: 
nuclear@8: 	while(fgets(buf, sizeof buf, fp)) {
nuclear@8: 		char *line = trim(buf);
nuclear@8: 
nuclear@8: 		if(!*line || *line == '#')
nuclear@8: 			continue;
nuclear@8: 
nuclear@8: 		if(mode_slices) {
nuclear@8: 			/* we're in the part that contains filenames... append to the list */
nuclear@11: 			slist->push_back(std::string(prefix) + std::string(line));
nuclear@11: 			sz[2]++;
nuclear@8: 		} else {
nuclear@8: 			/* TODO we're in the options part... parse */
nuclear@9: 			if(strcmp(line, "SLICES") == 0) {
nuclear@9: 				mode_slices = 1;
nuclear@9: 			}
nuclear@8: 		}
nuclear@8: 	}
nuclear@8: 	fclose(fp);
nuclear@8: 
nuclear@11: 	if(slist->empty()) {
nuclear@8: 		fprintf(stderr, "file: %s contains no slices\n", fname);
nuclear@11: 		return false;
nuclear@8: 	}
nuclear@11: 	return true;
nuclear@8: }
nuclear@8: 
nuclear@8: static char *trim(char *str)
nuclear@8: {
nuclear@8: 	char *tmp = str + strlen(str) - 1;
nuclear@8: 
nuclear@8: 	while(isspace(*tmp))
nuclear@8: 		tmp--;
nuclear@8: 	tmp[1] = 0;
nuclear@8: 
nuclear@8: 	tmp = str;
nuclear@8: 	while(isspace(*tmp))
nuclear@8: 		tmp++;
nuclear@8: 	return tmp;
nuclear@8: }