Improve zoom accuracy.
[scilab.git] / scilab / modules / graphics / src / c / math_graphics.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 1998-2001 - ENPC - Jean-Philippe Chancelier
4  * Copyright (C) 2006 - INRIA - Jean-Baptiste Silvy
5  * 
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at    
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  */
13
14 /*------------------------------------------------------------------------
15  *    Graphic library
16  --------------------------------------------------------------------------*/
17
18 #include "math_graphics.h"
19
20 /* 
21  * we use spConfig.h for machine constants 
22  * XXX : spConfig should be merged and unified 
23  *       with other machine constant scilab code 
24  */
25
26 #define spINSIDE_SPARSE
27 #include "../../sparse/includes/spConfig.h"
28
29 double Mini(double *vect, integer n)
30 {
31   int i;
32   double vmin;
33   vmin = LARGEST_REAL;
34   for (i = 0 ; i < n ; i++)
35     /*    if ( isinf(vect[i])== 0 && isnan(vect[i])==0 && vect[i] < vmin)  */
36     if ( finite(vect[i])== 1 && vect[i] < vmin) 
37       vmin=vect[i];
38   return(vmin);
39 }
40
41 double Maxi(double *vect,integer n)
42 {
43   int i;
44   double maxi;
45   maxi= - LARGEST_REAL;
46   for (i =0 ; i < n ; i++)
47     /* if ( isinf(vect[i])== 0 && isnan(vect[i])==0 && vect[i] > maxi) */
48     if ( finite(vect[i])== 1 && vect[i] > maxi) 
49       maxi=vect[i];
50   return(maxi);
51 }
52
53 /*----------------------------------------------------------------------------*/
54
55 /* perform the rotation of point from to point dest  */
56 void rotate2D( double from[2], double center[2], double angle, double dest[2] )
57 {
58   rotate2Dim( from, center, cos( angle ), sin( angle ), dest ) ;
59 }
60
61 /*----------------------------------------------------------------------------*/
62 /* perform the rotation of point from to point to. */
63 /* the angle is directly given with its sine and cosine for speed */
64 void rotate2Dim( double from[2]   ,
65                  double center[2] ,
66                  double cosAngle  ,
67                  double sinAngle  ,
68                  double dest[2]    )
69 {
70   double diff[2] ;
71
72   /* put the center to (0,0) */
73   diff[0] = from[0] - center[0] ;
74   diff[1] = from[1] - center[1] ;
75
76   /* turn and translate back */
77   dest[0] = diff[0] * cosAngle - diff[1] * sinAngle + center[0] ;
78   dest[1] = diff[0] * sinAngle + diff[1] * cosAngle + center[1] ;
79 }
80 /*----------------------------------------------------------------------------*/
81 /* perform the rotation of point from to point dest given in int with angle in radian  */
82 void iRotate2D( int from[2], int center[2], double angle, int dest[2] )
83 {
84   iRotate2Dim( from, center, cos( angle ), sin( angle ), dest ) ;
85 }
86
87 /*----------------------------------------------------------------------------*/
88 /* perform the rotation of point from to point to. */
89 /* the angle is directly given with its sine and cosine for speed */
90 void iRotate2Dim( int    from[2]   ,
91                   int    center[2] ,
92                   double cosAngle  ,
93                   double sinAngle  ,
94                   int    dest[2]    )
95 {
96   int diff[2] ;
97
98   /* put the center to (0,0) */
99   diff[0] = from[0] - center[0] ;
100   diff[1] = from[1] - center[1] ;
101
102   /* turn and translate back */
103   dest[0] = round( diff[0] * cosAngle - diff[1] * sinAngle + center[0] ) ;
104   dest[1] = round( diff[0] * sinAngle + diff[1] * cosAngle + center[1] ) ;
105 }
106 /*----------------------------------------------------------------------------*/
107 /* perform an homethety point from to point dest. The 2 factors stand for the ration */
108 /* along the 2 coordinates */
109 void homothety2D( double from[2], double center[2], double factors[2], double dest[2] )
110 {
111   dest[0] = center[0] + factors[0] * ( from[0] - center[0] ) ;
112   dest[1] = center[1] + factors[1] * ( from[1] - center[1] ) ;
113 }
114 /*----------------------------------------------------------------------------*/
115 /* perform an homethety point from to point dest given in pixels. */
116 /* The 2 factors stand for the ration along the 2 coordinates */
117 void iHomothety2D( int from[2], int center[2], double factors[2], int dest[2] )
118 {
119   dest[0] = round( center[0] + factors[0] * ( from[0] - center[0] ) ) ;
120   dest[1] = round( center[1] + factors[1] * ( from[1] - center[1] ) ) ;
121 }
122 /*----------------------------------------------------------------------------*/
123 /* perform the translation of point from to point to with vector trans */
124 void translate2D( double from[2], double trans[2], double dest[2] )
125 {
126   dest[0] = from[0] + trans[0] ;
127   dest[1] = from[1] + trans[1] ;
128 }
129 /*----------------------------------------------------------------------------*/
130 void iTranslate2D( int from[2], int trans[2], int dest[2] )
131 {
132   dest[0] = from[0] + trans[0] ;
133   dest[1] = from[1] + trans[1] ;
134 }
135 /*----------------------------------------------------------------------------*/
136 void vectSubstract2D(const double vect1[2], const double vect2[], double res[2])
137 {
138   res[0] = vect1[0] - vect2[0];
139   res[1] = vect1[1] - vect2[1];
140 }
141 /*----------------------------------------------------------------------------*/
142 void vectAdd2D(const double v1[2], const double v2[2], double res[2])
143 {
144   res[0] = v1[0] + v2[0];
145   res[1] = v1[1] + v2[1];
146 }
147 /*----------------------------------------------------------------------------*/
148 void scalarMult2D(const double v[2], const double scalar, double res[2])
149 {
150   res[0] = scalar * v[0];
151   res[1] = scalar * v[1];
152 }
153 /*----------------------------------------------------------------------------*/
154 void normalize2d( double vect[2] )
155 {
156   double norm = NORM_2D(vect) ;
157   vect[0] /= norm ;
158   vect[1] /= norm ;
159 }
160 /*----------------------------------------------------------------------------*/
161 void iNormalize2d( int vect[2] )
162 {
163   double norm = NORM_2D(vect) ;
164   vect[0] = round( vect[0] / norm ) ;
165   vect[1] = round( vect[1] / norm ) ;
166 }
167 /*----------------------------------------------------------------------------*/
168 BOOL isPointInTriangle(const double point[2], const double a[2],
169                        const double b[2], const double c[2])
170 {
171   return (   areOnSameSideOfLine(point, a, b, c)
172           && areOnSameSideOfLine(point, b, a, c)
173           && areOnSameSideOfLine(point, c, a, b));
174 }
175 /*----------------------------------------------------------------------------*/
176 BOOL areOnSameSideOfLine(const double p1[2], const double p2[2],
177                          const double a[2], const double b[2])
178 {
179   // point are on the same if and only if (AB^AP1).(AB^AP2) >= 0
180   double ab[3];
181   double ap1[3];
182   double ap2[3];
183   double cp1[3];
184   double cp2[3];
185
186   ab[0] = b[0] - a[0];
187   ab[1] = b[1] - a[1];
188   ab[2] = 0.0;
189  
190   ap1[0] = p1[0] - a[0];
191   ap1[1] = p1[1] - a[1];
192   ap1[2] = 0.0;
193
194   ap2[0] = p2[0] - a[0];
195   ap2[1] = p2[1] - a[1];
196   ap2[2] = 0.0;
197
198   crossProduct(ab, ap1, cp1);
199   crossProduct(ab, ap2, cp2);
200
201   return (DOT_PROD_3D(cp1, cp2) >= 0.0);
202
203 }
204 /*----------------------------------------------------------------------------*/
205 void crossProduct( const double v1[3], const double v2[3], double res[3] )
206 {
207   /* save data to be able to use v1 o v2 as res */
208   double v10 = v1[0];
209   double v20 = v2[0];
210   double v11 = v1[1];
211   double v21 = v2[1];
212   res[0] = v11 * v2[2] - v1[2] * v21 ;
213   res[1] = v1[2] * v20 - v10 * v2[2] ;
214   res[2] = v10 * v21 - v11 * v20 ;
215 }
216 /*----------------------------------------------------------------------------*/
217 void vectSubstract3D(const double v1[3] ,const double v2[3], double res[3])
218 {
219   res[0] = v1[0] - v2[0];
220   res[1] = v1[1] - v2[1];
221   res[2] = v1[2] - v2[2];
222 }
223 /*----------------------------------------------------------------------------*/
224 void vectAdd3D(const double v1[3], const double v2[3], double res[3])
225 {
226   res[0] = v1[0] + v2[0];
227   res[1] = v1[1] + v2[1];
228   res[2] = v1[2] + v2[2];
229 }
230 /*----------------------------------------------------------------------------*/
231 void scalarMult3D(const double v[3], double scalar, double res[3])
232 {
233   res[0] = scalar * v[0];
234   res[1] = scalar * v[1];
235   res[2] = scalar * v[2];
236 }
237 /*----------------------------------------------------------------------------*/
238 void normalize3D( double vect[3] )
239 {
240   double norm = NORM_3D(vect) ;
241   vect[0] /= norm ;
242   vect[1] /= norm ;
243   vect[2] /= norm ;
244 }
245 /*----------------------------------------------------------------------------*/
246 void setToIdentity(double mat4D[4][4])
247 {
248   int i;
249   int j;
250   for (i = 0; i < 4; i++)
251   {
252     for (j = 0; j < 4; j++)
253     {
254       mat4D[i][j] = 0.0;
255     }
256     mat4D[i][i] = 1.0;
257   }
258 }
259 /*----------------------------------------------------------------------------*/
260 void mat4DMult(const double mat4D[4][4], const double vect3D[3], double res[3])
261 {
262   int i;
263   double res4D[4];
264   // w coordinate of the vector is supposed to be 1;
265   for (i = 0; i < 4; i++) {
266     res4D[i]  =  vect3D[0] * mat4D[i][0] + vect3D[1] * mat4D[i][1]
267                + vect3D[2] * mat4D[i][2] + mat4D[i][3];
268   }
269   res[0] = res4D[0] / res4D[3];
270   res[1] = res4D[1] / res4D[3];
271   res[2] = res4D[2] / res4D[3];
272 }
273 /*----------------------------------------------------------------------------*/
274 /* check if two values can be considered equal given an accurracy */
275 int safeEqual( double val1, double val2, double accuracy )
276 {
277   /* the val1 == val2 is put to avoid division by 0 */
278   return ( val1 == val2 ) || ( Abs( val1 - val2 ) < accuracy * Max( Abs(val1), Abs(val2 ) ) ) ;
279 }
280 /*----------------------------------------------------------------------------*/