luc/sweet
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

add new functions

Browse Source
This commit is contained in:
Luc 2019-03-09 16:13:14 +01:00
parent 3ecdb295b5
commit c30b00c955
7 changed files with 350 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Makefile
View File

@ -10,15 +10,9 @@ OBJS_WITH_SHORT = sweet_math_short.o sweet_matrix_short.o
CFLAGS = -pedantic -fPIC -Wall CFLAGS = -pedantic -fPIC -Wall
LD = -shared LD = -shared
ifeq ($(TARGET),short)
all: $(PROG) all: $(PROG)
$(PROG): $(OBJS) $(OBJS_WITH_SHORT) $(PROG): $(OBJS) $(OBJS_WITH_SHORT)
$(CC) $(LD) $(OBJS) $(OBJS_WITH_SHORT) -o $(PROG) $(CC) $(LD) $(OBJS) $(OBJS_WITH_SHORT) -o $(PROG)
else
all: $(PROG)
$(PROG): $(OBJS)
$(CC) $(LD) $(OBJS) -o $(PROG)
endif
.c.o: .c.o:
$(CC) $(CFLAGS) -c $*.c $(CC) $(CFLAGS) -c $*.c

54
sweet_math.c
View File

@ -17,6 +17,7 @@
*/ */
#include <stdio.h> #include <stdio.h>
#include <stdlib.h>
#include <math.h> #include <math.h>
#include <stdarg.h> #include <stdarg.h>
#include "sweet_math.h" #include "sweet_math.h"
@ -78,10 +79,10 @@ sweet_math_quadratic_polynomial (vec2 * r, float a, float b, float c)
#define FRAC_1_ON_3 0.3333333333333333333333333333333333333333333333333333333333333333333 #define FRAC_1_ON_3 0.3333333333333333333333333333333333333333333333333333333333333333333
#define FRAC_2_ON_3 0.6666666666666666666666666666666666666666666666666666666666666666666 #define FRAC_2_ON_3 0.6666666666666666666666666666666666666666666666666666666666666666667
#define FRAC_1_ON_6 0.1666666666666666666666666666666666666666666666666666666666666666666 #define FRAC_1_ON_6 0.1666666666666666666666666666666666666666666666666666666666666666667
#define FRAC_1_ON_9 0.1111111111111111111111111111111111111111111111111111111111111111111 #define FRAC_1_ON_9 0.1111111111111111111111111111111111111111111111111111111111111111111
#define FRAC_2_ON_27 0.0740740740740740740740740740740740740740740740740740740740740740740 #define FRAC_2_ON_27 0.0740740740740740740740740740740740740740740740740740740740740740741
#define FRAC_1_ON_27 0.0370370370370370370370370370370370370370370370370370370370370370370 #define FRAC_1_ON_27 0.0370370370370370370370370370370370370370370370370370370370370370370
int int
@ -117,6 +118,8 @@ sweet_math_cubic_polynomial (vec3 * roots, float t, float a, float b, float c)
if (D < 0) if (D < 0)
{ {
roots->x = (r + s) - a_over_3; roots->x = (r + s) - a_over_3;
roots->y = -r - a_over_3;
roots->z = roots->y;
return 1; return 1;
} }
else if (sweet_math_approx_equals (D, 0, EPSILON)) else if (sweet_math_approx_equals (D, 0, EPSILON))
@ -486,6 +489,33 @@ sweet_vector_rescale4 (vec4 v, float size)
return v; return v;
} }
vec2
sweet_vector_inverse2 (vec2 v)
{
v.x = - v.x;
v.y = - v.y;
return v;
}
vec3
sweet_vector_inverse3 (vec3 v)
{
v.x = - v.x;
v.y = - v.y;
v.z = - v.z;
return v;
}
vec4
sweet_vector_inverse4 (vec4 v)
{
v.x = - v.x;
v.y = - v.y;
v.z = - v.z;
v.w = - v.w;
return v;
}
vec2 vec2
sweet_vector_add2 (vec2 v1, vec2 v2) sweet_vector_add2 (vec2 v1, vec2 v2)
{ {
@ -598,6 +628,24 @@ sweet_vector_angle3 (vec3 a, vec3 b)
return acos (dot); return acos (dot);
} }
vec3
sweet_vector_normalized_random ()
{
vec3 v;
v.x = (float)rand() / (float)RAND_MAX;
v.y = (float)rand() / (float)RAND_MAX;
v.z = (float)rand() / (float)RAND_MAX;
return sweet_vector_normalize3 (v);
}
vec3
sweet_vector_random_ortho(vec3 v)
{
vec3 u = sweet_vector_normalized_random ();
v = sweet_vector_cross(u, v);
return sweet_vector_normalize3(v);
}
vec2 vec2
sweet_vector_product2 (vec2 v1, vec2 v2) sweet_vector_product2 (vec2 v1, vec2 v2)
{ {

11
sweet_math.h
View File

@ -185,6 +185,13 @@ vec2 sweet_vector_rescale2 (vec2 vector, float size);
vec3 sweet_vector_rescale3 (vec3 vector, float size); vec3 sweet_vector_rescale3 (vec3 vector, float size);
vec4 sweet_vector_rescale4 (vec4 vector, float size); vec4 sweet_vector_rescale4 (vec4 vector, float size);
/** Inverse */
/** @param Vector as vector*/
/** @return - vector */
vec2 sweet_vector_inverse2 (vec2 vector);
vec3 sweet_vector_inverse3 (vec3 vector);
vec4 sweet_vector_inverse4 (vec4 vector);
/** Add */ /** Add */
/** @param Vector1 as vector */ /** @param Vector1 as vector */
/** @param Vector2 as vector */ /** @param Vector2 as vector */
@ -216,6 +223,10 @@ vec4 sweet_vector_sub4 (vec4 vector1, vec4 vector2);
float sweet_vector_angle2 (vec2 vector1, vec2 vector2); float sweet_vector_angle2 (vec2 vector1, vec2 vector2);
float sweet_vector_angle3 (vec3 vector1, vec3 vector2); float sweet_vector_angle3 (vec3 vector1, vec3 vector2);
/** Random */
vec3 sweet_vector_normalized_random ();
vec3 sweet_vector_random_ortho(vec3 v);
/** Mult */ /** Mult */
/** @param Vector1 as vector */ /** @param Vector1 as vector */
/** @param Vector2 as vector */ /** @param Vector2 as vector */

3
sweet_math_short.h
View File

@ -81,6 +81,9 @@ vec4 vec4_new_v4 (vec4 v);
#define vec2_rescale sweet_vector_rescale2 #define vec2_rescale sweet_vector_rescale2
#define vec3_rescale sweet_vector_rescale3 #define vec3_rescale sweet_vector_rescale3
#define vec4_rescale sweet_vector_rescale4 #define vec4_rescale sweet_vector_rescale4
#define vec2_inverse sweet_vector_inverse2
#define vec3_inverse sweet_vector_inverse3
#define vec4_inverse sweet_vector_inverse4
#define vec2_add sweet_vector_add2 #define vec2_add sweet_vector_add2
#define vec3_add sweet_vector_add3 #define vec3_add sweet_vector_add3
#define vec4_add sweet_vector_add4 #define vec4_add sweet_vector_add4

383
sweet_matrix.c
View File

@ -16,8 +16,6 @@
* License along with this library. If not, see <http://www.gnu.org/licenses/>. * License along with this library. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <string.h> #include <string.h>
#include <math.h> #include <math.h>
#include "sweet_matrix.h" #include "sweet_matrix.h"
@ -86,7 +84,6 @@ sweet_matrix_null4 (void)
return m; return m;
} }
mat2 mat2
sweet_matrix_identity2 (void) sweet_matrix_identity2 (void)
{ {
@ -149,13 +146,12 @@ sweet_matrix_identity4 (void)
return m; return m;
} }
mat2 mat2
sweet_matrix_ortho_basis2 (vec2 u, vec2 v) sweet_matrix_ortho_basis2 (vec2 u, vec2 v)
{ {
mat2 m; mat2 m;
m.v[0] = u.x; m.v[0] = u.x;
m.v[3] = v.x; m.v[1] = v.x;
m.v[2] = u.y; m.v[2] = u.y;
m.v[3] = v.y; m.v[3] = v.y;
@ -182,28 +178,55 @@ sweet_matrix_ortho_basis3 (vec3 u, vec3 v, vec3 w)
return m; return m;
} }
mat4
sweet_matrix_ortho_basis3h (vec3 p, vec3 u, vec3 v, vec3 w)
{
mat4 m;
m.v[0] = u.x;
m.v[1] = v.x;
m.v[2] = w.x;
m.v[3] = 0;
m.v[4] = u.y;
m.v[5] = v.y;
m.v[6] = w.y;
m.v[7] = 0;
m.v[8] = u.z;
m.v[9] = v.z;
m.v[10] = w.z;
m.v[11] = 0;
m.v[12] = -p.x * u.x - p.y * u.y - p.z * u.z;
m.v[13] = -p.x * v.x - p.y * v.y - p.z * v.z;
m.v[14] = -p.x * w.x - p.y * w.y - p.z * w.z;
m.v[15] = 1;
return m;
}
mat4 mat4
sweet_matrix_ortho_basis4 (vec4 u, vec4 v, vec4 w, vec4 t) sweet_matrix_ortho_basis4 (vec4 u, vec4 v, vec4 w, vec4 t)
{ {
mat4 m; mat4 m;
m.v[0] = u.x; m.v[0] = u.x;
m.v[1] = v.y; m.v[1] = v.x;
m.v[2] = w.z; m.v[2] = w.x;
m.v[3] = t.w; m.v[3] = t.x;
m.v[4] = u.x; m.v[4] = u.y;
m.v[5] = v.y; m.v[5] = v.y;
m.v[6] = w.z; m.v[6] = w.y;
m.v[7] = t.w; m.v[7] = t.y;
m.v[8] = u.x; m.v[8] = u.z;
m.v[9] = v.y; m.v[9] = v.z;
m.v[10] = w.z; m.v[10] = w.z;
m.v[11] = t.w; m.v[11] = t.z;
m.v[12] = u.x; m.v[12] = u.w;
m.v[13] = v.y; m.v[13] = v.w;
m.v[14] = w.z; m.v[14] = w.w;
m.v[15] = t.w; m.v[15] = t.w;
return m; return m;
@ -331,6 +354,14 @@ sweet_matrix_look_at (vec3 pos, vec3 dir, vec3 up)
return m; return m;
} }
mat4
sweet_matrix_look_at_ortho (vec3 pos, vec3 dir, vec3 up, vec3 right)
{
mat4 m;
dir = sweet_vector_scale3 (dir, -1);
return sweet_matrix_ortho_basis3h (pos, right, up, dir);
}
mat3 mat3
sweet_matrix_quat_rotation3 (quaternion q) sweet_matrix_quat_rotation3 (quaternion q)
{ {
@ -495,13 +526,13 @@ sweet_matrix_translation2h (vec2 v)
m.v[1] = 0.0; m.v[1] = 0.0;
m.v[2] = 0.0; m.v[2] = 0.0;
m.v[4] = 0.0; m.v[3] = 0.0;
m.v[5] = 1.0; m.v[4] = 1.0;
m.v[6] = 0.0; m.v[5] = 0.0;
m.v[8] = v.x; m.v[6] = v.x;
m.v[9] = v.y; m.v[7] = v.y;
m.v[10] = 1.0; m.v[8] = 1.0;
return m; return m;
} }
@ -987,6 +1018,62 @@ sweet_matrix_product4 (mat4 * p, mat4 * a, mat4 * b)
p->v[15] = av3 * bv12 + av7 * bv13 + av11 * bv14 + a->v[15] * b->v[15]; p->v[15] = av3 * bv12 + av7 * bv13 + av11 * bv14 + a->v[15] * b->v[15];
} }
void
sweet_matrix_add_diag3 (mat3 * out, mat3 * in, float x, float y, float z)
{
out->v[0] = in->v[0] + x;
out->v[1] = in->v[1];
out->v[2] = in->v[2];
out->v[3] = in->v[3];
out->v[4] = in->v[4] + y;
out->v[5] = in->v[5];
out->v[6] = in->v[6];
out->v[6] = in->v[7];
out->v[8] = in->v[8] + z;
}
void
sweet_matrix_substract_diag3 (mat3 * out, mat3 * in, float x, float y, float z)
{
out->v[0] = in->v[0] - x;
out->v[1] = in->v[1];
out->v[2] = in->v[2];
out->v[3] = in->v[3];
out->v[4] = in->v[4] - y;
out->v[5] = in->v[5];
out->v[6] = in->v[6];
out->v[6] = in->v[7];
out->v[8] = in->v[8] - z;
}
void
sweet_matrix_add3 (mat3 * out, mat3 * m1, mat3 * m2)
{
out->v[0] = m1->v[0] + m2->v[0];
out->v[1] = m1->v[1] + m2->v[1];
out->v[2] = m1->v[2] + m2->v[2];
out->v[3] = m1->v[3] + m2->v[3];
out->v[4] = m1->v[4] + m2->v[4];
out->v[5] = m1->v[5] + m2->v[5];
out->v[6] = m1->v[6] + m2->v[6];
out->v[7] = m1->v[7] + m2->v[7];
out->v[8] = m1->v[8] + m2->v[8];
}
void
sweet_matrix_substract3 (mat3 * out, mat3 * m1, mat3 * m2)
{
out->v[0] = m1->v[0] - m2->v[0];
out->v[1] = m1->v[1] - m2->v[1];
out->v[2] = m1->v[2] - m2->v[2];
out->v[3] = m1->v[3] - m2->v[3];
out->v[4] = m1->v[4] - m2->v[4];
out->v[5] = m1->v[5] - m2->v[5];
out->v[6] = m1->v[6] - m2->v[6];
out->v[7] = m1->v[7] - m2->v[7];
out->v[8] = m1->v[8] - m2->v[8];
}
void void
sweet_matrix_sub3 (mat2 * sub, mat3 * m) sweet_matrix_sub3 (mat2 * sub, mat3 * m)
{ {
@ -1089,18 +1176,31 @@ sweet_matrix_eigen_value3 (float * h1, float * h2, float * h3, mat3 * m)
float b; float b;
float c; float c;
int n; int n;
vec3 v; vec3 v = sweet_vector_new3 (0,0,0);
a = m->v[0] + m->v[4] + m->v[8]; a = m->v[0] + m->v[4] + m->v[8];
b = m->v[1] * m->v[3] + m->v[2] * m->v[6] + m->v[7] * m->v[5] - b = m->v[1] * m->v[3] + m->v[2] * m->v[6] + m->v[7] * m->v[5] -
m->v[0] * m->v[4] - m->v[0] * m->v[8] - m->v[4] * m->v[8]; m->v[0] * m->v[4] - m->v[0] * m->v[8] - m->v[4] * m->v[8];
c = m->v[0] * m->v[4] * m->v[8] - m->v[0] * m->v[7] * m->v[5] - c = m->v[0] * (m->v[4] * m->v[8] - m->v[5] * m->v[7]) -
m->v[1] * m->v[3] * m->v[8] + m->v[1] * m->v[6] * m->v[5] + m->v[1] * (m->v[3] * m->v[8] - m->v[5] * m->v[6]) +
m->v[2] * m->v[3] * m->v[7] - m->v[2] * m->v[6] * m->v[4]; m->v[2] * (m->v[3] * m->v[7] - m->v[4] * m->v[6]);
n = sweet_math_cubic_polynomial (&v, -1, a, b, c); n = sweet_math_cubic_polynomial (&v, -1, a, b, c);
if (v.y > v.x)
{
float t = v.x; v.x = v.y; v.y = t;
}
if (v.z > v.y)
{
float t = v.y; v.y = v.z; v.z = t;
}
if (v.y > v.x)
{
float t = v.x; v.x = v.y; v.y = t;
}
*h1 = v.x; *h1 = v.x;
*h2 = v.y; *h2 = v.y;
*h3 = v.z; *h3 = v.z;
@ -1205,7 +1305,7 @@ static void substract_row (mat3 * a, float * s, int l1, int l2, int i)
} }
#define EPSILON 0.00005 #define EPSILON 1e-05f
int int
sweet_matrix_gaussian_elimination3 (mat3 * a, vec3 * v) sweet_matrix_gaussian_elimination3 (mat3 * a, vec3 * v)
@ -1254,107 +1354,178 @@ sweet_matrix_solve3 (mat3 * m, vec3 * v)
return 0; return 0;
} }
static int is_diagonal (mat3 * m) static void build_Q0 (mat3 * q, mat3 * m)
{ {
if (!sweet_math_approx_equals (m->v[1], 0.0, EPSILON)) { return 0; } float a = sqrt (m->v[0] * m->v[0] + m->v[1] * m->v[1] + m->v[2] * m->v[2]);
if (!sweet_math_approx_equals (m->v[2], 0.0, EPSILON)) { return 0; } if (m->v[0] > 0) { a = - a; }
if (!sweet_math_approx_equals (m->v[3], 0.0, EPSILON)) { return 0; }
if (!sweet_math_approx_equals (m->v[5], 0.0, EPSILON)) { return 0; } vec3 v = sweet_vector_new3 (m->v[0] - a, m->v[1], m->v[2]);
if (!sweet_math_approx_equals (m->v[6], 0.0, EPSILON)) { return 0; } float norm = sweet_vector_norm3 (v);
if (!sweet_math_approx_equals (m->v[7], 0.0, EPSILON)) { return 0; } if (norm > 0)
return 1; v = sweet_vector_scale3 (v, 1.0 / norm);
q->v[0] = 1 - 2 * v.x * v.x;
q->v[1] = - 2 * v.y * v.x;
q->v[2] = - 2 * v.z * v.x;
q->v[3] = - 2 * v.x * v.y;
q->v[4] = 1 - 2 * v.y * v.y;
q->v[5] = - 2 * v.z * v.y;
q->v[6] = - 2 * v.x * v.z;
q->v[7] = - 2 * v.y * v.z;
q->v[8] = 1 - 2 * v.z * v.z;
} }
static int col_n_zero (mat3 * m, int i) static void build_Q1 (mat3 * q, mat3 * m)
{ {
i = i * 3; float a = sqrt (m->v[4] * m->v[4] + m->v[5] * m->v[5]);
return (sweet_math_approx_equals (m->v[i], 0, EPSILON) && if (m->v[4] > 0) { a = -a; }
sweet_math_approx_equals (m->v[i+1], 0, EPSILON) &&
sweet_math_approx_equals (m->v[i+2], 0, EPSILON)); vec2 v = sweet_vector_new2 (m->v[4] - a, m->v[5]);
float norm = sweet_vector_norm2 (v);
if (norm > 0)
v = sweet_vector_scale2 (v, 1.0 / norm);
q->v[0] = 1;
q->v[1] = 0;
q->v[2] = 0;
q->v[3] = 0;
q->v[4] = 1 - 2 * v.x * v.x;
q->v[5] = - 2 * v.y * v.x;
q->v[6] = 0;
q->v[7] = - 2 * v.x * v.y;
q->v[8] = 1 - 2 * v.y * v.y;
} }
static int row_n_zero (mat3 * m, int i) void sweet_matrix_QR3 (mat3 * Q, mat3 * R, mat3 * m)
{ {
return (sweet_math_approx_equals (m->v[i], 0, EPSILON) && mat3 tmp;
sweet_math_approx_equals (m->v[i+3], 0, EPSILON) && mat3 Q0;
sweet_math_approx_equals (m->v[i+6], 0, EPSILON)); mat3 Q1;
build_Q0 (&Q0, m);
sweet_matrix_product3 (&tmp, &Q0, m);
build_Q1 (&Q1, &tmp);
sweet_matrix_product3 (Q, &Q0, &Q1);
sweet_matrix_transpose3 (&tmp, Q);
sweet_matrix_product3 (R, &tmp, m);
} }
static void get_eigen_vector (vec3 * v, mat3 a, float h) void sweet_matrix_hess3 (mat3 * H, mat3 * Q, mat3 * m)
{ {
a.v[0] = a.v[0] - h; mat3 tmp;
a.v[4] = a.v[4] - h;
a.v[8] = a.v[8] - h;
v->x = 0; float a = sqrt (m->v[1] * m->v[1] + m->v[2] * m->v[2]);
v->y = 0; if (m->v[1] > 0)
v->z = 0; a = - a;
sweet_matrix_gaussian_elimination3 (&a, v); vec2 v = sweet_vector_new2 (m->v[1] - a, m->v[2]);
float norm = sweet_vector_norm2 (v);
if (norm > 0)
v = sweet_vector_scale2 (v, 1.0/norm);
if (col_n_zero (&a, 0)) Q->v[0] = 1;
v->x = 1; Q->v[1] = 0;
else if (col_n_zero(&a, 1)) Q->v[2] = 0;
v->y = 1; Q->v[3] = 0;
else if (col_n_zero(&a, 2)) Q->v[4] = 1 - 2 * v.x * v.x;
v->z = 1; Q->v[5] = - 2 * v.x * v.y;
else Q->v[6] = 0;
Q->v[7] = - 2 * v.x * v.y;
Q->v[8] = 1 - 2 * v.y * v.y;
sweet_matrix_product3 (&tmp, Q, m);
sweet_matrix_product3 (H, &tmp, Q);
}
static void shift (mat3 * m, float s)
{
m->v[0] = m->v[0] + s;
m->v[4] = m->v[4] + s;
m->v[8] = m->v[8] + s;
}
void print_matrix (char * title, mat3 * m)
{
printf ("%s:\n[%.10f %.10f %.10f]\n[%.10f %.10f %.10f]\n[%.10f %.10f %.10f]\n",
title,
m->v[0], m->v[3], m->v[6],
m->v[1], m->v[4], m->v[7],
m->v[2], m->v[5], m->v[8]);
}
static void sort3 (float * val, vec3 * vec)
{
if (val[1] > val[0])
{ {
if (sweet_math_approx_equals (a.v[2], 0, EPSILON) && float t = val[0];
sweet_math_approx_equals (a.v[5], 0, EPSILON)) vec3 v = vec[0];
{ val[0] = val[1];
if (sweet_math_approx_equals (a.v[8], 0, EPSILON)) val[1] = t;
{ vec[0] = vec[1];
v->z = 1; vec[1] = v;
v->y = - a.v[7];
v->x = (-v->y * a.v[3] - a.v[6]);
}
else
{
v->z = 0;
v->y = 1;
v->x = (-a.v[3]);
}
}
else if (row_n_zero (&a, 1))
{
v->y = 1;
v->z = 0;
v->x = (-v->y * a.v[3]);
}
} }
*v = sweet_vector_normalize3 (*v); if (val[2] > val[1])
{
float t = val[1];
vec3 v = vec[1];
val[1] = val[2];
val[2] = t;
vec[1] = vec[2];
vec[2] = v;
}
/*
if (val[1] > val[0])
{
float t = val[0];
vec3 v = vec[0];
val[0] = val[1];
val[1] = t;
vec[0] = vec[1];
vec[1] = v;
}*/
} }
int void sweet_matrix_eigen_vectors3 (vec3 * eigen_vectors, float * eigen_values, mat3 * A, int it)
sweet_matrix_eigen_vector3 (vec3 * vectors, mat3 * m)
{ {
/* (A - hI)V = 0 */ int i;
mat3 a; mat3 H, S, Q, W, R;
int n; sweet_matrix_hess3 (&H, &W, A);
float h1, h2, h3;
if (is_diagonal (m)) float s = H.v[8];
shift (&H, -H.v[8]);
sweet_matrix_QR3 (&Q, &R, &H);
sweet_matrix_product3 (&H, &R, &Q);
S = Q;
for (i = 0; i < it; i++)
{ {
vectors[0] = sweet_vector_new3 (1,0,0); sweet_matrix_QR3 (&Q, &R, &H);
vectors[1] = sweet_vector_new3 (0,1,0); sweet_matrix_product3 (&H, &R, &Q);
vectors[2] = sweet_vector_new3 (0,0,1); sweet_matrix_product3 (&S, &S, &Q);
return 3;
} }
n = sweet_matrix_eigen_value3 (&h1, &h2, &h3, m); shift (&H, s);
sweet_matrix_product3 (&S, &W, &S);
/* vector associated with h1 */ eigen_values[0] = H.v[0];
get_eigen_vector (&vectors[0], *m, h1); eigen_values[1] = H.v[4];
eigen_values[2] = H.v[8];
eigen_vectors[0] = sweet_vector_new3 (S.v[0], S.v[1], S.v[2]);
eigen_vectors[1] = sweet_vector_new3 (S.v[3], S.v[4], S.v[5]);
eigen_vectors[2] = sweet_vector_new3 (S.v[6], S.v[7], S.v[8]);
sort3 (eigen_values, eigen_vectors);
/* vector associated with h2 */ /*
get_eigen_vector (&vectors[1], *m, h2); printf ("%f %f %f\n", eigen_vectors[0].x, eigen_vectors[0].y, eigen_vectors[0].z);
printf ("%f %f %f\n", eigen_vectors[1].x, eigen_vectors[1].y, eigen_vectors[1].z);
printf ("%f %f %f\n", eigen_vectors[2].x, eigen_vectors[2].y, eigen_vectors[2].z);
/* vector associated with h3 */ printf ("%f %f %f\n", eigen_values[0], eigen_values[1], eigen_values[2]);
if (n > 2) */
get_eigen_vector (&vectors[2], *m, h3);
else
vectors[2] = sweet_vector_cross (vectors[0], vectors[1]);
return n;
} }

