1.1 --- a/src/volume.c	Sat Apr 07 16:07:12 2012 +0300
     1.2 +++ b/src/volume.c	Sun Apr 08 07:11:54 2012 +0300
     1.3 @@ -10,19 +10,24 @@
     1.4  #endif
     1.5  
     1.6  #include <imago2.h>
     1.7 +#include "volume.h"
     1.8  
     1.9  struct slice {
    1.10  	char *name;
    1.11  	struct slice *next;
    1.12  };
    1.13  
    1.14 +static void calc_gradients(float *voxels, int xsz, int ysz, int zsz);
    1.15  static struct slice *read_voldesc(const char *fname, struct volume *vol);
    1.16  static char *trim(char *str);
    1.17  
    1.18  struct volume *load_volume(const char *fname)
    1.19  {
    1.20 +	int i;
    1.21  	struct slice *slist;
    1.22  	struct volume *vol = 0;
    1.23 +	float *voxels = 0, *vptr;
    1.24 +	struct img_pixmap img;
    1.25  
    1.26  	if(!(vol = malloc(sizeof *vol))) {
    1.27  		perror("failed to allocate volume");
    1.28 @@ -35,6 +40,60 @@
    1.29  		goto err;
    1.30  	}
    1.31  
    1.32 +	/* load the first image to determine slice dimensions */
    1.33 +	img_init(&img);
    1.34 +	if(img_load(&img, slist->name) == -1) {
    1.35 +		fprintf(stderr, "failed to load volume slice: %s\n", slist->name);
    1.36 +		goto err;
    1.37 +	}
    1.38 +
    1.39 +	vol->sz[0] = img.width;
    1.40 +	vol->sz[1] = img.height;
    1.41 +	printf("volume dimensions: %dx%dx%d\n", img.width, img.height, vol->sz[2]);
    1.42 +
    1.43 +	/* allocate the whole volume at once */
    1.44 +	if(!(voxels = malloc(vol->sz[0] * vol->sz[1] * vol->sz[2] * 4 * sizeof *voxels))) {
    1.45 +		img_destroy(&img);
    1.46 +		goto err;
    1.47 +	}
    1.48 +	vptr = voxels;
    1.49 +	img_destroy(&img);
    1.50 +
    1.51 +	/* put the volume data into the alpha component */
    1.52 +	for(i=0; i<vol->sz[2]; i++) {
    1.53 +		int x, y, xsz, ysz;
    1.54 +		float *pixels, *src;
    1.55 +		struct slice *slice = slist;
    1.56 +		slist = slist->next;
    1.57 +
    1.58 +		if(!(pixels = img_load_pixels(slice->name, &xsz, &ysz, IMG_FMT_GREYF))) {
    1.59 +			fprintf(stderr, "failed to load volume slice: %s\n", slice->name);
    1.60 +			goto err;
    1.61 +		}
    1.62 +
    1.63 +		if(xsz != vol->sz[0] || ysz != vol->sz[1]) {
    1.64 +			fprintf(stderr, "inconsistent dimensions for slice %s: %dx%d (%dx%d expected)\n",
    1.65 +					slice->name, xsz, ysz, vol->sz[0], vol->sz[1]);
    1.66 +			goto err;
    1.67 +		}
    1.68 +		free(slice->name);
    1.69 +		free(slice);
    1.70 +
    1.71 +		src = pixels;
    1.72 +		for(y=0; y<ysz; y++) {
    1.73 +			for(x=0; x<xsz; x++) {
    1.74 +				vptr[0] = vptr[1] = vptr[2] = 0.0f;
    1.75 +				vptr[3] = *src++;
    1.76 +				vptr += 4;
    1.77 +			}
    1.78 +		}
    1.79 +		img_free_pixels(pixels);
    1.80 +	}
    1.81 +
    1.82 +	/* calculate gradients */
    1.83 +	/*calc_gradients(voxels, vol->sz[0], vol->sz[1], vol->sz[2]);*/
    1.84 +
    1.85 +	/* create the volume texture */
    1.86  	glGenTextures(1, &vol->tex_vol);
    1.87  	glBindTexture(GL_TEXTURE_3D, vol->tex_vol);
    1.88  	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    1.89 @@ -42,26 +101,12 @@
    1.90  	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    1.91  	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
    1.92  	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
    1.93 +	glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA32F, vol->sz[0], vol->sz[1], vol->sz[2], 0,
    1.94 +			GL_RGBA, GL_FLOAT, voxels);
    1.95  
    1.96 -	for(i=0; i<vol->sz[2]; i++) {
    1.97 -		struct img_pixmap img;
    1.98 -		struct slice *slice = head;
    1.99 -		head = head->next;
   1.100 -
   1.101 -		img_init(&img);
   1.102 -		if(img_load(&img, slice->name) == -1) {
   1.103 -			fprintf(stderr, "failed to load volume slice: %s\n", slice->name);
   1.104 -			goto err;
   1.105 -		}
   1.106 -	}
   1.107 -
   1.108 -	glGenTextures(1, &vol->tex_grad);
   1.109 -	glBindTexture(GL_TEXTURE_3D, vol->tex_grad);
   1.110 -	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
   1.111 -	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
   1.112 -	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
   1.113 -	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
   1.114 -	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
   1.115 +	/* ... and destroy our copy */
   1.116 +	free(voxels);
   1.117 +	return vol;
   1.118  
   1.119  err:
   1.120  	while(slist) {
   1.121 @@ -70,20 +115,54 @@
   1.122  		free(tmp->name);
   1.123  		free(tmp);
   1.124  	}
   1.125 +	free_volume(vol);
   1.126 +	free(voxels);
   1.127 +	return 0;
   1.128 +}
   1.129 +
   1.130 +void free_volume(struct volume *vol)
   1.131 +{
   1.132  	if(vol) {
   1.133 -		if(vol->tex_vol)
   1.134 +		if(vol->tex_vol) {
   1.135  			glDeleteTextures(1, &vol->tex_vol);
   1.136 -		if(vol->tex_grad)
   1.137 -			glDeleteTextures(1, &vol->tex_grad);
   1.138 +		}
   1.139  		free(vol);
   1.140  	}
   1.141 -	return 0;
   1.142 +}
   1.143 +
   1.144 +struct vec4 { float x, y, z, w; };
   1.145 +
   1.146 +static void calc_gradients(float *voxels, int xsz, int ysz, int zsz)
   1.147 +{
   1.148 +	int x, y, z, slice_pixels = xsz * ysz;
   1.149 +	struct vec4 *vptr = (struct vec4*)voxels;
   1.150 +
   1.151 +	for(z=0; z<zsz; z++) {
   1.152 +		for(y=0; y<ysz; y++) {
   1.153 +			for(x=0; x<xsz; x++) {
   1.154 +				float x0, x1, y0, y1, z0, z1;
   1.155 +
   1.156 +				x0 = x > 0 ? (vptr - 1)->w : vptr->w;
   1.157 +				y0 = y > 0 ? (vptr - xsz)->w : vptr->w;
   1.158 +				z0 = z > 0 ? (vptr - slice_pixels)->w : vptr->w;
   1.159 +				x1 = x < xsz - 1 ? (vptr + 1)->w : vptr->w;
   1.160 +				y1 = y < ysz - 1 ? (vptr + xsz)->w : vptr->w;
   1.161 +				z1 = z < zsz - 1 ? (vptr + slice_pixels)->w : vptr->w;
   1.162 +
   1.163 +				vptr->x = x1 - x0;
   1.164 +				vptr->y = y1 - y0;
   1.165 +				vptr->z = z1 - z0;
   1.166 +
   1.167 +				vptr++;
   1.168 +			}
   1.169 +		}
   1.170 +	}
   1.171  }
   1.172  
   1.173  static struct slice *read_voldesc(const char *fname, struct volume *vol)
   1.174  {
   1.175  	FILE *fp = 0;
   1.176 -	char buf[512];
   1.177 +	char buf[512], *slash, *prefix = "";
   1.178  	int mode_slices = 0;
   1.179  	struct slice *head = 0, *tail;
   1.180  
   1.181 @@ -97,6 +176,13 @@
   1.182  		goto err;
   1.183  	}
   1.184  
   1.185 +	if((slash = strrchr(fname, '/'))) {
   1.186 +		int len = slash - fname + 1;
   1.187 +		prefix = alloca(len + 1);
   1.188 +		memcpy(prefix, fname, len);
   1.189 +		prefix[len] = 0;
   1.190 +	}
   1.191 +
   1.192  	while(fgets(buf, sizeof buf, fp)) {
   1.193  		char *line = trim(buf);
   1.194  
   1.195 @@ -106,12 +192,12 @@
   1.196  		if(mode_slices) {
   1.197  			/* we're in the part that contains filenames... append to the list */
   1.198  			struct slice *node = malloc(sizeof *node);
   1.199 -			if(!node || !(node->name = malloc(strlen(line) + 1))) {
   1.200 +			if(!node || !(node->name = malloc(strlen(prefix) + strlen(line) + 1))) {
   1.201  				perror("failed to allocate list node");
   1.202  				free(node);
   1.203  				goto err;
   1.204  			}
   1.205 -			strcpy(node->name, line);
   1.206 +			sprintf(node->name, "%s%s", prefix, line);
   1.207  			node->next = 0;
   1.208  
   1.209  			if(head) {
   1.210 @@ -123,6 +209,10 @@
   1.211  			vol->sz[2]++;
   1.212  		} else {
   1.213  			/* TODO we're in the options part... parse */
   1.214 +
   1.215 +			if(strcmp(line, "SLICES") == 0) {
   1.216 +				mode_slices = 1;
   1.217 +			}
   1.218  		}
   1.219  	}
   1.220  	fclose(fp);