qvolray

annotate src/volume.cc @ 21:4c62be57fc1a

find changesets by author, revision, files, or words in the commit message
foo
author John Tsiombikas <nuclear@member.fsf.org>
date Wed, 11 Apr 2012 16:59:45 +0300
parents 535762131d34
children 450d4c50470f
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@11 5 #include <list>
nuclear@11 6 #include <string>
nuclear@17 7 #include <GL/glew.h>
nuclear@8 8 #include <imago2.h>
nuclear@9 9 #include "volume.h"
nuclear@8 10
nuclear@9 11 static void calc_gradients(float *voxels, int xsz, int ysz, int zsz);
nuclear@8 12 static struct slice *read_voldesc(const char *fname, struct volume *vol);
nuclear@8 13 static char *trim(char *str);
nuclear@8 14
nuclear@11 15 Volume::Volume()
nuclear@8 16 {
nuclear@11 17 sz[0] = sz[1] = sz[2] = 0;
nuclear@11 18 zaspect = 1.0;
nuclear@11 19 tex = 0;
nuclear@11 20 }
nuclear@11 21
nuclear@11 22 Volume::~Volume()
nuclear@11 23 {
nuclear@11 24 if(tex) {
nuclear@11 25 glDeleteTextures(1, &tex);
nuclear@11 26 }
nuclear@11 27 }
nuclear@11 28
nuclear@18 29 bool Volume::create(int xsz, int ysz, int zsz, float *data)
nuclear@18 30 {
nuclear@18 31 if(!tex)
nuclear@18 32 glGenTextures(1, &tex);
nuclear@18 33 glBindTexture(GL_TEXTURE_3D, tex);
nuclear@18 34 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
nuclear@18 35 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
nuclear@18 36 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
nuclear@18 37 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
nuclear@18 38 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
nuclear@18 39 glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA32F_ARB, xsz, ysz, zsz, 0, GL_RGBA, GL_FLOAT, data);
nuclear@18 40 }
nuclear@18 41
nuclear@11 42 bool Volume::load(const char *fname)
nuclear@11 43 {
nuclear@11 44 std::list<std::string> slist;
nuclear@9 45 float *voxels = 0, *vptr;
nuclear@9 46 struct img_pixmap img;
nuclear@8 47
nuclear@11 48 if(!(read_voldesc(fname, &slist))) {
nuclear@11 49 return false;
nuclear@8 50 }
nuclear@8 51
nuclear@9 52 /* load the first image to determine slice dimensions */
nuclear@9 53 img_init(&img);
nuclear@13 54 if(img_load(&img, slist.begin()->c_str()) == -1) {
nuclear@13 55 fprintf(stderr, "failed to load volume slice: %s\n", slist.begin()->c_str());
nuclear@11 56 return false;
nuclear@9 57 }
nuclear@9 58
nuclear@11 59 sz[0] = img.width;
nuclear@11 60 sz[1] = img.height;
nuclear@13 61 printf("volume dimensions: %dx%dx%d\n", img.width, img.height, sz[2]);
nuclear@9 62
nuclear@9 63 /* allocate the whole volume at once */
nuclear@11 64 voxels = new float[sz[0] * sz[1] * sz[2] * 4];
nuclear@9 65 vptr = voxels;
nuclear@9 66 img_destroy(&img);
nuclear@9 67
nuclear@9 68 /* put the volume data into the alpha component */
nuclear@17 69 //for(auto slice : slist) { // fucking vs2010 doesn't support range-for
nuclear@17 70 for(auto iter = slist.begin(); iter != slist.end(); iter++) {
nuclear@17 71 auto slice = *iter;
nuclear@17 72
nuclear@9 73 int x, y, xsz, ysz;
nuclear@9 74 float *pixels, *src;
nuclear@9 75
nuclear@13 76 if(!(pixels = (float*)img_load_pixels(slice.c_str(), &xsz, &ysz, IMG_FMT_GREYF))) {
nuclear@13 77 fprintf(stderr, "failed to load volume slice: %s\n", slice.c_str());
nuclear@11 78 delete [] voxels;
nuclear@11 79 return false;
nuclear@9 80 }
nuclear@9 81
nuclear@11 82 if(xsz != sz[0] || ysz != sz[1]) {
nuclear@9 83 fprintf(stderr, "inconsistent dimensions for slice %s: %dx%d (%dx%d expected)\n",
nuclear@13 84 slice.c_str(), xsz, ysz, sz[0], sz[1]);
nuclear@11 85 delete [] voxels;
nuclear@11 86 return false;
nuclear@9 87 }
nuclear@9 88
nuclear@9 89 src = pixels;
nuclear@9 90 for(y=0; y<ysz; y++) {
nuclear@9 91 for(x=0; x<xsz; x++) {
nuclear@9 92 vptr[0] = vptr[1] = vptr[2] = 0.0f;
nuclear@9 93 vptr[3] = *src++;
nuclear@9 94 vptr += 4;
nuclear@9 95 }
nuclear@9 96 }
nuclear@9 97 img_free_pixels(pixels);
nuclear@9 98 }
nuclear@9 99
nuclear@9 100 /* calculate gradients */
nuclear@11 101 calc_gradients(voxels, sz[0], sz[1], sz[2]);
nuclear@9 102
nuclear@9 103 /* create the volume texture */
nuclear@18 104 create(sz[0], sz[1], sz[2], voxels);
nuclear@8 105
nuclear@9 106 /* ... and destroy our copy */
nuclear@11 107 delete [] voxels;
nuclear@11 108 return true;
nuclear@9 109 }
nuclear@9 110
nuclear@11 111 unsigned int Volume::get_texture() const
nuclear@9 112 {
nuclear@11 113 return tex;
nuclear@11 114 }
nuclear@11 115
nuclear@11 116 void Volume::bind(int tex_unit) const
nuclear@11 117 {
nuclear@11 118 glActiveTexture(GL_TEXTURE0 + tex_unit);
nuclear@11 119 glBindTexture(GL_TEXTURE_3D, tex);
nuclear@9 120 }
nuclear@9 121
nuclear@9 122 struct vec4 { float x, y, z, w; };
nuclear@9 123
nuclear@9 124 static void calc_gradients(float *voxels, int xsz, int ysz, int zsz)
nuclear@9 125 {
nuclear@9 126 int x, y, z, slice_pixels = xsz * ysz;
nuclear@11 127 vec4 *vptr = (vec4*)voxels;
nuclear@9 128
nuclear@9 129 for(z=0; z<zsz; z++) {
nuclear@9 130 for(y=0; y<ysz; y++) {
nuclear@9 131 for(x=0; x<xsz; x++) {
nuclear@9 132 float x0, x1, y0, y1, z0, z1;
nuclear@9 133
nuclear@9 134 x0 = x > 0 ? (vptr - 1)->w : vptr->w;
nuclear@9 135 y0 = y > 0 ? (vptr - xsz)->w : vptr->w;
nuclear@9 136 z0 = z > 0 ? (vptr - slice_pixels)->w : vptr->w;
nuclear@9 137 x1 = x < xsz - 1 ? (vptr + 1)->w : vptr->w;
nuclear@9 138 y1 = y < ysz - 1 ? (vptr + xsz)->w : vptr->w;
nuclear@9 139 z1 = z < zsz - 1 ? (vptr + slice_pixels)->w : vptr->w;
nuclear@9 140
nuclear@9 141 vptr->x = x1 - x0;
nuclear@9 142 vptr->y = y1 - y0;
nuclear@9 143 vptr->z = z1 - z0;
nuclear@9 144
nuclear@9 145 vptr++;
nuclear@9 146 }
nuclear@9 147 }
nuclear@9 148 }
nuclear@8 149 }
nuclear@8 150
nuclear@13 151 bool Volume::read_voldesc(const char *fname, std::list<std::string> *slist)
nuclear@8 152 {
nuclear@8 153 FILE *fp = 0;
nuclear@9 154 char buf[512], *slash, *prefix = "";
nuclear@8 155 int mode_slices = 0;
nuclear@11 156
nuclear@11 157 slist->clear();
nuclear@8 158
nuclear@8 159 if(!(fp = fopen(fname, "r"))) {
nuclear@8 160 fprintf(stderr, "failed to open volume descriptor: %s\n", fname);
nuclear@11 161 return false;
nuclear@8 162 }
nuclear@8 163 fgets(buf, sizeof buf, fp);
nuclear@8 164 if(strstr(buf, "VOLDESC") != buf) {
nuclear@8 165 fprintf(stderr, "invalid file format while trying to read volume descriptor: %s\n", fname);
nuclear@11 166 fclose(fp);
nuclear@11 167 return false;
nuclear@8 168 }
nuclear@8 169
nuclear@13 170 if((slash = strrchr((char*)fname, '/'))) {
nuclear@9 171 int len = slash - fname + 1;
nuclear@13 172 prefix = (char*)alloca(len + 1);
nuclear@9 173 memcpy(prefix, fname, len);
nuclear@9 174 prefix[len] = 0;
nuclear@9 175 }
nuclear@9 176
nuclear@8 177 while(fgets(buf, sizeof buf, fp)) {
nuclear@8 178 char *line = trim(buf);
nuclear@8 179
nuclear@8 180 if(!*line || *line == '#')
nuclear@8 181 continue;
nuclear@8 182
nuclear@8 183 if(mode_slices) {
nuclear@8 184 /* we're in the part that contains filenames... append to the list */
nuclear@11 185 slist->push_back(std::string(prefix) + std::string(line));
nuclear@11 186 sz[2]++;
nuclear@8 187 } else {
nuclear@8 188 /* TODO we're in the options part... parse */
nuclear@9 189 if(strcmp(line, "SLICES") == 0) {
nuclear@9 190 mode_slices = 1;
nuclear@9 191 }
nuclear@8 192 }
nuclear@8 193 }
nuclear@8 194 fclose(fp);
nuclear@8 195
nuclear@11 196 if(slist->empty()) {
nuclear@8 197 fprintf(stderr, "file: %s contains no slices\n", fname);
nuclear@11 198 return false;
nuclear@8 199 }
nuclear@11 200 return true;
nuclear@8 201 }
nuclear@8 202
nuclear@8 203 static char *trim(char *str)
nuclear@8 204 {
nuclear@8 205 char *tmp = str + strlen(str) - 1;
nuclear@8 206
nuclear@8 207 while(isspace(*tmp))
nuclear@8 208 tmp--;
nuclear@8 209 tmp[1] = 0;
nuclear@8 210
nuclear@8 211 tmp = str;
nuclear@8 212 while(isspace(*tmp))
nuclear@8 213 tmp++;
nuclear@8 214 return tmp;
nuclear@8 215 }