volmetrics
changeset 0:88d390af583f
Image loader
author | Eleni Maria Stea <eleni@mutantstargoat.com> |
---|---|
date | Sat, 11 Jan 2014 17:22:36 +0200 |
parents | |
children | cca2e05dbabe |
files | Makefile src/image.h src/main.cc src/matrix.cc src/matrix.h src/vector.cc src/vector.h src/volume.cc src/volume.h |
diffstat | 9 files changed, 958 insertions(+), 0 deletions(-) [+] |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/Makefile Sat Jan 11 17:22:36 2014 +0200 1.3 @@ -0,0 +1,24 @@ 1.4 +src = $(wildcard src/*.cc src/math/*.cc src/shaders/*.cc) 1.5 +obj = $(src:.cc=.o) 1.6 +dep = $(obj:.o=.d) 1.7 +bin = test1 1.8 + 1.9 +dbg = -g 1.10 +opt = -O0 1.11 +inc = -Isrc -I/usr/local/include -I/usr/local/lib 1.12 + 1.13 +CXX = g++ 1.14 +CXXFLAGS = -pedantic -Wall $(dbg) $(opt) $(inc) 1.15 +LDFLAGS = -lGL -lGLU -lGLEW -lglut -limago -lm 1.16 + 1.17 +$(bin): $(obj) 1.18 + $(CXX) -o $@ $(obj) $(LDFLAGS) 1.19 + 1.20 +-include $(dep) 1.21 + 1.22 +%.d: %.cc 1.23 + @$(CPP) $(CXXFLAGS) $< -MM -MT $(@:.d=.o) >$@ 1.24 + 1.25 +.PHONY: clean 1.26 +clean: 1.27 + rm -f $(obj) $(bin) $(dep)
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 2.2 +++ b/src/image.h Sat Jan 11 17:22:36 2014 +0200 2.3 @@ -0,0 +1,23 @@ 2.4 +#ifndef IMAGE_H_ 2.5 +#define IMAGE_H_ 2.6 + 2.7 +class Image { 2.8 +private: 2.9 + float *pixels; 2.10 + int width; 2.11 + int height; 2.12 +public: 2.13 + Image(); 2.14 + ~Image(); 2.15 + 2.16 + bool load(const char *fname); 2.17 + 2.18 + float *get_pixels(); 2.19 + const float *get_pixels() const; 2.20 + void set_pixels(const float *pixels, int width, int height); 2.21 + 2.22 + int get_width() const; 2.23 + int get_height() const; 2.24 +}; 2.25 + 2.26 +#endif //IMAGE_H_
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 3.2 +++ b/src/main.cc Sat Jan 11 17:22:36 2014 +0200 3.3 @@ -0,0 +1,87 @@ 3.4 +#include <GL/glew.h> 3.5 +#include <GL/glut.h> 3.6 + 3.7 +#include <stdio.h> 3.8 +#include <assert.h> 3.9 + 3.10 +//static void init(void); 3.11 +static void display(void); 3.12 +static void reshape(int x, int y); 3.13 +static void keyboard(unsigned char key, int x, int y); 3.14 +static void mouse(int button, int state, int x, int y); 3.15 +static void motion(int x, int y); 3.16 + 3.17 +static int win_xsz, win_ysz; 3.18 + 3.19 +int main(int argc, char **argv) 3.20 +{ 3.21 + glutInit(&argc, argv); 3.22 + glutInitWindowSize(1280, 720); 3.23 + glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE); 3.24 + 3.25 + //TODO parse arguments 3.26 + 3.27 + glutCreateWindow("My Colonoscopie OEO!"); 3.28 + 3.29 + glutDisplayFunc(display); 3.30 + glutReshapeFunc(reshape); 3.31 + glutKeyboardFunc(keyboard); 3.32 + glutMouseFunc(mouse); 3.33 + glutMotionFunc(motion); 3.34 + 3.35 + glewInit(); 3.36 + 3.37 + //call init 3.38 + 3.39 + glutMainLoop(); 3.40 + return 0; 3.41 +} 3.42 + 3.43 +int init(void) 3.44 +{ 3.45 + //TODO 3.46 + return 0; 3.47 +} 3.48 + 3.49 +void display(void) 3.50 +{ 3.51 + //render 3.52 + 3.53 + glutSwapBuffers(); 3.54 + assert(glGetError() == GL_NO_ERROR); 3.55 +} 3.56 + 3.57 +void reshape(int x, int y) 3.58 +{ 3.59 + glViewport(0, 0, x, y); 3.60 + if(x != win_xsz || y != win_ysz) { 3.61 + //TODO raytex_needs_recalc = 1; 3.62 + win_xsz = x; 3.63 + win_ysz = y; 3.64 + } 3.65 +} 3.66 + 3.67 +void keyboard(unsigned char key, int x, int y) 3.68 +{ 3.69 + switch(key) { 3.70 + case 27: 3.71 + exit(0); 3.72 + default: 3.73 + break; 3.74 + } 3.75 +} 3.76 + 3.77 +static int prev_x, prev_y; 3.78 +void mouse(int bn, int state, int x, int y) 3.79 +{ 3.80 + prev_x = x; 3.81 + prev_y = y; 3.82 +} 3.83 + 3.84 +void motion(int x, int y) 3.85 +{ 3.86 + int dx = x - prev_x; 3.87 + int dy = y - prev_y; 3.88 + prev_x = x; 3.89 + prev_y = y; 3.90 +}
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 4.2 +++ b/src/matrix.cc Sat Jan 11 17:22:36 2014 +0200 4.3 @@ -0,0 +1,393 @@ 4.4 +#include <math.h> 4.5 +#include <stdio.h> 4.6 +#include <string.h> 4.7 +#include "matrix.h" 4.8 + 4.9 +Matrix4x4 Matrix4x4::identity = Matrix4x4(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1); 4.10 + 4.11 +Matrix4x4::Matrix4x4() 4.12 +{ 4.13 + *this = identity; 4.14 +} 4.15 + 4.16 +Matrix4x4::Matrix4x4( float m11, float m12, float m13, float m14, 4.17 + float m21, float m22, float m23, float m24, 4.18 + float m31, float m32, float m33, float m34, 4.19 + float m41, float m42, float m43, float m44) 4.20 +{ 4.21 + m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[0][3] = m14; 4.22 + m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[1][3] = m24; 4.23 + m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; m[2][3] = m34; 4.24 + m[3][0] = m41; m[3][1] = m42; m[3][2] = m43; m[3][3] = m44; 4.25 +} 4.26 + 4.27 +Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2) 4.28 +{ 4.29 + Matrix4x4 res; 4.30 + const float *op1 = m1.m[0], *op2 = m2.m[0]; 4.31 + float *dest = res.m[0]; 4.32 + 4.33 + for(int i=0; i<16; i++) { 4.34 + *dest++ = *op1++ + *op2++; 4.35 + } 4.36 + return res; 4.37 +} 4.38 + 4.39 +Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2) 4.40 +{ 4.41 + Matrix4x4 res; 4.42 + const float *op1 = m1.m[0], *op2 = m2.m[0]; 4.43 + float *dest = res.m[0]; 4.44 + 4.45 + for(int i=0; i<16; i++) { 4.46 + *dest++ = *op1++ - *op2++; 4.47 + } 4.48 + return res; 4.49 +} 4.50 + 4.51 +Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2) { 4.52 + Matrix4x4 res; 4.53 + 4.54 + for(int i=0; i<4; i++) { 4.55 + for(int j=0; j<4; j++) { 4.56 + res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j]; 4.57 + } 4.58 + } 4.59 + 4.60 + return res; 4.61 +} 4.62 + 4.63 +void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2) 4.64 +{ 4.65 + float *op1 = m1.m[0]; 4.66 + const float *op2 = m2.m[0]; 4.67 + 4.68 + for(int i=0; i<16; i++) { 4.69 + *op1++ += *op2++; 4.70 + } 4.71 +} 4.72 + 4.73 +void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2) 4.74 +{ 4.75 + float *op1 = m1.m[0]; 4.76 + const float *op2 = m2.m[0]; 4.77 + 4.78 + for(int i=0; i<16; i++) { 4.79 + *op1++ -= *op2++; 4.80 + } 4.81 +} 4.82 + 4.83 +void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2) 4.84 +{ 4.85 + Matrix4x4 res; 4.86 + for(int i=0; i<4; i++) { 4.87 + for(int j=0; j<4; j++) { 4.88 + res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j]; 4.89 + } 4.90 + } 4.91 + memcpy(m1.m, res.m, 16 * sizeof(float)); 4.92 +} 4.93 + 4.94 +Matrix4x4 operator *(const Matrix4x4 &mat, float scalar) 4.95 +{ 4.96 + Matrix4x4 res; 4.97 + const float *mptr = mat.m[0]; 4.98 + float *dptr = res.m[0]; 4.99 + 4.100 + for(int i=0; i<16; i++) { 4.101 + *dptr++ = *mptr++ * scalar; 4.102 + } 4.103 + return res; 4.104 +} 4.105 + 4.106 +Matrix4x4 operator *(float scalar, const Matrix4x4 &mat) 4.107 +{ 4.108 + Matrix4x4 res; 4.109 + const float *mptr = mat.m[0]; 4.110 + float *dptr = res.m[0]; 4.111 + 4.112 + for(int i=0; i<16; i++) { 4.113 + *dptr++ = *mptr++ * scalar; 4.114 + } 4.115 + return res; 4.116 +} 4.117 + 4.118 +void operator *=(Matrix4x4 &mat, float scalar) 4.119 +{ 4.120 + float *mptr = mat.m[0]; 4.121 + 4.122 + for(int i=0; i<16; i++) { 4.123 + *mptr++ *= scalar; 4.124 + } 4.125 +} 4.126 + 4.127 +void Matrix4x4::translate(const Vector3 &trans) 4.128 +{ 4.129 + Matrix4x4 tmat(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1); 4.130 + *this *= tmat; 4.131 +} 4.132 + 4.133 +void Matrix4x4::set_translation(const Vector3 &trans) 4.134 +{ 4.135 + *this = Matrix4x4(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1); 4.136 +} 4.137 + 4.138 +void Matrix4x4::rotate(const Vector3 &euler) 4.139 +{ 4.140 + float sinx = sin(euler.x); 4.141 + float cosx = cos(euler.x); 4.142 + float siny = sin(euler.y); 4.143 + float cosy = cos(euler.y); 4.144 + float sinz = sin(euler.z); 4.145 + float cosz = cos(euler.z); 4.146 + 4.147 + Matrix4x4 xrot = Matrix4x4( 4.148 + 1, 0, 0, 0, 4.149 + 0, cosx, sinx, 0, 4.150 + 0, -sinx, cosx, 0, 4.151 + 0, 0, 0, 1); 4.152 + 4.153 + Matrix4x4 yrot = Matrix4x4( 4.154 + cosy, 0, -siny, 0, 4.155 + 0, 1, 0, 0, 4.156 + siny, 0, cosy, 0, 4.157 + 0, 0, 0, 1); 4.158 + 4.159 + Matrix4x4 zrot = Matrix4x4( 4.160 + cosz, sinz, 0, 0, 4.161 + -sinz, cosz, 0, 0, 4.162 + 0, 0, 1, 0, 4.163 + 0, 0, 0, 1); 4.164 + 4.165 + *this *= (xrot * yrot * zrot); 4.166 +} 4.167 + 4.168 +void Matrix4x4::rotate(const Vector3 &axis, float angle) 4.169 +{ 4.170 + float sina = (float)sin(angle); 4.171 + float cosa = (float)cos(angle); 4.172 + float invcosa = 1-cosa; 4.173 + float nxsq = axis.x * axis.x; 4.174 + float nysq = axis.y * axis.y; 4.175 + float nzsq = axis.z * axis.z; 4.176 + 4.177 + Matrix4x4 xform; 4.178 + xform[0][0] = nxsq + (1-nxsq) * cosa; 4.179 + xform[0][1] = axis.x * axis.y * invcosa - axis.z * sina; 4.180 + xform[0][2] = axis.x * axis.z * invcosa + axis.y * sina; 4.181 + xform[1][0] = axis.x * axis.y * invcosa + axis.z * sina; 4.182 + xform[1][1] = nysq + (1-nysq) * cosa; 4.183 + xform[1][2] = axis.y * axis.z * invcosa - axis.x * sina; 4.184 + xform[2][0] = axis.x * axis.z * invcosa - axis.y * sina; 4.185 + xform[2][1] = axis.y * axis.z * invcosa + axis.x * sina; 4.186 + xform[2][2] = nzsq + (1-nzsq) * cosa; 4.187 + 4.188 + *this *= Matrix4x4(xform); 4.189 +} 4.190 + 4.191 +void Matrix4x4::set_rotation(const Vector3 &axis, float angle) 4.192 +{ 4.193 + float sina = (float)sin(angle); 4.194 + float cosa = (float)cos(angle); 4.195 + float invcosa = 1-cosa; 4.196 + float nxsq = axis.x * axis.x; 4.197 + float nysq = axis.y * axis.y; 4.198 + float nzsq = axis.z * axis.z; 4.199 + 4.200 + reset_identity(); 4.201 + m[0][0] = nxsq + (1-nxsq) * cosa; 4.202 + m[0][1] = axis.x * axis.y * invcosa - axis.z * sina; 4.203 + m[0][2] = axis.x * axis.z * invcosa + axis.y * sina; 4.204 + m[1][0] = axis.x * axis.y * invcosa + axis.z * sina; 4.205 + m[1][1] = nysq + (1-nysq) * cosa; 4.206 + m[1][2] = axis.y * axis.z * invcosa - axis.x * sina; 4.207 + m[2][0] = axis.x * axis.z * invcosa - axis.y * sina; 4.208 + m[2][1] = axis.y * axis.z * invcosa + axis.x * sina; 4.209 + m[2][2] = nzsq + (1-nzsq) * cosa; 4.210 +} 4.211 + 4.212 +void Matrix4x4::scale(const Vector3 &scale_vec) 4.213 +{ 4.214 + Matrix4x4 smat( scale_vec.x, 0, 0, 0, 4.215 + 0, scale_vec.y, 0, 0, 4.216 + 0, 0, scale_vec.z, 0, 4.217 + 0, 0, 0, 1.0); 4.218 + *this *= smat; 4.219 +} 4.220 + 4.221 +void Matrix4x4::set_scaling(const Vector3 &scale_vec) 4.222 +{ 4.223 + *this = Matrix4x4( scale_vec.x, 0, 0, 0, 4.224 + 0, scale_vec.y, 0, 0, 4.225 + 0, 0, scale_vec.z, 0, 4.226 + 0, 0, 0, 1); 4.227 +} 4.228 + 4.229 +void Matrix4x4::set_column_vector(const Vector3 &vec, unsigned int col_index) 4.230 +{ 4.231 + m[0][col_index] = vec.x; 4.232 + m[1][col_index] = vec.y; 4.233 + m[2][col_index] = vec.z; 4.234 + m[3][col_index] = 1; 4.235 +} 4.236 + 4.237 +void Matrix4x4::set_row_vector(const Vector3 &vec, unsigned int row_index) 4.238 +{ 4.239 + m[row_index][0] = vec.x; 4.240 + m[row_index][1] = vec.y; 4.241 + m[row_index][2] = vec.z; 4.242 + m[row_index][3] = 1; 4.243 +} 4.244 + 4.245 +Vector3 Matrix4x4::get_column_vector(unsigned int col_index) const 4.246 +{ 4.247 + return Vector3(m[0][col_index], m[1][col_index], m[2][col_index]); 4.248 +} 4.249 + 4.250 +Vector3 Matrix4x4::get_row_vector(unsigned int row_index) const 4.251 +{ 4.252 + return Vector3(m[row_index][0], m[row_index][1], m[row_index][2]); 4.253 +} 4.254 + 4.255 +void Matrix4x4::transpose() 4.256 +{ 4.257 + Matrix4x4 tmp = *this; 4.258 + for(int i=0; i<4; i++) { 4.259 + for(int j=0; j<4; j++) { 4.260 + m[i][j] = tmp[j][i]; 4.261 + } 4.262 + } 4.263 +} 4.264 + 4.265 +Matrix4x4 Matrix4x4::transposed() const 4.266 +{ 4.267 + Matrix4x4 res; 4.268 + for(int i=0; i<4; i++) { 4.269 + for(int j=0; j<4; j++) { 4.270 + res[i][j] = m[j][i]; 4.271 + } 4.272 + } 4.273 + return res; 4.274 +} 4.275 + 4.276 +float Matrix4x4::determinant() const 4.277 +{ 4.278 + float det11 = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - 4.279 + (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + 4.280 + (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); 4.281 + 4.282 + float det12 = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - 4.283 + (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + 4.284 + (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); 4.285 + 4.286 + float det13 = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - 4.287 + (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + 4.288 + (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); 4.289 + 4.290 + float det14 = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - 4.291 + (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + 4.292 + (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); 4.293 + 4.294 + return m[0][0] * det11 - m[0][1] * det12 + m[0][2] * det13 - m[0][3] * det14; 4.295 +} 4.296 + 4.297 + 4.298 +Matrix4x4 Matrix4x4::adjoint() const 4.299 +{ 4.300 + Matrix4x4 coef; 4.301 + 4.302 + coef.m[0][0] = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - 4.303 + (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + 4.304 + (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); 4.305 + coef.m[0][1] = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - 4.306 + (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + 4.307 + (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); 4.308 + coef.m[0][2] = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - 4.309 + (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + 4.310 + (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); 4.311 + coef.m[0][3] = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - 4.312 + (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + 4.313 + (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); 4.314 + 4.315 + coef.m[1][0] = (m[0][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - 4.316 + (m[0][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + 4.317 + (m[0][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); 4.318 + coef.m[1][1] = (m[0][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - 4.319 + (m[0][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + 4.320 + (m[0][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); 4.321 + coef.m[1][2] = (m[0][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - 4.322 + (m[0][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + 4.323 + (m[0][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); 4.324 + coef.m[1][3] = (m[0][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - 4.325 + (m[0][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + 4.326 + (m[0][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); 4.327 + 4.328 + coef.m[2][0] = (m[0][1] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) - 4.329 + (m[0][2] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) + 4.330 + (m[0][3] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])); 4.331 + coef.m[2][1] = (m[0][0] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) - 4.332 + (m[0][2] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) + 4.333 + (m[0][3] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])); 4.334 + coef.m[2][2] = (m[0][0] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) - 4.335 + (m[0][1] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) + 4.336 + (m[0][3] * (m[1][0] * m[3][1] - m[3][0] * m[1][1])); 4.337 + coef.m[2][3] = (m[0][0] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])) - 4.338 + (m[0][1] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])) + 4.339 + (m[0][2] * (m[1][0] * m[3][1] - m[3][0] * m[1][1])); 4.340 + 4.341 + coef.m[3][0] = (m[0][1] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) - 4.342 + (m[0][2] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) + 4.343 + (m[0][3] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])); 4.344 + coef.m[3][1] = (m[0][0] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) - 4.345 + (m[0][2] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) + 4.346 + (m[0][3] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])); 4.347 + coef.m[3][2] = (m[0][0] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) - 4.348 + (m[0][1] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) + 4.349 + (m[0][3] * (m[1][0] * m[2][1] - m[2][0] * m[1][1])); 4.350 + coef.m[3][3] = (m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])) - 4.351 + (m[0][1] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])) + 4.352 + (m[0][2] * (m[1][0] * m[2][1] - m[2][0] * m[1][1])); 4.353 + 4.354 + coef.transpose(); 4.355 + 4.356 + for(int i=0; i<4; i++) { 4.357 + for(int j=0; j<4; j++) { 4.358 + coef.m[i][j] = j%2 ? -coef.m[i][j] : coef.m[i][j]; 4.359 + if(i%2) coef.m[i][j] = -coef.m[i][j]; 4.360 + } 4.361 + } 4.362 + 4.363 + return coef; 4.364 +} 4.365 + 4.366 +Matrix4x4 Matrix4x4::inverse() const 4.367 +{ 4.368 + Matrix4x4 adj = adjoint(); 4.369 + 4.370 + return adj * (1.0f / determinant()); 4.371 +} 4.372 + 4.373 +inline float *Matrix4x4::operator [](int index) 4.374 +{ 4.375 + return m[index]; 4.376 +} 4.377 + 4.378 +inline const float *Matrix4x4::operator [](int index) const 4.379 +{ 4.380 + return m[index]; 4.381 +} 4.382 + 4.383 +inline void Matrix4x4::reset_identity() 4.384 +{ 4.385 + *this = identity; 4.386 +} 4.387 + 4.388 +std::ostream &operator <<(std::ostream &out, const Matrix4x4 &mat) 4.389 +{ 4.390 + for(int i=0; i<4; i++) { 4.391 + char str[100]; 4.392 + sprintf(str, "[ %12.5f %12.5f %12.5f %12.5f ]\n", (float)mat.m[i][0], (float)mat.m[i][1], (float)mat.m[i][2], (float)mat.m[i][3]); 4.393 + out << str; 4.394 + } 4.395 + return out; 4.396 +}
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 5.2 +++ b/src/matrix.h Sat Jan 11 17:22:36 2014 +0200 5.3 @@ -0,0 +1,81 @@ 5.4 +#ifndef MATRIX_H_ 5.5 +#define MATRIX_H_ 5.6 + 5.7 +#include <iostream> 5.8 +#include "vector.h" 5.9 + 5.10 +class Matrix4x4 { 5.11 +public: 5.12 + float m[4][4]; 5.13 + 5.14 + static Matrix4x4 identity; 5.15 + 5.16 + Matrix4x4(); 5.17 + Matrix4x4( float m11, float m12, float m13, float m14, 5.18 + float m21, float m22, float m23, float m24, 5.19 + float m31, float m32, float m33, float m34, 5.20 + float m41, float m42, float m43, float m44); 5.21 + 5.22 + /* binary operations matrix (op) matrix */ 5.23 + friend Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2); 5.24 + friend Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2); 5.25 + friend Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2); 5.26 + 5.27 + friend void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2); 5.28 + friend void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2); 5.29 + friend inline void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2); 5.30 + 5.31 + /* binary operations matrix (op) scalar and scalar (op) matrix */ 5.32 + friend Matrix4x4 operator *(const Matrix4x4 &mat, float scalar); 5.33 + friend Matrix4x4 operator *(float scalar, const Matrix4x4 &mat); 5.34 + 5.35 + friend void operator *=(Matrix4x4 &mat, float scalar); 5.36 + 5.37 + inline float *operator [](int index); 5.38 + inline const float *operator [](int index) const; 5.39 + 5.40 + inline void reset_identity(); 5.41 + 5.42 + void translate(const Vector3 &trans); 5.43 + void set_translation(const Vector3 &trans); 5.44 + 5.45 + void rotate(const Vector3 &euler); /* 3d rotation with euler angles */ 5.46 + void rotate(const Vector3 &axis, float angle); /* 3d axis/angle rotation */ 5.47 + void set_rotation(const Vector3 &euler_angles); 5.48 + void set_rotation(const Vector3 &axis, float angle); 5.49 + 5.50 + void scale(const Vector3 &scale_vec); 5.51 + void set_scaling(const Vector3 &scale_vec); 5.52 + 5.53 + void set_column_vector(const Vector3 &vec, unsigned int col_index); 5.54 + void set_row_vector(const Vector3 &vec, unsigned int row_index); 5.55 + Vector3 get_column_vector(unsigned int col_index) const; 5.56 + Vector3 get_row_vector(unsigned int row_index) const; 5.57 + 5.58 + void transpose(); 5.59 + Matrix4x4 transposed() const; 5.60 + float determinant() const; 5.61 + Matrix4x4 adjoint() const; 5.62 + Matrix4x4 inverse() const; 5.63 +}; 5.64 + 5.65 +/* binary operations matrix (op) matrix */ 5.66 +Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2); 5.67 +Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2); 5.68 +Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2); 5.69 + 5.70 +void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2); 5.71 +void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2); 5.72 +inline void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2); 5.73 + 5.74 +/* binary operations matrix (op) scalar and scalar (op) matrix */ 5.75 +Matrix4x4 operator *(const Matrix4x4 &mat, float scalar); 5.76 +Matrix4x4 operator *(float scalar, const Matrix4x4 &mat); 5.77 + 5.78 +void operator *=(Matrix4x4 &mat, float scalar); 5.79 + 5.80 +std::ostream &operator <<(std::ostream &out, const Matrix4x4 &mat); 5.81 + 5.82 + 5.83 +#endif 5.84 +
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 6.2 +++ b/src/vector.cc Sat Jan 11 17:22:36 2014 +0200 6.3 @@ -0,0 +1,213 @@ 6.4 +#include <math.h> 6.5 +#include <stdio.h> 6.6 +#include "vector.h" 6.7 +#include "matrix.h" 6.8 + 6.9 +Vector3::Vector3() { 6.10 + x = 0; 6.11 + y = 0; 6.12 + z = 0; 6.13 +} 6.14 + 6.15 +Vector3::Vector3(float x, float y, float z) { 6.16 + this->x = x; 6.17 + this->y = y; 6.18 + this->z = z; 6.19 +} 6.20 + 6.21 +void Vector3::transform(const Matrix4x4 &tm) { 6.22 + float x1 = tm.m[0][0]*x + tm.m[0][1]*y + tm.m[0][2]*z + tm.m[0][3]; 6.23 + float y1 = tm.m[1][0]*x + tm.m[1][1]*y + tm.m[1][2]*z + tm.m[1][3]; 6.24 + float z1 = tm.m[2][0]*x + tm.m[2][1]*y + tm.m[2][2]*z + tm.m[2][3]; 6.25 + x = x1; 6.26 + y = y1; 6.27 + z = z1; 6.28 +} 6.29 + 6.30 +void Vector3::printv() { 6.31 + printf("%f\t%f\t%f\n", x, y, z); 6.32 +} 6.33 + 6.34 + 6.35 +bool operator < (const Vector3 &a, const Vector3 &b) { 6.36 + return a.x < b.x && a.y < b.y && a.z < b.z; 6.37 +} 6.38 + 6.39 +bool operator > (const Vector3 &a, const Vector3 &b) { 6.40 + return a.x > b.x && a.y > b.y && a.z > b.z; 6.41 +} 6.42 + 6.43 +bool operator == (const Vector3 &a, const Vector3 &b) { 6.44 + return a.x == b.x && a.y == b.y && a.z == b.z; 6.45 +} 6.46 + 6.47 +Vector3 operator + (const Vector3 &a, const Vector3 &b) { 6.48 + return Vector3(a.x + b.x, a.y + b.y, a.z + b.z); 6.49 +} 6.50 + 6.51 +Vector3 operator - (const Vector3 &a, const Vector3 &b) { 6.52 + return Vector3(a.x - b.x, a.y - b.y, a.z - b.z); 6.53 +} 6.54 + 6.55 +Vector3 operator - (const Vector3 &a) { 6.56 + return Vector3(-a.x, -a.y, -a.z); 6.57 +} 6.58 + 6.59 +Vector3 operator * (const Vector3 &a, const Vector3 &b) { 6.60 + return Vector3(a.x * b.x, a.y * b.y, a.z * b.z); 6.61 +} 6.62 + 6.63 +Vector3 operator * (const Vector3 &a, float b) { 6.64 + return Vector3(a.x*b, a.y*b, a.z*b); 6.65 +} 6.66 + 6.67 +Vector3 operator * (float b, const Vector3 &a) { 6.68 + return Vector3(a.x*b, a.y*b, a.z*b); 6.69 +} 6.70 + 6.71 +Vector3 operator / (const Vector3 &a, float b) { 6.72 + return Vector3(a.x / b, a.y / b, a.z / b); 6.73 +} 6.74 + 6.75 +const Vector3 &operator += (Vector3 &a, const Vector3 &b) { 6.76 + a.x += b.x; 6.77 + a.y += b.y; 6.78 + a.z += b.z; 6.79 + return a; 6.80 +} 6.81 + 6.82 +float length(const Vector3 &a) { 6.83 + return sqrt(a.x*a.x + a.y*a.y + a.z*a.z); 6.84 +} 6.85 + 6.86 +float dot(const Vector3 &a, const Vector3 &b) { 6.87 + return a.x*b.x + a.y*b.y + a.z*b.z; 6.88 +} 6.89 + 6.90 +Vector3 cross(const Vector3 &a, const Vector3 &b) { 6.91 + return Vector3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); 6.92 +} 6.93 + 6.94 +Vector3 normalize(const Vector3 &vec) { 6.95 + float mag = sqrt(vec.x*vec.x + vec.y*vec.y + vec.z*vec.z); 6.96 + return vec / mag; 6.97 +} 6.98 + 6.99 +Vector3 reflect(const Vector3 &v, const Vector3 &n) { 6.100 + return 2.0 * dot(v, n) * n - v; 6.101 +} 6.102 + 6.103 +////// 6.104 +// 6.105 +///// 6.106 + 6.107 +Vector4::Vector4() { 6.108 + x = 0; 6.109 + y = 0; 6.110 + z = 0; 6.111 + w = 0; 6.112 +} 6.113 + 6.114 +Vector4::Vector4(float x, float y, float z, float w) { 6.115 + this->x = x; 6.116 + this->y = y; 6.117 + this->z = z; 6.118 + this->w = w; 6.119 +} 6.120 + 6.121 +void Vector4::transform(const Matrix4x4 &tm) { 6.122 + float x1 = tm.m[0][0]*x + tm.m[0][1]*y + tm.m[0][2]*z + tm.m[0][3]; 6.123 + float y1 = tm.m[1][0]*x + tm.m[1][1]*y + tm.m[1][2]*z + tm.m[1][3]; 6.124 + float z1 = tm.m[2][0]*x + tm.m[2][1]*y + tm.m[2][2]*z + tm.m[2][3]; 6.125 + float w1 = tm.m[3][0]*x + tm.m[3][1]*y + tm.m[3][2]*z + tm.m[3][3]; 6.126 + x = x1; 6.127 + y = y1; 6.128 + z = z1; 6.129 + w = w1; 6.130 +} 6.131 + 6.132 +void Vector4::printv() { 6.133 + printf("%f\t%f\t%f\t%f\n", x, y, z, w); 6.134 +} 6.135 + 6.136 + 6.137 +bool operator < (const Vector4 &a, const Vector4 &b) { 6.138 + return a.x < b.x && a.y < b.y && a.z < b.z && a.w < b.w; 6.139 +} 6.140 + 6.141 +bool operator > (const Vector4 &a, const Vector4 &b) { 6.142 + return a.x > b.x && a.y > b.y && a.z > b.z && a.w > b.w; 6.143 +} 6.144 + 6.145 +bool operator == (const Vector4 &a, const Vector4 &b) { 6.146 + return a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w; 6.147 +} 6.148 + 6.149 +Vector4 operator + (const Vector4 &a, const Vector4 &b) { 6.150 + return Vector4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w); 6.151 +} 6.152 + 6.153 +Vector4 operator - (const Vector4 &a, const Vector4 &b) { 6.154 + return Vector4(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w); 6.155 +} 6.156 + 6.157 +Vector4 operator - (const Vector4 &a) { 6.158 + return Vector4(-a.x, -a.y, -a.z, -a.w); 6.159 +} 6.160 + 6.161 +Vector4 operator * (const Vector4 &a, const Vector4 &b) { 6.162 + return Vector4(a.x * b.x, a.y * b.y, a.z * b.z, a.w * b.w); 6.163 +} 6.164 + 6.165 +Vector4 operator * (const Vector4 &a, float b) { 6.166 + return Vector4(a.x*b, a.y*b, a.z*b, a.w*b); 6.167 +} 6.168 + 6.169 +Vector4 operator * (float b, const Vector4 &a) { 6.170 + return Vector4(a.x * b, a.y * b, a.z * b, a.w * b); 6.171 +} 6.172 + 6.173 +Vector4 operator / (const Vector4 &a, float b) { 6.174 + return Vector4(a.x / b, a.y / b, a.z / b, a.w / b); 6.175 +} 6.176 + 6.177 +const Vector4 &operator += (Vector4 &a, const Vector4 &b) { 6.178 + a.x += b.x; 6.179 + a.y += b.y; 6.180 + a.z += b.z; 6.181 + a.w += b.w; 6.182 + return a; 6.183 +} 6.184 + 6.185 +float length(const Vector4 &a) { 6.186 + return sqrt(a.x*a.x + a.y*a.y + a.z*a.z + a.w*a.w); 6.187 +} 6.188 + 6.189 +float dot(const Vector4 &a, const Vector4 &b) { 6.190 + return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; 6.191 +} 6.192 + 6.193 +Vector4 cross(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3) { 6.194 + float a = (v2.x * v3.y) - (v2.y * v3.x); 6.195 + float b = (v2.x * v3.z) - (v2.z * v3.x); 6.196 + float c = (v2.x * v3.w) - (v2.w * v3.x); 6.197 + float d = (v2.y * v3.z) - (v2.z * v3.y); 6.198 + float e = (v2.y * v3.w) - (v2.w * v3.y); 6.199 + float f = (v2.z * v3.w) - (v2.w * v3.z); 6.200 + 6.201 + Vector4 result; 6.202 + result.x = (v1.y * f) - (v1.z * e) + (v1.w * d); 6.203 + result.y = - (v1.x * f) + (v1.z * c) - (v1.w * b); 6.204 + result.z = (v1.x * e) - (v1.y * c) + (v1.w * a); 6.205 + result.w = - (v1.x * d) + (v1.y * b) - (v1.z * a); 6.206 + return result; 6.207 +} 6.208 + 6.209 +Vector4 normalize(const Vector4 &vec) { 6.210 + float mag = sqrt(vec.x * vec.x + vec.y * vec.y + vec.z * vec.z + vec.w * vec.w); 6.211 + return vec / mag; 6.212 +} 6.213 + 6.214 +Vector4 reflect(const Vector4 &v, const Vector4 &n) { 6.215 + return 2.0 * dot(v, n) * n - v; 6.216 +}
7.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 7.2 +++ b/src/vector.h Sat Jan 11 17:22:36 2014 +0200 7.3 @@ -0,0 +1,68 @@ 7.4 +#ifndef VECTOR_H_ 7.5 +#define VECTOR_H_ 7.6 + 7.7 +class Matrix4x4; 7.8 + 7.9 +class Vector3 { 7.10 +public: 7.11 + float x,y,z; 7.12 + Vector3(); 7.13 + Vector3(float x, float y, float z); 7.14 + void transform(const Matrix4x4 &tm); 7.15 + void printv(); 7.16 +}; 7.17 + 7.18 +bool operator < (const Vector3 &a, const Vector3 &b); 7.19 +bool operator > (const Vector3 &a, const Vector3 &b); 7.20 +bool operator == (const Vector3 &a, const Vector3 &b); 7.21 + 7.22 +Vector3 operator + (const Vector3 &a, const Vector3 &b); 7.23 +Vector3 operator - (const Vector3 &a, const Vector3 &b); 7.24 +Vector3 operator - (const Vector3 &a); 7.25 +Vector3 operator * (const Vector3 &a, const Vector3 &b); 7.26 +Vector3 operator * (const Vector3 &a, float b); 7.27 +Vector3 operator * (float b, const Vector3 &a); 7.28 +Vector3 operator / (const Vector3 &a, float b); 7.29 + 7.30 +const Vector3 &operator += (Vector3 &a, const Vector3 &b); 7.31 + 7.32 +float length (const Vector3 &a); 7.33 +float dot (const Vector3 &a, const Vector3 &b); 7.34 +Vector3 cross (const Vector3 &a, const Vector3 &b); 7.35 +Vector3 normalize (const Vector3 &a); 7.36 + 7.37 +Vector3 reflect(const Vector3 &v, const Vector3 &n); 7.38 + 7.39 +class Vector4 { 7.40 + public: 7.41 + float x, y, z, w; 7.42 + Vector4(); 7.43 + Vector4(float x, float y, float z, float w); 7.44 + void transform(const Matrix4x4 &tm); 7.45 + void printv(); 7.46 +}; 7.47 + 7.48 +bool operator < (const Vector4 &a, const Vector4 &b); 7.49 +bool operator > (const Vector4 &a, const Vector4 &b); 7.50 +bool operator == (const Vector4 &a, const Vector4 &b); 7.51 + 7.52 +Vector4 operator + (const Vector4 &a, const Vector4 &b); 7.53 +Vector4 operator - (const Vector4 &a, const Vector4 &b); 7.54 +Vector4 operator - (const Vector4 &a); 7.55 +Vector4 operator * (const Vector4 &a, const Vector4 &b); 7.56 +Vector4 operator * (const Vector4 &a, float b); 7.57 +Vector4 operator * (float b, const Vector4 &a); 7.58 +Vector4 operator / (const Vector4 &a, float b); 7.59 + 7.60 +const Vector4 &operator += (Vector4 &a, const Vector4 &b); 7.61 + 7.62 +float length (const Vector4 &a); 7.63 +float dot (const Vector4 &a, const Vector4 &b); 7.64 +Vector4 cross (const Vector4 &v1, const Vector4 &v2, const Vector4 &v3); 7.65 +Vector4 normalize (const Vector4 &a); 7.66 + 7.67 +Vector4 reflect(const Vector4 &v, const Vector4 &n); 7.68 + 7.69 + 7.70 +#endif 7.71 +
8.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 8.2 +++ b/src/volume.cc Sat Jan 11 17:22:36 2014 +0200 8.3 @@ -0,0 +1,65 @@ 8.4 +#include <imago2.h> 8.5 + 8.6 +#include "image.h" 8.7 + 8.8 +Image::Image() 8.9 +{ 8.10 + pixels = 0; 8.11 + width = 0; 8.12 + height = 0; 8.13 +} 8.14 + 8.15 +Image::~Image() 8.16 +{ 8.17 + delete [] pixels; 8.18 +} 8.19 + 8.20 +bool Image::load(const char *fname) 8.21 +{ 8.22 + int new_width, new_height; 8.23 + float *new_pixels = (float*)img_load_pixels(fname, &new_width, &new_height, IMG_FMT_GREYF); 8.24 + 8.25 + if(!new_pixels) { 8.26 + fprintf(stderr, "Failed to load image: %s\n", fname); 8.27 + return false; 8.28 + } 8.29 + 8.30 + set_pixels(new_pixels, new_width, new_height); 8.31 + img_free_pixels(new_pixels); 8.32 + return true; 8.33 +} 8.34 + 8.35 +float *Image::get_pixels() 8.36 +{ 8.37 + return pixels; 8.38 +} 8.39 + 8.40 +const float *Image::get_pixels() const 8.41 +{ 8.42 + return pixels; 8.43 +} 8.44 + 8.45 +void Image::set_pixels(const float *pixels, int width, int height) 8.46 +{ 8.47 + if(!pixels) 8.48 + return; 8.49 + 8.50 + delete [] this->pixels; 8.51 + 8.52 + this->pixels = new float[width * height]; 8.53 + for(int i=0; i<width * height; i++) { 8.54 + this->pixels[i] = pixels[i]; 8.55 + } 8.56 + this->width = width; 8.57 + this->height = height; 8.58 +} 8.59 + 8.60 +int Image::get_width() const 8.61 +{ 8.62 + return width; 8.63 +} 8.64 + 8.65 +int Image::get_height() const 8.66 +{ 8.67 + return height; 8.68 +}