1 /* EXTRAITS DE LA LICENCE
2 Copyright CEA, contributeurs : Luc BILLARD et Damien
3 CALISTE, laboratoire L_Sim, (2001-2005)
4
5 Adresse mèl :
6 BILLARD, non joignable par mèl ;
7 CALISTE, damien P caliste AT cea P fr.
8
9 Ce logiciel est un programme informatique servant à visualiser des
10 structures atomiques dans un rendu pseudo-3D.
11
12 Ce logiciel est régi par la licence CeCILL soumise au droit français et
13 respectant les principes de diffusion des logiciels libres. Vous pouvez
14 utiliser, modifier et/ou redistribuer ce programme sous les conditions
15 de la licence CeCILL telle que diffusée par le CEA, le CNRS et l'INRIA
16 sur le site "http://www.cecill.info".
17
18 Le fait que vous puissiez accéder à cet en-tête signifie que vous avez
19 pris connaissance de la licence CeCILL, et que vous en avez accepté les
20 termes (cf. le fichier Documentation/licence.fr.txt fourni avec ce logiciel).
21 */
22
23 /* LICENCE SUM UP
24 Copyright CEA, contributors : Luc BILLARD et Damien
25 CALISTE, laboratoire L_Sim, (2001-2005)
26
27 E-mail address:
28 BILLARD, not reachable any more ;
29 CALISTE, damien P caliste AT cea P fr.
30
31 This software is a computer program whose purpose is to visualize atomic
32 configurations in 3D.
33
34 This software is governed by the CeCILL license under French law and
35 abiding by the rules of distribution of free software. You can use,
36 modify and/ or redistribute the software under the terms of the CeCILL
37 license as circulated by CEA, CNRS and INRIA at the following URL
38 "http://www.cecill.info".
39
40 The fact that you are presently reading this means that you have had
41 knowledge of the CeCILL license and that you accept its terms. You can
42 find a copy of this licence shipped with this software at Documentation/licence.en.txt.
43 */
44
45 #include "toolMatrix.h"
46
47 #include <visu_tools.h>
48 #include <math.h>
49 #include <visu_configFile.h>
50 #include "toolConfigFile.h"
51
52 /**
53 * SECTION:toolMatrix
54 * @short_description: Defines basic handlings on matrix.
55 *
56 * <para>Some very basic linear algebra are redefined here. It also
57 * gives access to coordinates conversion, essentially between
58 * cartesian and spherical.</para>
59 */
60
61 /**
62 * VisuBoxVertices:
63 * @vertices: (array fixed-size=8) (element-type ToolVector):
64 *
65 * Structure used for bindings.
66 *
67 * Since: 3.7
68 */
69
70 /**
71 * VisuBoxCell:
72 * @box: (array fixed-size=6) (element-type gdouble):
73 *
74 * Structure used for bindings.
75 *
76 * Since: 3.7
77 */
78
79 /**
80 * ToolVector:
81 *
82 * An opaque object, to exchange three floats as #GObject properties.
83 *
84 * Since: 3.8
85 */
86
87 /**
88 * ToolGridSize:
89 * @grid: (array fixed-size=3) (element-type guint):
90 *
91 * Structure used for bindings.
92 *
93 * Since: 3.7
94 */
95
96 /**
97 * tool_matrix_reducePrimitiveVectors:
98 * @reduced: (out caller-allocates) (array fixed-size=6): a storage for 6 floating point values ;
99 * @full: (in) (array fixed-size=9): a full 3x3 matrix to be transformed.
100 *
101 * This routine transforms the given matrix @full into a reduced array
102 * used by V_Sim to store box definition.
103 *
104 * Returns: FALSE if the given matrix is planar.
105 */
tool_matrix_reducePrimitiveVectors(double reduced[6],double full[3][3])106 gboolean tool_matrix_reducePrimitiveVectors(double reduced[6], double full[3][3])
107 {
108 double X[3];
109 double Y[3];
110 double Z[3];
111 double u[3], x[3];
112 int i, j, k;
113 double deltaIJ;
114 double norm;
115
116 g_return_val_if_fail(reduced && full, FALSE);
117
118 DBG_fprintf(stderr, "Matrix: transform full to reduced matrix.\n");
119 DBG_fprintf(stderr, "Matrix: full is %8g %8g %8g\n"
120 " %8g %8g %8g\n"
121 " %8g %8g %8g\n",
122 full[0][0], full[0][1], full[0][2],
123 full[1][0], full[1][1], full[1][2],
124 full[2][0], full[2][1], full[2][2]);
125 /* Compute the X vector of the new basis, colinear with old x. */
126 for (i = 0; i < 3; i++)
127 {
128 X[i] = full[0][i];
129 x[i] = full[0][i];
130 }
131
132 /* Compute the Y vector of the new basis, orthogonal to X and
133 coplanar with X and old y vector. */
134 u[0] = full[0][1] * full[1][2] - full[0][2] * full[1][1];
135 u[1] = full[0][2] * full[1][0] - full[0][0] * full[1][2];
136 u[2] = full[0][0] * full[1][1] - full[0][1] * full[1][0];
137 /* DBG_fprintf(stderr, "x : %f %f %f\n", x[0], x[1], x[2]); */
138 /* DBG_fprintf(stderr, "x vect y : %f %f %f\n", u[0], u[1], u[2]); */
139 deltaIJ = x[0] * u[1] - x[1] * u[0];
140 if (deltaIJ != 0.)
141 {
142 i = 0;
143 j = 1;
144 k = 2;
145 DBG_fprintf(stderr, " Using deltaIJ scheme with (i, j, k)"
146 " = (%d, %d, %d)\n", i, j, k);
147 }
148 else
149 {
150 deltaIJ = x[0] * u[2] - x[2] * u[0];
151 if (deltaIJ != 0.)
152 {
153 i = 0;
154 j = 2;
155 k = 1;
156 DBG_fprintf(stderr, " Using deltaIJ scheme with (i, j, k)"
157 " = (%d, %d, %d)\n", i, j, k);
158 }
159 else
160 {
161 deltaIJ = x[1] * u[2] - x[2] * u[1];
162 if (deltaIJ != 0.)
163 {
164 i = 1;
165 j = 2;
166 k = 0;
167 DBG_fprintf(stderr, " Using deltaIJ scheme with (i, j, k)"
168 " = (%d, %d, %d)\n", i, j, k);
169 }
170 else
171 return FALSE;
172 }
173 }
174 Y[k] = -1.;
175 Y[i] = (x[k] * u[j] - x[j] * u[k]) / deltaIJ;
176 Y[j] = (x[i] * u[k] - x[k] * u[i]) / deltaIJ;
177 /* We need to turn Y if y.Y is negative. */
178 norm = 0.;
179 for (i = 0; i < 3; i++)
180 norm += full[1][i] * Y[i];
181 if (norm < 0.)
182 for (i = 0; i < 3; i++)
183 Y[i] *= -1.;
184
185 /* Compute the new Z vector in order to form a direct orthogonal
186 basis with X and Y. */
187 Z[0] = X[1] * Y[2] - X[2] * Y[1];
188 Z[1] = X[2] * Y[0] - X[0] * Y[2];
189 Z[2] = X[0] * Y[1] - X[1] * Y[0];
190
191 /* Normalise the new basis (X, Y, Z). */
192 norm = 0.;
193 for (i = 0; i < 3; i++)
194 norm += X[i] * X[i];
195 norm = sqrt(norm);
196 for (i = 0; i < 3; i++)
197 X[i] /= norm;
198 norm = 0.;
199 for (i = 0; i < 3; i++)
200 norm += Y[i] * Y[i];
201 norm = sqrt(norm);
202 for (i = 0; i < 3; i++)
203 Y[i] /= norm;
204 norm = 0.;
205 for (i = 0; i < 3; i++)
206 norm += Z[i] * Z[i];
207 norm = sqrt(norm);
208 for (i = 0; i < 3; i++)
209 Z[i] /= norm;
210
211 /* DBG_fprintf(stderr, "X : %f %f %f\n", X[0], X[1], X[2]); */
212 /* DBG_fprintf(stderr, "Y : %f %f %f\n", Y[0], Y[1], Y[2]); */
213 /* DBG_fprintf(stderr, "Z : %f %f %f\n", Z[0], Z[1], Z[2]); */
214
215 /* Compute the reduce value for the basis. */
216 DBG_fprintf(stderr, " Write to reduced (%p).\n", (gpointer)reduced);
217 reduced[0] = 0.;
218 for (i = 0; i < 3; i++)
219 reduced[0] += X[i] * full[0][i];
220
221 reduced[1] = 0.;
222 for (i = 0; i < 3; i++)
223 reduced[1] += X[i] * full[1][i];
224
225 reduced[2] = 0.;
226 for (i = 0; i < 3; i++)
227 reduced[2] += Y[i] * full[1][i];
228
229 reduced[3] = 0.;
230 for (i = 0; i < 3; i++)
231 reduced[3] += X[i] * full[2][i];
232
233 reduced[4] = 0.;
234 for (i = 0; i < 3; i++)
235 reduced[4] += Y[i] * full[2][i];
236
237 reduced[5] = 0.;
238 for (i = 0; i < 3; i++)
239 reduced[5] += Z[i] * full[2][i];
240 DBG_fprintf(stderr, " Write OK.\n");
241
242 return TRUE;
243 }
244
245 /**
246 * tool_matrix_dtof:
247 * @mf: a matrix in single precision.
248 * @md: a matrix in double precision.
249 *
250 * Cast @md into @mf.
251 *
252 * Since: 3.7
253 **/
tool_matrix_dtof(float mf[3][3],double md[3][3])254 void tool_matrix_dtof(float mf[3][3], double md[3][3])
255 {
256 int i, j;
257
258 for (i = 0; i < 3; i++)
259 for (j = 0; j < 3; j++)
260 mf[i][j] = md[i][j];
261 }
262
263 /**
264 * tool_matrix_setIdentity:
265 * @mat: (array fixed-size=9): a matrix location.
266 *
267 * Initialise @mat with the identity.
268 *
269 * Since: 3.7
270 **/
tool_matrix_setIdentity(float mat[3][3])271 void tool_matrix_setIdentity(float mat[3][3])
272 {
273 mat[0][0] = 1.f;
274 mat[0][1] = 0.f;
275 mat[0][2] = 0.f;
276 mat[1][0] = 0.f;
277 mat[1][1] = 1.f;
278 mat[1][2] = 0.f;
279 mat[2][0] = 0.f;
280 mat[2][1] = 0.f;
281 mat[2][2] = 1.f;
282 }
283 /**
284 * tool_matrix_set:
285 * @mat: a matrix.
286 * @orig: a matrix.
287 *
288 * Copy @orig into @mat.
289 *
290 * Since: 3.8
291 **/
tool_matrix_set(float mat[3][3],float orig[3][3])292 void tool_matrix_set(float mat[3][3], float orig[3][3])
293 {
294 mat[0][0] = orig[0][0];
295 mat[0][1] = orig[0][1];
296 mat[0][2] = orig[0][2];
297 mat[1][0] = orig[1][0];
298 mat[1][1] = orig[1][1];
299 mat[1][2] = orig[1][2];
300 mat[2][0] = orig[2][0];
301 mat[2][1] = orig[2][1];
302 mat[2][2] = orig[2][2];
303 }
304 /**
305 * tool_matrix_rotate:
306 * @mat: a matrix.
307 * @angle: an angle in degrees.
308 * @dir: a direction.
309 *
310 * Create a rotation matrix along axis @dir of the given @angle.
311 *
312 * Since: 3.8
313 **/
tool_matrix_rotate(float mat[3][3],float angle,ToolXyzDir dir)314 void tool_matrix_rotate(float mat[3][3], float angle, ToolXyzDir dir)
315 {
316 float rot[3][3], work[3][3];
317 int perm[3][3] = {{1, 2, 0}, {2, 0, 1}, {0, 1, 2}};
318
319 rot[perm[dir][0]][perm[dir][0]] = cos(angle * TOOL_PI180);
320 rot[perm[dir][0]][perm[dir][1]] = -sin(angle * TOOL_PI180);
321 rot[perm[dir][0]][perm[dir][2]] = 0.f;
322 rot[perm[dir][1]][perm[dir][0]] = sin(angle * TOOL_PI180);
323 rot[perm[dir][1]][perm[dir][1]] = cos(angle * TOOL_PI180);
324 rot[perm[dir][1]][perm[dir][2]] = 0.f;
325 rot[perm[dir][2]][perm[dir][0]] = 0.f;
326 rot[perm[dir][2]][perm[dir][1]] = 0.f;
327 rot[perm[dir][2]][perm[dir][2]] = 1.f;
328
329 tool_matrix_set(work, mat);
330 tool_matrix_productMatrix(mat, work, rot);
331 }
332
333 /**
334 * tool_matrix_productMatrix:
335 * @matRes: an array of floating point values of size 3x3 ;
336 * @matA: an array of floating point values of size 3x3 ;
337 * @matB: an array of floating point values of size 3x3.
338 *
339 * Compute the mathematical product between @matA and @matB and
340 * put the result matrix in @matRes.
341 *
342 * Since: 3.2
343 */
tool_matrix_productMatrix(float matRes[3][3],float matA[3][3],float matB[3][3])344 void tool_matrix_productMatrix(float matRes[3][3], float matA[3][3], float matB[3][3])
345 {
346 int i, j, k;
347
348 for (i = 0; i < 3; i++)
349 for (j = 0; j < 3; j++)
350 {
351 matRes[i][j] = 0.;
352 for (k = 0; k < 3; k++)
353 matRes[i][j] += matA[i][k] * matB[k][j];
354 }
355 }
356 /**
357 * tool_matrix_productVector:
358 * @vectRes: an array of floating point values of size 3 ;
359 * @mat: an array of floating point values of size 3x3 ;
360 * @vect: an array of floating point values of size 3.
361 *
362 * Compute the mathematical product between @matA and @vect and
363 * put the result vector in @vectRes.
364 *
365 * Since: 3.2
366 */
tool_matrix_productVector(float vectRes[3],float mat[3][3],float vect[3])367 void tool_matrix_productVector(float vectRes[3], float mat[3][3], float vect[3])
368 {
369 int i, j;
370
371 for (i = 0; i < 3; i++)
372 {
373 vectRes[i] = 0.;
374 for (j = 0; j < 3; j++)
375 vectRes[i] += mat[i][j] * vect[j];
376 }
377 }
378 /**
379 * tool_matrix_determinant:
380 * @mat: a matrix.
381 *
382 * Calculate the determinant of matrix @mat.
383 *
384 * Since: 3.6
385 *
386 * Returns: the determinant value.
387 */
tool_matrix_determinant(float mat[3][3])388 float tool_matrix_determinant(float mat[3][3])
389 {
390 DBG_fprintf(stderr, "Tool Matrix: %g is det( %8g %8g %8g\n"
391 " %8g %8g %8g\n"
392 " %8g %8g %8g )\n",
393 mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
394 mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
395 mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]),
396 mat[0][0], mat[0][1], mat[0][2],
397 mat[1][0], mat[1][1], mat[1][2],
398 mat[2][0], mat[2][1], mat[2][2]);
399
400 return mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) -
401 mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) +
402 mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
403 }
404 /**
405 * tool_matrix_invert:
406 * @inv: a matrix (out values).
407 * @mat: a matrix.
408 *
409 * Calculate the inverse matrix of matrix @mat and store it in @inv.
410 *
411 * Since: 3.6
412 *
413 * Returns: FALSE if @mat is singular.
414 */
tool_matrix_invert(float inv[3][3],float mat[3][3])415 gboolean tool_matrix_invert(float inv[3][3], float mat[3][3])
416 {
417 float det;
418
419 det = tool_matrix_determinant(mat);
420 if (det == 0.f)
421 return FALSE;
422
423 det = 1.f / det;
424
425 inv[0][0] = det * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]);
426 inv[0][1] = det * (mat[0][2] * mat[2][1] - mat[0][1] * mat[2][2]);
427 inv[0][2] = det * (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]);
428
429 inv[1][0] = det * (mat[1][2] * mat[2][0] - mat[1][0] * mat[2][2]);
430 inv[1][1] = det * (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]);
431 inv[1][2] = det * (mat[0][2] * mat[1][0] - mat[0][0] * mat[1][2]);
432
433 inv[2][0] = det * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
434 inv[2][1] = det * (mat[0][1] * mat[2][0] - mat[0][0] * mat[2][1]);
435 inv[2][2] = det * (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
436
437 return TRUE;
438 }
439 /**
440 * tool_matrix_getRotationFromFull:
441 * @rot: a rotation matrix (out values).
442 * @full: the description of basis set in full development.
443 * @box: the description of basis set in align X axis.
444 *
445 * There is a rotation matrix to transform from full cartesian
446 * coordinates into reduced box cartesian coordinates.
447 *
448 * Since: 3.6
449 *
450 * Returns: TRUE if @full does not describe properly a 3D box.
451 */
tool_matrix_getRotationFromFull(float rot[3][3],double full[3][3],double box[6])452 gboolean tool_matrix_getRotationFromFull(float rot[3][3],
453 double full[3][3], double box[6])
454 {
455 float boxMat[3][3], fileMat[3][3], fileMatInv[3][3];
456
457 /* We create the rotation matrix that pass from the cartesian
458 coordinates of the file to the cartesian coordinates of the box. */
459 boxMat[0][0] = box[0];
460 boxMat[0][1] = box[1];
461 boxMat[0][2] = box[3];
462 boxMat[1][0] = 0.f;
463 boxMat[1][1] = box[2];
464 boxMat[1][2] = box[4];
465 boxMat[2][0] = 0.f;
466 boxMat[2][1] = 0.f;
467 boxMat[2][2] = box[5];
468 fileMat[0][0] = (float)full[0][0];
469 fileMat[0][1] = (float)full[1][0];
470 fileMat[0][2] = (float)full[2][0];
471 fileMat[1][0] = (float)full[0][1];
472 fileMat[1][1] = (float)full[1][1];
473 fileMat[1][2] = (float)full[2][1];
474 fileMat[2][0] = (float)full[0][2];
475 fileMat[2][1] = (float)full[1][2];
476 fileMat[2][2] = (float)full[2][2];
477 if (!tool_matrix_invert(fileMatInv, fileMat))
478 return FALSE;
479
480 tool_matrix_productMatrix(rot, boxMat, fileMatInv);
481 DBG_fprintf(stderr, "Tool Matrix: rotation matrix %8g %8g %8g\n"
482 " %8g %8g %8g\n"
483 " %8g %8g %8g\n",
484 rot[0][0], rot[0][1], rot[0][2],
485 rot[1][0], rot[1][1], rot[1][2],
486 rot[2][0], rot[2][1], rot[2][2]);
487 return TRUE;
488 }
489
490 #ifndef RAD2DEG
491 #define RAD2DEG(x) (57.29577951308232311 * x)
492 #endif
493
494 #ifndef DEG2RAD
495 #define DEG2RAD(x) (0.01745329251994329509 * x)
496 #endif
497
498 /**
499 * tool_matrix_cartesianToSpherical:
500 * @spherical: an allocated array of 3 floating point values to store the result ;
501 * @cartesian: an allocated array of 3 floating point values to read the input.
502 *
503 * A method to transform cartesian coordinates in spherical
504 * coordinates (radius, phi and theta).
505 *
506 * Since: 3.3
507 */
tool_matrix_cartesianToSpherical(float spherical[3],const float cartesian[3])508 void tool_matrix_cartesianToSpherical(float spherical[3], const float cartesian[3])
509 {
510 /* s[0] = rho, s[1] = theta, s[2] = phi
511 c[0] = x, c[1] = y, c[2] = z */
512 double rho;
513 double theta;
514 double phi;
515
516 if (cartesian[0] == 0 && cartesian[1] == 0 && cartesian[2] == 0)
517 {
518 spherical[TOOL_MATRIX_SPHERICAL_MODULUS] = 0;
519 spherical[TOOL_MATRIX_SPHERICAL_THETA] = 0;
520 spherical[TOOL_MATRIX_SPHERICAL_PHI] = 0;
521 return;
522 }
523
524 rho = sqrt(cartesian[0] * cartesian[0] +
525 cartesian[1] * cartesian[1] +
526 cartesian[2] * cartesian[2]);
527
528 if (cartesian[0] == 0 && cartesian[1] == 0)
529 theta = (cartesian[2] > 0) ? 0 : G_PI;
530 else
531 theta = acos(CLAMP(cartesian[2] / rho, -1., 1.));
532
533 if (cartesian[0] != 0)
534 {
535 phi = atan(cartesian[1] / cartesian[0]) + G_PI*((cartesian[0] < 0) ? 1 : 0);
536 }
537 else
538 {
539 if (cartesian[1] == 0) /* facultatif*/
540 phi = 0; /* facultatif*/
541 else if(cartesian[1] > 0)
542 phi = G_PI_2;
543 else
544 phi = -G_PI_2;
545 }
546
547 spherical[TOOL_MATRIX_SPHERICAL_MODULUS] = rho;
548 spherical[TOOL_MATRIX_SPHERICAL_THETA] = RAD2DEG(theta);
549 spherical[TOOL_MATRIX_SPHERICAL_PHI] = tool_modulo_float(RAD2DEG(phi), 360);
550 }
551
552 /**
553 * tool_matrix_sphericalToCartesian:
554 * @cartesian: an allocated array of 3 floating point values to store the result ;
555 * @spherical: an allocated array of 3 floating point values to read the input.
556 *
557 * A method to transform spherical coordinates (radius, phi and theta)
558 * to cartesian coordinates.
559 *
560 * Since: 3.3
561 */
tool_matrix_sphericalToCartesian(float cartesian[3],const float spherical[3])562 void tool_matrix_sphericalToCartesian(float cartesian[3], const float spherical[3])
563 {
564 cartesian[0] = spherical[TOOL_MATRIX_SPHERICAL_MODULUS] *
565 sin(DEG2RAD(spherical[TOOL_MATRIX_SPHERICAL_THETA])) *
566 cos(DEG2RAD(spherical[TOOL_MATRIX_SPHERICAL_PHI]));
567 cartesian[1] = spherical[TOOL_MATRIX_SPHERICAL_MODULUS] *
568 sin(DEG2RAD(spherical[TOOL_MATRIX_SPHERICAL_THETA])) *
569 sin(DEG2RAD(spherical[TOOL_MATRIX_SPHERICAL_PHI]));
570 cartesian[2] = spherical[TOOL_MATRIX_SPHERICAL_MODULUS] *
571 cos(DEG2RAD(spherical[TOOL_MATRIX_SPHERICAL_THETA]));
572 }
573
574 /**
575 * tool_matrix_getInter2D:
576 * @lambda: a location to store a float.
577 * @a: a point.
578 * @b: another point.
579 * @A: a point.
580 * @B: another point.
581 *
582 * Get the intersection coeeficient of lines [ab] and [AB].
583 *
584 * Returns: TRUE if [ab] and [AB] have an intersection.
585 */
tool_matrix_getInter2D(float * lambda,float a[2],float b[2],float A[2],float B[2])586 gboolean tool_matrix_getInter2D(float *lambda,
587 float a[2], float b[2], float A[2], float B[2])
588 {
589 float denom;
590
591 denom = (b[0] - a[0]) * (B[1] - A[1]) - (b[1] - a[1]) * (B[0] - A[0]);
592 if (denom == 0.f)
593 return FALSE;
594 *lambda = (A[0] - a[0]) * (B[1] - A[1]) - (A[1] - a[1]) * (B[0] - A[0]);
595 *lambda /= denom;
596 /* fprintf(stderr, "%g\n", *lambda); */
597 return TRUE;
598 }
599 /**
600 * tool_matrix_getInter2DFromList: (skip)
601 * @i: a location to store a point.
602 * @lambda: a location to store a float.
603 * @a: a point.
604 * @b: another point.
605 * @set: a list of points.
606 *
607 * Same as tool_matrix_getInter2D(), but from a list of points.
608 *
609 * Returns: TRUE if an intersection exists.
610 */
tool_matrix_getInter2DFromList(float i[2],float * lambda,float a[2],float b[2],GList * set)611 gboolean tool_matrix_getInter2DFromList(float i[2], float *lambda,
612 float a[2], float b[2], GList *set)
613 {
614 float *pt1, *pt2;
615 float l, min;
616
617 i[0] = a[0];
618 i[1] = a[1];
619
620 min = 1.2f;
621 for (pt1 = (float*)(g_list_last(set)->data); set; set = g_list_next(set))
622 {
623 pt2 = (float*)set->data;
624 if (tool_matrix_getInter2D(&l, a, b, pt1, pt2))
625 min = (l >= 0.f)?MIN(min, l):min;
626 pt1 = pt2;
627 }
628 if (min > 1.00001f)
629 return FALSE;
630 if (lambda)
631 *lambda = min;
632 i[0] = (b[0] - a[0]) * min + a[0];
633 i[1] = (b[1] - a[1]) * min + a[1];
634 /* fprintf(stderr, "%g -> %gx%g\n", min, i[0], i[1]); */
635 return TRUE;
636 }
637
638 #define FLAG_PARAMETER_THRESHOLD "scale_log_threshold"
639 #define DESC_PARAMETER_THRESHOLD "Value of the threshold used in the zero centred TOOL_MATRIX_SCALING_LOG scaling function ; a positive float (1e-3)"
640 static float threshold = 1e-3;
641 static void exportParameters(GString *data, VisuData *dataObj);
642
643 /**
644 * tool_matrix_getScaledLinear:
645 * @x: the initial value ;
646 * @minmax: the boundaries for the @x argument ;
647 *
648 * Transform @x into [0;1] with a linear scale.
649 *
650 * Returns: a value into [0;1].
651 *
652 * Since: 3.5
653 */
tool_matrix_getScaledLinear(float x,float minmax[2])654 float tool_matrix_getScaledLinear(float x, float minmax[2])
655 {
656 return ((CLAMP(x, minmax[0], minmax[1]) - minmax[0]) /
657 (minmax[1] - minmax[0]));
658 }
659 /**
660 * tool_matrix_getScaledLog:
661 * @x: the initial value ;
662 * @minmax: the boundaries for the @x argument.
663 *
664 * Transform @x into [0;1] with a log scale.
665 *
666 * Returns: a value into [0;1].
667 *
668 * Since: 3.5
669 */
tool_matrix_getScaledLog(float x,float minmax[2])670 float tool_matrix_getScaledLog(float x, float minmax[2])
671 {
672 /* float v; */
673 float lMinMax[2];
674
675 lMinMax[0] = log10(MAX(1e-12, minmax[0]));
676 lMinMax[1] = log10(MAX(1e-12, minmax[1]));
677 return tool_matrix_getScaledLinear(log10(MAX(1e-12, x)), lMinMax);
678 /* return (v == 0.)?0.: - (log10(v) - param) / param; */
679 }
680 /**
681 * tool_matrix_getScaledZeroCentredLog:
682 * @x: the initial value ;
683 * @minmax: the boundaries for the @x argument.
684 *
685 * Transform @x into [0;1] with a log scale with zero centred values.
686 *
687 * Returns: a value into [0;1].
688 *
689 * Since: 3.5
690 */
tool_matrix_getScaledZeroCentredLog(float x,float minmax[2])691 float tool_matrix_getScaledZeroCentredLog(float x, float minmax[2])
692 {
693 float v, m;
694
695 m = MAX(minmax[1], -minmax[0]);
696 v = CLAMP(x, -m, m);
697 return 0.5 + (v < 0.?-1.:1.) * (log(m * threshold) -
698 log(MAX(ABS(v), m * threshold))) /
699 (2. * log(threshold));
700 }
701 /**
702 * tool_matrix_getScaledLinearInv:
703 * @x: the initial value ;
704 * @minmax: the boundaries for the @x argument.
705 *
706 * Reverse function for tool_matrix_getScaledLinear().
707 *
708 * Returns: a value into [0;1].
709 *
710 * Since: 3.5
711 */
tool_matrix_getScaledLinearInv(float x,float minmax[2])712 float tool_matrix_getScaledLinearInv(float x, float minmax[2])
713 {
714 return (minmax[0] + CLAMP(x, 0., 1.) * (minmax[1] - minmax[0]));
715 }
716 /**
717 * tool_matrix_getScaledLogInv:
718 * @x: the initial value ;
719 * @minmax: the boundaries for the @x argument.
720 *
721 * Reverse function for tool_matrix_getScaledLog().
722 *
723 * Returns: a value into [0;1].
724 *
725 * Since: 3.5
726 */
tool_matrix_getScaledLogInv(float x,float minmax[2])727 float tool_matrix_getScaledLogInv(float x, float minmax[2])
728 {
729 return MAX(1e-12, minmax[0]) * pow(MAX(1e-12, minmax[1]) / MAX(1e-12, minmax[0]), CLAMP(x, 0., 1.));
730 /* return (minmax[0] + (minmax[1] - minmax[0]) * exp((1. - CLAMP(x, 0., 1.)) * param)); */
731 }
732 /**
733 * tool_matrix_getScaledZeroCentredLogInv:
734 * @x: the initial value ;
735 * @minmax: the boundaries for the @x argument.
736 *
737 * Reverse function for tool_matrix_getScaledZeroCentredLog().
738 *
739 * Returns: a value into [0;1].
740 *
741 * Since: 3.5
742 */
tool_matrix_getScaledZeroCentredLogInv(float x,float minmax[2])743 float tool_matrix_getScaledZeroCentredLogInv(float x, float minmax[2])
744 {
745 float s, m, out;
746
747 DBG_fprintf(stderr, "Matrix: get inv ZCL %g (%g-%g).\n", x, minmax[0], minmax[1]);
748 s = (x < 0.5)?-1.:1.;
749 m = MAX(minmax[1], -minmax[0]);
750 out = s * m * threshold * exp(s * (1. - 2. * CLAMP(x, 0., 1.)) * log(threshold));
751 DBG_fprintf(stderr, " | %g (%g)\n", out, tool_matrix_getScaledZeroCentredLog(out, minmax));
752 return out;
753 }
754
exportParameters(GString * data,VisuData * dataObj _U_)755 static void exportParameters(GString *data, VisuData *dataObj _U_)
756 {
757 g_string_append_printf(data, "# %s\n", DESC_PARAMETER_THRESHOLD);
758 g_string_append_printf(data, "%s: %f\n\n", FLAG_PARAMETER_THRESHOLD,
759 threshold);
760 }
761 /**
762 * tool_matrix_init: (skip)
763 *
764 * This method is used by V_Sim internally and should not be called.
765 *
766 * Since: 3.5
767 */
tool_matrix_init(void)768 void tool_matrix_init(void)
769 {
770 float rg[2] = {G_MINFLOAT, G_MAXFLOAT};
771 VisuConfigFileEntry *entry;
772
773 /* Set private variables. */
774 entry = visu_config_file_addFloatArrayEntry(VISU_CONFIG_FILE_PARAMETER,
775 FLAG_PARAMETER_THRESHOLD,
776 DESC_PARAMETER_THRESHOLD,
777 1, &threshold, rg, FALSE);
778 visu_config_file_entry_setVersion(entry, 3.5f);
779 visu_config_file_addExportFunction(VISU_CONFIG_FILE_PARAMETER,
780 exportParameters);
781 }
782
vector_copy(gconstpointer boxed)783 static gpointer vector_copy(gconstpointer boxed)
784 {
785 gfloat *vector, *orig = (gfloat*)boxed;
786
787 if (!boxed)
788 return (gpointer)0;
789
790 vector = g_malloc(sizeof(gfloat) * 3);
791 vector[0] = orig[0];
792 vector[1] = orig[1];
793 vector[2] = orig[2];
794 return (gpointer)vector;
795 }
796 /**
797 * tool_vector_get_type:
798 *
799 * Create and retrieve a #GType for a #ToolVector object.
800 *
801 * Since: 3.8
802 *
803 * Returns: a new type for #ToolVector structures.
804 */
tool_vector_get_type(void)805 GType tool_vector_get_type(void)
806 {
807 static GType g_define_type_id = 0;
808
809 if (g_define_type_id == 0)
810 g_define_type_id = g_boxed_type_register_static("ToolVector",
811 (GBoxedCopyFunc)vector_copy, g_free);
812 return g_define_type_id;
813 }
814
815 /**
816 * tool_vector_new:
817 * @orig: (array fixed-size=3): a vector.
818 *
819 * Creates a new #ToolVector by copying @orig.
820 *
821 * Since: 3.8
822 *
823 * Returns: (transfer full) (type ToolVector):
824 **/
tool_vector_new(const float orig[3])825 float* tool_vector_new(const float orig[3])
826 {
827 return vector_copy(orig);
828 }
829
830 /**
831 * tool_vector_set:
832 * @dest: a vector.
833 * @orig: a vector.
834 *
835 * Copy @orig into @dest and also test if @dest was already equals to @orig.
836 *
837 * Since: 3.8
838 *
839 * Returns: TRUE if @dest is different from @orig.
840 **/
tool_vector_set(float dest[3],const float orig[3])841 gboolean tool_vector_set(float dest[3], const float orig[3])
842 {
843 if (dest[0] == orig[0] && dest[1] == orig[1] && dest[2] == orig[2])
844 return FALSE;
845
846 dest[0] = orig[0];
847 dest[1] = orig[1];
848 dest[2] = orig[2];
849 return TRUE;
850 }
851
852 /**
853 * tool_vector_spherical:
854 * @cart: (array fixed-size=3): a vector in cartesian coordinates.
855 * @at: a spherical direction.
856 *
857 * Convert @cart in spherical coordinates. It is equivalent to
858 * tool_matrix_cartesianToSpherical() but returns the value in a given
859 * direction only. This is intended for bindings.
860 *
861 * Since: 3.8
862 *
863 * Returns: the coodinate in spherical.
864 **/
tool_vector_spherical(const float cart[3],ToolMatrixSphericalCoord at)865 float tool_vector_spherical(const float cart[3], ToolMatrixSphericalCoord at)
866 {
867 float sph[3];
868
869 tool_matrix_cartesianToSpherical(sph, cart);
870 return sph[at];
871 }
872
minmax_copy(gpointer boxed)873 static gpointer minmax_copy(gpointer boxed)
874 {
875 gfloat *mm, *orig = (gfloat*)boxed;
876
877 if (!boxed)
878 return (gpointer)0;
879
880 mm = g_malloc(sizeof(gfloat) * 2);
881 mm[0] = orig[0];
882 mm[1] = orig[1];
883 return (gpointer)mm;
884 }
885 /**
886 * tool_minmax_get_type:
887 *
888 * Create and retrieve a #GType for a #ToolMinmax object.
889 *
890 * Since: 3.8
891 *
892 * Returns: a new type for #ToolMinmax structures.
893 */
tool_minmax_get_type(void)894 GType tool_minmax_get_type(void)
895 {
896 static GType g_define_type_id = 0;
897
898 if (g_define_type_id == 0)
899 g_define_type_id = g_boxed_type_register_static("ToolMinmax", minmax_copy, g_free);
900 return g_define_type_id;
901 }
902
903 /**
904 * tool_minmax:
905 * @global: (array fixed-size=2): the global min and max values.
906 * @minmax: (array fixed-size=2): some min and max values.
907 *
908 * Take the minimum and maximum of @global and @minmax and put them
909 * into @global. @global should be initialised.
910 *
911 * Since: 3.8
912 **/
tool_minmax(float global[2],const float minmax[2])913 void tool_minmax(float global[2], const float minmax[2])
914 {
915 global[0] = MIN(minmax[0], global[0]);
916 global[1] = MAX(minmax[1], global[1]);
917 }
918 /**
919 * tool_minmax_fromDbl:
920 * @global: (array fixed-size=2): the global min and max values.
921 * @minmax: (array fixed-size=2): some min and max values.
922 *
923 * Same as tool_minmax() for double inputs.
924 *
925 * Since: 3.8
926 **/
tool_minmax_fromDbl(float global[2],const double minmax[2])927 void tool_minmax_fromDbl(float global[2], const double minmax[2])
928 {
929 global[0] = MIN(minmax[0], global[0]);
930 global[1] = MAX(minmax[1], global[1]);
931 }
932