6
sweet_matrix.h
View File

@ -41,6 +41,7 @@ mat4 sweet_matrix_identity4 (void);
/** @return transform matrix to the new basis */ /** @return transform matrix to the new basis */
mat2 sweet_matrix_ortho_basis2 (vec2 u, vec2 v); mat2 sweet_matrix_ortho_basis2 (vec2 u, vec2 v);
mat3 sweet_matrix_ortho_basis3 (vec3 u, vec3 v, vec3 w); mat3 sweet_matrix_ortho_basis3 (vec3 u, vec3 v, vec3 w);
mat4 sweet_matrix_ortho_basis3h (vec3 p, vec3 u, vec3 v, vec3 w);
mat4 sweet_matrix_ortho_basis4 (vec4 u, vec4 v, vec4 w, vec4 t); mat4 sweet_matrix_ortho_basis4 (vec4 u, vec4 v, vec4 w, vec4 t);
/** Frustum */ /** Frustum */
@ -78,6 +79,8 @@ mat4 sweet_matrix_ortho (float left, float right,
/** @return camera matrix */ /** @return camera matrix */
mat4 sweet_matrix_look_at (vec3 pos, vec3 dir, vec3 up); mat4 sweet_matrix_look_at (vec3 pos, vec3 dir, vec3 up);
mat4 sweet_matrix_look_at_ortho (vec3 pos, vec3 dir, vec3 up, vec3 right);
/** Rotation from quaterion */ /** Rotation from quaterion */
/** @param quaternion representing rotation */ /** @param quaternion representing rotation */
/** @return 3d Rotation matrix */ /** @return 3d Rotation matrix */
@ -170,7 +173,8 @@ int sweet_matrix_gaussian_elimination3 (mat3 * m, vec3 * v);
int sweet_matrix_solve3 (mat3 * a, vec3 * v); int sweet_matrix_solve3 (mat3 * a, vec3 * v);
/** Eigen Vector */ /** Eigen Vector */
int sweet_matrix_eigen_vector3 (vec3 * vectors, mat3 * m); void sweet_matrix_eigen_vectors3 (vec3 * eigen_vectors, float * eigen_values, mat3 * A, int it);
#endif #endif

4
sweet_matrix_short.h
View File

@ -51,11 +51,13 @@ mat4 mat4_new_4v (vec4 e1, vec4 e2, vec4 e3, vec4 e4);
#define identity4 sweet_matrix_identity4 #define identity4 sweet_matrix_identity4
#define ortho_basis2 sweet_matrix_ortho_basis2 #define ortho_basis2 sweet_matrix_ortho_basis2
#define ortho_basis3 sweet_matrix_ortho_basis3 #define ortho_basis3 sweet_matrix_ortho_basis3
#define ortho_basis3h sweet_matrix_ortho_basis3h
#define ortho_basis4 sweet_matrix_ortho_basis4 #define ortho_basis4 sweet_matrix_ortho_basis4
#define matrix_frustum sweet_matrix_frustum #define matrix_frustum sweet_matrix_frustum
#define matrix_perspective sweet_matrix_perspective #define matrix_perspective sweet_matrix_perspective
#define matrix_ortho sweet_matrix_ortho #define matrix_ortho sweet_matrix_ortho
#define matrix_look_at sweet_matrix_look_at #define matrix_look_at sweet_matrix_look_at
#define matrix_look_at_ortho sweet_matrix_look_at_ortho
#define matrix_quat_rotation3 sweet_matrix_quat_rotation3 #define matrix_quat_rotation3 sweet_matrix_quat_rotation3
#define matrix_quat_rotation3h sweet_matrix_quat_rotation3h #define matrix_quat_rotation3h sweet_matrix_quat_rotation3h
#define rotation2 sweet_matrix_rotation2 #define rotation2 sweet_matrix_rotation2
@ -91,7 +93,7 @@ mat4 mat4_new_4v (vec4 e1, vec4 e2, vec4 e3, vec4 e4);
#define mat3_eigen_value sweet_matrix_eigen_value3 #define mat3_eigen_value sweet_matrix_eigen_value3
#define mat3_gaussian_elimination sweet_matrix_gaussian_elimination3 #define mat3_gaussian_elimination sweet_matrix_gaussian_elimination3
#define mat3_solve sweet_matrix_solve3 #define mat3_solve sweet_matrix_solve3
#define mat3_eigen_vector sweet_matrix_eigen_vector3 #define mat3_eigen_vector sweet_matrix_eigen_vectors3
#endif #endif