qvolray

annotate src/volume.c @ 9:a6765984e057

find changesets by author, revision, files, or words in the commit message
moved the volume loading to volume.c
author John Tsiombikas <nuclear@member.fsf.org>
date Sun, 08 Apr 2012 07:11:54 +0300
parents ae10631bb11b
children 9d2396738b60
rev   line source
nuclear@8 1 #include <stdio.h>
nuclear@8 2 #include <stdlib.h>
nuclear@8 3 #include <string.h>
nuclear@8 4 #include <ctype.h>
nuclear@8 5
nuclear@8 6 #ifndef __APPLE__
nuclear@8 7 #include <GL/gl.h>
nuclear@8 8 #else
nuclear@8 9 #include <OpenGL/gl.h>
nuclear@8 10 #endif
nuclear@8 11
nuclear@8 12 #include <imago2.h>
nuclear@9 13 #include "volume.h"
nuclear@8 14
nuclear@8 15 struct slice {
nuclear@8 16 char *name;
nuclear@8 17 struct slice *next;
nuclear@8 18 };
nuclear@8 19
nuclear@9 20 static void calc_gradients(float *voxels, int xsz, int ysz, int zsz);
nuclear@8 21 static struct slice *read_voldesc(const char *fname, struct volume *vol);
nuclear@8 22 static char *trim(char *str);
nuclear@8 23
nuclear@8 24 struct volume *load_volume(const char *fname)
nuclear@8 25 {
nuclear@9 26 int i;
nuclear@8 27 struct slice *slist;
nuclear@8 28 struct volume *vol = 0;
nuclear@9 29 float *voxels = 0, *vptr;
nuclear@9 30 struct img_pixmap img;
nuclear@8 31
nuclear@8 32 if(!(vol = malloc(sizeof *vol))) {
nuclear@8 33 perror("failed to allocate volume");
nuclear@8 34 return 0;
nuclear@8 35 }
nuclear@8 36 memset(vol, 0, sizeof *vol);
nuclear@8 37 vol->zaspect = 1;
nuclear@8 38
nuclear@8 39 if(!(slist = read_voldesc(fname, vol))) {
nuclear@8 40 goto err;
nuclear@8 41 }
nuclear@8 42
nuclear@9 43 /* load the first image to determine slice dimensions */
nuclear@9 44 img_init(&img);
nuclear@9 45 if(img_load(&img, slist->name) == -1) {
nuclear@9 46 fprintf(stderr, "failed to load volume slice: %s\n", slist->name);
nuclear@9 47 goto err;
nuclear@9 48 }
nuclear@9 49
nuclear@9 50 vol->sz[0] = img.width;
nuclear@9 51 vol->sz[1] = img.height;
nuclear@9 52 printf("volume dimensions: %dx%dx%d\n", img.width, img.height, vol->sz[2]);
nuclear@9 53
nuclear@9 54 /* allocate the whole volume at once */
nuclear@9 55 if(!(voxels = malloc(vol->sz[0] * vol->sz[1] * vol->sz[2] * 4 * sizeof *voxels))) {
nuclear@9 56 img_destroy(&img);
nuclear@9 57 goto err;
nuclear@9 58 }
nuclear@9 59 vptr = voxels;
nuclear@9 60 img_destroy(&img);
nuclear@9 61
nuclear@9 62 /* put the volume data into the alpha component */
nuclear@9 63 for(i=0; i<vol->sz[2]; i++) {
nuclear@9 64 int x, y, xsz, ysz;
nuclear@9 65 float *pixels, *src;
nuclear@9 66 struct slice *slice = slist;
nuclear@9 67 slist = slist->next;
nuclear@9 68
nuclear@9 69 if(!(pixels = img_load_pixels(slice->name, &xsz, &ysz, IMG_FMT_GREYF))) {
nuclear@9 70 fprintf(stderr, "failed to load volume slice: %s\n", slice->name);
nuclear@9 71 goto err;
nuclear@9 72 }
nuclear@9 73
nuclear@9 74 if(xsz != vol->sz[0] || ysz != vol->sz[1]) {
nuclear@9 75 fprintf(stderr, "inconsistent dimensions for slice %s: %dx%d (%dx%d expected)\n",
nuclear@9 76 slice->name, xsz, ysz, vol->sz[0], vol->sz[1]);
nuclear@9 77 goto err;
nuclear@9 78 }
nuclear@9 79 free(slice->name);
nuclear@9 80 free(slice);
nuclear@9 81
nuclear@9 82 src = pixels;
nuclear@9 83 for(y=0; y<ysz; y++) {
nuclear@9 84 for(x=0; x<xsz; x++) {
nuclear@9 85 vptr[0] = vptr[1] = vptr[2] = 0.0f;
nuclear@9 86 vptr[3] = *src++;
nuclear@9 87 vptr += 4;
nuclear@9 88 }
nuclear@9 89 }
nuclear@9 90 img_free_pixels(pixels);
nuclear@9 91 }
nuclear@9 92
nuclear@9 93 /* calculate gradients */
nuclear@9 94 /*calc_gradients(voxels, vol->sz[0], vol->sz[1], vol->sz[2]);*/
nuclear@9 95
nuclear@9 96 /* create the volume texture */
nuclear@8 97 glGenTextures(1, &vol->tex_vol);
nuclear@8 98 glBindTexture(GL_TEXTURE_3D, vol->tex_vol);
nuclear@8 99 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
nuclear@8 100 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
nuclear@8 101 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
nuclear@8 102 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
nuclear@8 103 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
nuclear@9 104 glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA32F, vol->sz[0], vol->sz[1], vol->sz[2], 0,
nuclear@9 105 GL_RGBA, GL_FLOAT, voxels);
nuclear@8 106
nuclear@9 107 /* ... and destroy our copy */
nuclear@9 108 free(voxels);
nuclear@9 109 return vol;
nuclear@8 110
nuclear@8 111 err:
nuclear@8 112 while(slist) {
nuclear@8 113 struct slice *tmp = slist;
nuclear@8 114 slist = slist->next;
nuclear@8 115 free(tmp->name);
nuclear@8 116 free(tmp);
nuclear@8 117 }
nuclear@9 118 free_volume(vol);
nuclear@9 119 free(voxels);
nuclear@9 120 return 0;
nuclear@9 121 }
nuclear@9 122
nuclear@9 123 void free_volume(struct volume *vol)
nuclear@9 124 {
nuclear@8 125 if(vol) {
nuclear@9 126 if(vol->tex_vol) {
nuclear@8 127 glDeleteTextures(1, &vol->tex_vol);
nuclear@9 128 }
nuclear@8 129 free(vol);
nuclear@8 130 }
nuclear@9 131 }
nuclear@9 132
nuclear@9 133 struct vec4 { float x, y, z, w; };
nuclear@9 134
nuclear@9 135 static void calc_gradients(float *voxels, int xsz, int ysz, int zsz)
nuclear@9 136 {
nuclear@9 137 int x, y, z, slice_pixels = xsz * ysz;
nuclear@9 138 struct vec4 *vptr = (struct vec4*)voxels;
nuclear@9 139
nuclear@9 140 for(z=0; z<zsz; z++) {
nuclear@9 141 for(y=0; y<ysz; y++) {
nuclear@9 142 for(x=0; x<xsz; x++) {
nuclear@9 143 float x0, x1, y0, y1, z0, z1;
nuclear@9 144
nuclear@9 145 x0 = x > 0 ? (vptr - 1)->w : vptr->w;
nuclear@9 146 y0 = y > 0 ? (vptr - xsz)->w : vptr->w;
nuclear@9 147 z0 = z > 0 ? (vptr - slice_pixels)->w : vptr->w;
nuclear@9 148 x1 = x < xsz - 1 ? (vptr + 1)->w : vptr->w;
nuclear@9 149 y1 = y < ysz - 1 ? (vptr + xsz)->w : vptr->w;
nuclear@9 150 z1 = z < zsz - 1 ? (vptr + slice_pixels)->w : vptr->w;
nuclear@9 151
nuclear@9 152 vptr->x = x1 - x0;
nuclear@9 153 vptr->y = y1 - y0;
nuclear@9 154 vptr->z = z1 - z0;
nuclear@9 155
nuclear@9 156 vptr++;
nuclear@9 157 }
nuclear@9 158 }
nuclear@9 159 }
nuclear@8 160 }
nuclear@8 161
nuclear@8 162 static struct slice *read_voldesc(const char *fname, struct volume *vol)
nuclear@8 163 {
nuclear@8 164 FILE *fp = 0;
nuclear@9 165 char buf[512], *slash, *prefix = "";
nuclear@8 166 int mode_slices = 0;
nuclear@8 167 struct slice *head = 0, *tail;
nuclear@8 168
nuclear@8 169 if(!(fp = fopen(fname, "r"))) {
nuclear@8 170 fprintf(stderr, "failed to open volume descriptor: %s\n", fname);
nuclear@8 171 return 0;
nuclear@8 172 }
nuclear@8 173 fgets(buf, sizeof buf, fp);
nuclear@8 174 if(strstr(buf, "VOLDESC") != buf) {
nuclear@8 175 fprintf(stderr, "invalid file format while trying to read volume descriptor: %s\n", fname);
nuclear@8 176 goto err;
nuclear@8 177 }
nuclear@8 178
nuclear@9 179 if((slash = strrchr(fname, '/'))) {
nuclear@9 180 int len = slash - fname + 1;
nuclear@9 181 prefix = alloca(len + 1);
nuclear@9 182 memcpy(prefix, fname, len);
nuclear@9 183 prefix[len] = 0;
nuclear@9 184 }
nuclear@9 185
nuclear@8 186 while(fgets(buf, sizeof buf, fp)) {
nuclear@8 187 char *line = trim(buf);
nuclear@8 188
nuclear@8 189 if(!*line || *line == '#')
nuclear@8 190 continue;
nuclear@8 191
nuclear@8 192 if(mode_slices) {
nuclear@8 193 /* we're in the part that contains filenames... append to the list */
nuclear@8 194 struct slice *node = malloc(sizeof *node);
nuclear@9 195 if(!node || !(node->name = malloc(strlen(prefix) + strlen(line) + 1))) {
nuclear@8 196 perror("failed to allocate list node");
nuclear@8 197 free(node);
nuclear@8 198 goto err;
nuclear@8 199 }
nuclear@9 200 sprintf(node->name, "%s%s", prefix, line);
nuclear@8 201 node->next = 0;
nuclear@8 202
nuclear@8 203 if(head) {
nuclear@8 204 tail->next = node;
nuclear@8 205 tail = node;
nuclear@8 206 } else {
nuclear@8 207 head = tail = node;
nuclear@8 208 }
nuclear@8 209 vol->sz[2]++;
nuclear@8 210 } else {
nuclear@8 211 /* TODO we're in the options part... parse */
nuclear@9 212
nuclear@9 213 if(strcmp(line, "SLICES") == 0) {
nuclear@9 214 mode_slices = 1;
nuclear@9 215 }
nuclear@8 216 }
nuclear@8 217 }
nuclear@8 218 fclose(fp);
nuclear@8 219
nuclear@8 220 if(!vol->sz[2]) {
nuclear@8 221 fprintf(stderr, "file: %s contains no slices\n", fname);
nuclear@8 222 goto err;
nuclear@8 223 }
nuclear@8 224
nuclear@8 225 return head;
nuclear@8 226
nuclear@8 227 err:
nuclear@8 228 while(head) {
nuclear@8 229 struct slice *tmp = head;
nuclear@8 230 head = head->next;
nuclear@8 231 free(tmp->name);
nuclear@8 232 free(tmp);
nuclear@8 233 }
nuclear@8 234 if(fp) {
nuclear@8 235 fclose(fp);
nuclear@8 236 }
nuclear@8 237 return 0;
nuclear@8 238 }
nuclear@8 239
nuclear@8 240 static char *trim(char *str)
nuclear@8 241 {
nuclear@8 242 char *tmp = str + strlen(str) - 1;
nuclear@8 243
nuclear@8 244 while(isspace(*tmp))
nuclear@8 245 tmp--;
nuclear@8 246 tmp[1] = 0;
nuclear@8 247
nuclear@8 248 tmp = str;
nuclear@8 249 while(isspace(*tmp))
nuclear@8 250 tmp++;
nuclear@8 251 return tmp;
nuclear@8 252 }