1.1 --- a/src/dicomview.cc	Sat Dec 27 02:35:58 2014 +0200
     1.2 +++ b/src/dicomview.cc	Sat Dec 27 06:32:28 2014 +0200
     1.3 @@ -1,21 +1,36 @@
     1.4  #include "opengl.h"
     1.5  #include "dicomview.h"
     1.6 +#include "rend_fast.h"
     1.7  
     1.8 +static int win_width, win_height;
     1.9  static float cam_theta, cam_phi, cam_dist = 6;
    1.10  
    1.11 +static Renderer *rend;
    1.12 +
    1.13  extern "C" {
    1.14  
    1.15  int init()
    1.16  {
    1.17 +	rend = new RendererFast;
    1.18 +	if(!rend->init()) {
    1.19 +		return -1;
    1.20 +	}
    1.21 +
    1.22 +	return 0;
    1.23  }
    1.24  
    1.25  void cleanup()
    1.26  {
    1.27 +	rend->destroy();
    1.28 +	delete rend;
    1.29  }
    1.30  
    1.31  void ev_display()
    1.32  {
    1.33  	glClear(GL_COLOR_BUFFER_BIT);
    1.34 +
    1.35 +	rend->update(0);
    1.36 +	rend->render();
    1.37  }
    1.38  
    1.39  void ev_reshape(int x, int y)
    1.40 @@ -23,11 +38,25 @@
    1.41  	win_width = x;
    1.42  	win_height = y;
    1.43  
    1.44 -	glViweport(0, 0, x, y);
    1.45 +	glViewport(0, 0, x, y);
    1.46  }
    1.47  
    1.48 -void ev_keyboard(int key, int press, int x, int y);
    1.49 -void ev_mouse_button(int bn, int press, int x, int y);
    1.50 -void ev_mouse_motion(int x, int y);
    1.51 +void ev_keyboard(int key, int press, int x, int y)
    1.52 +{
    1.53 +	if(press) {
    1.54 +		switch(key) {
    1.55 +		case 27:
    1.56 +			quit();
    1.57 +		}
    1.58 +	}
    1.59 +}
    1.60 +
    1.61 +void ev_mouse_button(int bn, int press, int x, int y)
    1.62 +{
    1.63 +}
    1.64 +
    1.65 +void ev_mouse_motion(int x, int y)
    1.66 +{
    1.67 +}
    1.68  
    1.69  }	// extern "C"

     2.1 --- a/src/image.cc	Sat Dec 27 02:35:58 2014 +0200
     2.2 +++ b/src/image.cc	Sat Dec 27 06:32:28 2014 +0200
     2.3 @@ -1,3 +1,4 @@
     2.4 +#include <string.h>
     2.5  #include <imago2.h>
     2.6  #include "image.h"
     2.7  
     2.8 @@ -39,7 +40,7 @@
     2.9  	// TODO dicom loader ... ?
    2.10  
    2.11  	int x, y;
    2.12 -	float *pixels = img_load_pixels(fname, &x, &y, IMG_FMT_GREYF);
    2.13 +	float *pixels = (float*)img_load_pixels(fname, &x, &y, IMG_FMT_GREYF);
    2.14  	if(!pixels) {
    2.15  		return false;
    2.16  	}

     3.1 --- a/src/main.cc	Sat Dec 27 02:35:58 2014 +0200
     3.2 +++ b/src/main.cc	Sat Dec 27 06:32:28 2014 +0200
     3.3 @@ -19,8 +19,6 @@
     3.4  static void mouse(int bn, int state, int x, int y);
     3.5  static void motion(int x, int y);
     3.6  
     3.7 -static float cam_theta, cam_phi, cam_dist = 6;
     3.8 -
     3.9  int main(int argc, char **argv)
    3.10  {
    3.11  	glutInit(&argc, argv);
    3.12 @@ -30,7 +28,8 @@
    3.13  
    3.14  	glutDisplayFunc(display);
    3.15  	glutReshapeFunc(reshape);
    3.16 -	glutKeyboardFunc(keyb);
    3.17 +	glutKeyboardFunc(key_down);
    3.18 +	glutKeyboardUpFunc(key_up);
    3.19  	glutMouseFunc(mouse);
    3.20  	glutMotionFunc(motion);
    3.21  

     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/src/volume.cc	Sat Dec 27 06:32:28 2014 +0200
     4.3 @@ -0,0 +1,198 @@
     4.4 +#include <stdio.h>
     4.5 +#include <string.h>
     4.6 +#include <ctype.h>
     4.7 +#include <errno.h>
     4.8 +#include <math.h>
     4.9 +#include "volume.h"
    4.10 +
    4.11 +static char *strip_space(char *s);
    4.12 +
    4.13 +
    4.14 +Volume::~Volume()
    4.15 +{
    4.16 +}
    4.17 +
    4.18 +int Volume::num_samples(int dim) const
    4.19 +{
    4.20 +	return 0;
    4.21 +}
    4.22 +
    4.23 +void Volume::normalf(float *norm, float x, float y, float z, float delta)
    4.24 +{
    4.25 +	float dx, dy, dz;
    4.26 +	dx = dy = dz = delta;
    4.27 +
    4.28 +	if(num_samples(0) > 0) {
    4.29 +		// discrete volume
    4.30 +		dx /= (float)num_samples(0);
    4.31 +		dy /= (float)num_samples(1);
    4.32 +		dz /= (float)num_samples(2);
    4.33 +	}
    4.34 +
    4.35 +	norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z);
    4.36 +	norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z);
    4.37 +	norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz);
    4.38 +
    4.39 +	float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
    4.40 +	if(len != 0.0) {
    4.41 +		norm[0] /= len;
    4.42 +		norm[1] /= len;
    4.43 +		norm[2] /= len;
    4.44 +	}
    4.45 +}
    4.46 +
    4.47 +void Volume::normali(float *norm, int x, int y, int z)
    4.48 +{
    4.49 +	int sz[3];
    4.50 +
    4.51 +	if((sz[0] = num_samples(0)) <= 0) {
    4.52 +		// if it's a continuous volume, just call normalf
    4.53 +		normalf(norm, x, y, z);
    4.54 +		return;
    4.55 +	}
    4.56 +	sz[1] = num_samples(1);
    4.57 +	sz[2] = num_samples(2);
    4.58 +
    4.59 +	int prevx = x <= 0 ? 0 : x;
    4.60 +	int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x;
    4.61 +	int prevy = y <= 0 ? 0 : y;
    4.62 +	int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y;
    4.63 +	int prevz = z <= 0 ? 0 : z;
    4.64 +	int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z;
    4.65 +
    4.66 +	norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z);
    4.67 +	norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z);
    4.68 +	norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz);
    4.69 +
    4.70 +	float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
    4.71 +	if(len != 0.0) {
    4.72 +		norm[0] /= len;
    4.73 +		norm[1] /= len;
    4.74 +		norm[2] /= len;
    4.75 +	}
    4.76 +}
    4.77 +
    4.78 +// ---- VoxelVolume (discrete) implementation ----
    4.79 +VoxelVolume::VoxelVolume()
    4.80 +{
    4.81 +	size[0] = size[1] = size[2] = 0;
    4.82 +}
    4.83 +
    4.84 +bool VoxelVolume::load(const char *fname)
    4.85 +{
    4.86 +	FILE *fp = fopen(fname, "r");
    4.87 +	if(!fp) {
    4.88 +		fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno));
    4.89 +		return false;
    4.90 +	}
    4.91 +
    4.92 +	char buf[512];
    4.93 +	while(fgets(buf, sizeof buf, fp)) {
    4.94 +		char *line = strip_space(buf);
    4.95 +
    4.96 +		Image img;
    4.97 +		if(!img.load(line)) {
    4.98 +			slices.clear();
    4.99 +			return false;
   4.100 +		}
   4.101 +
   4.102 +		if(slices.empty()) {
   4.103 +			size[0] = img.width;
   4.104 +			size[1] = img.height;
   4.105 +		} else {
   4.106 +			if(img.width != size[0] || img.height != size[1]) {
   4.107 +				fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n",
   4.108 +						(int)slices.size(), line, img.width, img.height, size[0], size[1]);
   4.109 +				slices.clear();
   4.110 +				return false;
   4.111 +			}
   4.112 +		}
   4.113 +
   4.114 +		slices.push_back(img);
   4.115 +		img.pixels = 0;	// otherwise the destructor will free it
   4.116 +	}
   4.117 +
   4.118 +	size[2] = slices.size();
   4.119 +	fclose(fp);
   4.120 +	return true;
   4.121 +}
   4.122 +
   4.123 +int VoxelVolume::num_samples(int dim) const
   4.124 +{
   4.125 +	return size[dim];
   4.126 +}
   4.127 +
   4.128 +static inline int clamp(int x, int low, int high)
   4.129 +{
   4.130 +	return x < low ? low : (x > high ? high : x);
   4.131 +}
   4.132 +
   4.133 +static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty,
   4.134 +		float *pixels, int width, int height)
   4.135 +{
   4.136 +	float v00, v01, v10, v11;
   4.137 +
   4.138 +	v00 = pixels[y0 * width + x0];
   4.139 +	v01 = pixels[y1 * width + x0];
   4.140 +	v10 = pixels[y0 * width + x1];
   4.141 +	v11 = pixels[y1 * width + x1];
   4.142 +
   4.143 +	float v0 = v00 + (v01 - v00) * ty;
   4.144 +	float v1 = v10 + (v11 - v10) * ty;
   4.145 +
   4.146 +	return v0 + (v1 - v0) * tx;
   4.147 +}
   4.148 +
   4.149 +float VoxelVolume::valuef(float x, float y, float z) const
   4.150 +{
   4.151 +	if(slices.empty()) {
   4.152 +		return 0.0f;
   4.153 +	}
   4.154 +
   4.155 +	float floor_x = floor(x);
   4.156 +	float ceil_x = ceil(x);
   4.157 +	float tx = x - floor_x;
   4.158 +
   4.159 +	float floor_y = floor(y);
   4.160 +	float ceil_y = ceil(y);
   4.161 +	float ty = y - floor_y;
   4.162 +
   4.163 +	float floor_z = floor(z);
   4.164 +	float ceil_z = ceil(z);
   4.165 +	float tz = z - floor_z;
   4.166 +
   4.167 +	int x0 = clamp(floor_x, 0, size[0] - 1);
   4.168 +	int x1 = clamp(ceil_x, 0, size[0] - 1);
   4.169 +
   4.170 +	int y0 = clamp(floor_y, 0, size[1] - 1);
   4.171 +	int y1 = clamp(ceil_y, 0, size[1] - 1);
   4.172 +
   4.173 +	int s0 = clamp(floor_z, 0, size[2] - 1);
   4.174 +	int s1 = clamp(ceil_z, 0, size[2] - 1);
   4.175 +
   4.176 +	float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]);
   4.177 +	float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]);
   4.178 +
   4.179 +	return val_s0 + (val_s1 - val_s0) * tz;
   4.180 +}
   4.181 +
   4.182 +float VoxelVolume::valuei(int x, int y, int z) const
   4.183 +{
   4.184 +	x = clamp(x, 0, size[0] - 1);
   4.185 +	y = clamp(y, 0, size[1] - 1);
   4.186 +	z = clamp(z, 0, size[2] - 1);
   4.187 +
   4.188 +	return slices[z].pixels[y * size[0] + x];
   4.189 +}
   4.190 +
   4.191 +static char *strip_space(char *s)
   4.192 +{
   4.193 +	while(*s && isspace(*s)) s++;
   4.194 +	if(!*s) return 0;
   4.195 +
   4.196 +	char *end = s + strlen(s) - 1;
   4.197 +	while(end > s && isspace(*end)) end--;
   4.198 +	end[1] = 0;
   4.199 +
   4.200 +	return s;
   4.201 +}

     5.1 --- a/src/volume.h	Sat Dec 27 02:35:58 2014 +0200
     5.2 +++ b/src/volume.h	Sat Dec 27 06:32:28 2014 +0200
     5.3 @@ -1,16 +1,18 @@
     5.4  #ifndef VOLUME_H_
     5.5  #define VOLUME_H_
     5.6  
     5.7 +#include <vector>
     5.8  #include "image.h"
     5.9  
    5.10  class Volume {
    5.11  public:
    5.12  	virtual ~Volume();
    5.13  
    5.14 -	// returns -1 if the volume is continuously defined
    5.15 +	// returns 0 if the volume is continuously defined
    5.16  	virtual int num_samples(int dim) const;
    5.17  
    5.18 -	// central differences offset will be delta / num_samples
    5.19 +	// central differences offset will be delta / num_samples for discrete volumes
    5.20 +	// and delta for continuous volumes
    5.21  	virtual void normalf(float *norm, float x, float y, float z, float delta = 1.0);
    5.22  	virtual void normali(float *norm, int x, int y, int z);
    5.23  
    5.24 @@ -20,11 +22,11 @@
    5.25  
    5.26  class VoxelVolume : public Volume {
    5.27  protected:
    5.28 +	int size[3];
    5.29  	std::vector<Image> slices;
    5.30  
    5.31  public:
    5.32  	VoxelVolume();
    5.33 -	~VoxelVolume();
    5.34  
    5.35  	bool load(const char *fname);
    5.36