Merge branch 'roberto' into experimental

src/surface.c

@@ -1,6 +1,7 @@
 #include <complex.h>
 #include <math.h>
 #include <stdlib.h>
+#include <string.h>
 
 #define TEST
 
@@ -106,7 +107,7 @@ float *generate_data_surface(unsigned char *s)
 	int *cara;
 
 	parm.f = cube;
-	parm.m = 5;
+	parm.m = 6;
 	parm.n = parm.m;
 	parm.grid = 1;
 
@@ -183,79 +184,105 @@ float *generate_data_surface(unsigned char *s)
 }
 
 static void __calculate_normal(
-	float *p1, float *p2, float *p3, float *normal, unsigned char n)
+    float *p1, float *p2, float *p3, float *normal, unsigned char n)
 {
-	float alpha;
-	vec4 v1, v2, v3;
-	vec4 u1, u2, u3;
+    unsigned char i;
+    float alpha;
+    float *v1, *v2, *v3;
+    float *u1, *u2, *u3;
 
-	switch (n)
-	{
-	case 3:
+    v1 = malloc(n * sizeof(float));
+    v2 = malloc(n * sizeof(float));
+    v3 = malloc(n * sizeof(float));
+    u1 = malloc(n * sizeof(float));
+    u2 = malloc(n * sizeof(float));
+    u3 = malloc(n * sizeof(float));
 
-		glm_vec3_sub(p2, p1, v1);
-		glm_vec3_sub(p3, p1, v2);
+    /* 
+    	Calculate a normal vector of a plain using Gram-Schmidt process
+    */
+    {
+        for (i = 0; i < n; ++i) {
+            v1[i] = p2[i] - p1[i];  
+            v2[i] = p3[i] - p1[i]; 
+            v3[i] = p1[i];          
+        }
 
-		glm_vec3_cross(v1, v2, normal);
-		glm_vec3_normalize(normal);
-		return;
+        for (i = 0; i < n; ++i) {
+            u1[i] = v1[i]; 
+        }
 
-	case 4:
+        {
+            float proj[n];
+            float dot_v2_u1 = 0.0f, dot_u1_u1 = 0.0f;
+            for (i = 0; i < n; ++i) {
+                dot_v2_u1 += v2[i] * u1[i];
+                dot_u1_u1 += u1[i] * u1[i];
+            }
+            alpha = dot_v2_u1 / dot_u1_u1;
 
-		/*
-			In Grant-Shmidth we need 3 linearly independian vector that forms a
-		   basis, so we can have a ortonormal version of that basis, since, we
-		   must have v1 = p3 - p1 v2 = p2 - p1 Then v3 = p1, will most certantly
-		   be linerly independiant to v1 and v2.
-		*/
-		glm_vec4_sub(p2, p1, v1);
-		glm_vec4_sub(p3, p1, v2);
-		glm_vec4_copy(p1, v3);
+            for (i = 0; i < n; ++i) {
+                proj[i] = u1[i] * alpha;
+                u2[i] = v2[i] - proj[i]; 
+            }
+        }
 
-		/* Setup U1 */
-		{
-			glm_vec4_copy(v1, u1);
-		}
+        {
+            float proj1[n], proj2[n];
+            float dot_v3_u1 = 0.0f, dot_u1_u1 = 0.0f;
+            float dot_v3_u2 = 0.0f, dot_u2_u2 = 0.0f;
 
-		/* Setup U2 */
-		{
-			vec4 proj;
+            for (i = 0; i < n; ++i) {
+                dot_v3_u1 += v3[i] * u1[i];
+                dot_u1_u1 += u1[i] * u1[i];
+            }
+            for (i = 0; i < n; ++i) {
+                proj1[i] = u1[i] * (dot_v3_u1 / dot_u1_u1);
+            }
 
-			alpha = glm_vec4_dot(v2, u1) / glm_vec4_dot(u1, u1);
-			glm_vec4_scale(u1, alpha, proj);
+            for (i = 0; i < n; ++i) {
+                dot_v3_u2 += v3[i] * u2[i];
+                dot_u2_u2 += u2[i] * u2[i];
+            }
+            for (i = 0; i < n; ++i) {
+                proj2[i] = u2[i] * (dot_v3_u2 / dot_u2_u2);
+                u3[i] = v3[i] - proj1[i] - proj2[i]; 
+            }
+        }
 
-			glm_vec4_sub(v2, proj, u2);
-		}
+        float magnitude = 0.0f;
+        for (i = 0; i < n; ++i) {
+            magnitude += u3[i] * u3[i];
+        }
+        magnitude = sqrtf(magnitude);
 
-		/* Setup U3 */
-		{
-			vec4 proj1, proj2;
+        for (i = 0; i < n; ++i) {
+            normal[i] = u3[i] / magnitude;  
+        }
 
-			alpha = glm_vec4_dot(v3, u1) / glm_vec4_dot(u1, u1);
-			glm_vec4_scale(u1, alpha, proj1);
-
-			alpha = glm_vec4_dot(v3, u2) / glm_vec4_dot(u2, u2);
-			glm_vec4_scale(u2, alpha, proj2);
-
-			glm_vec4_sub(v3, proj1, u3);
-			glm_vec4_sub(u3, proj2, u3);
-		}
-
-		glm_vec4_copy(u3, normal);
-		glm_vec4_normalize(normal);
-		return;
-	}
+        free(v1);
+        free(v2);
+        free(v3);
+        free(u1);
+        free(u2);
+        free(u3);
+        return;
+    }
 }
 
+
 float *generate_normals_surface(float *d, unsigned char m)
 {
 	float *n;
 	n = malloc((*d + 1) * sizeof(float));
 	*n = *d;
+
+	float * norm_vec;
+	norm_vec=malloc(m*sizeof(float));
+
 	for (int i = 0; i < *d; i += 3 * m)
 	{
 
-		vec4 norm_vec;
 
 		__calculate_normal(
 			(d + 1) + i, (d + 1) + i + m, (d + 1) + i + 2 * m, norm_vec, m);
@@ -263,5 +290,7 @@ float *generate_normals_surface(float *d, unsigned char m)
 		glm_vec3_copy(norm_vec, (n + 1) + i + m);
 		glm_vec3_copy(norm_vec, (n + 1) + i + 2 * m);
 	}
+
+	free(norm_vec);
 	return n;
 }