* Bug 16365 fixed: median(m,'r'|'c') was wrong after 5dc990
[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  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14  *
15  */
16
17 /*------------------------------------------------------------------------
18  *    Graphic library
19  --------------------------------------------------------------------------*/
20
21 #include "math_graphics.h"
22
23 /*
24  * we use spConfig.h for machine constants
25  * XXX : spConfig should be merged and unified
26  *       with other machine constant scilab code
27  */
28
29 #define spINSIDE_SPARSE
30 #include "../../sparse/includes/spConfig.h"
31
32 double Mini(const double vect[], int n)
33 {
34     int i = 0;
35     double vmin = LARGEST_REAL;
36     for (i = 0 ; i < n ; i++)
37     {
38         /*    if (isinf(vect[i])== 0 && isnan(vect[i])==0 && vect[i] < vmin)  */
39         if (finite(vect[i]) == 1 && vect[i] < vmin)
40         {
41             vmin = vect[i];
42         }
43     }
44
45     return vmin;
46 }
47
48 double Maxi(const double vect[], int n)
49 {
50     int i = 0;
51     double maxi = - LARGEST_REAL;
52     for (i = 0 ; i < n ; i++)
53     {
54         /* if (isinf(vect[i])== 0 && isnan(vect[i])==0 && vect[i] > maxi) */
55         if (finite(vect[i]) == 1 && vect[i] > maxi)
56         {
57             maxi = vect[i];
58         }
59     }
60
61     return maxi;
62 }
63
64 void MiniMaxi(const double vect[], int n, double * const min, double * const max)
65 {
66     int i = 0;
67     double _min = LARGEST_REAL;
68     double _max = -LARGEST_REAL;
69     for (; i < n ; i++)
70     {
71         /*    if ( isinf(vect[i])== 0 && isnan(vect[i])==0 && vect[i] < vmin)  */
72         if (finite(vect[i]) == 1)
73         {
74             if (vect[i] < _min)
75             {
76                 _min = vect[i];
77             }
78             if (vect[i] > _max)
79             {
80                 _max = vect[i];
81             }
82         }
83     }
84
85     *min = _min;
86     *max = _max;
87 }
88
89 /*----------------------------------------------------------------------------*/
90
91 /* perform the rotation of point from to point dest  */
92 void rotate2D(double from[2], double center[2], double angle, double dest[2])
93 {
94     rotate2Dim(from, center, cos(angle), sin(angle), dest);
95 }
96
97 /*----------------------------------------------------------------------------*/
98 /* perform the rotation of point from to point to. */
99 /* the angle is directly given with its sine and cosine for speed */
100 void rotate2Dim(double from[2]   ,
101                 double center[2] ,
102                 double cosAngle  ,
103                 double sinAngle  ,
104                 double dest[2])
105 {
106     double diff[2];
107
108     /* put the center to (0,0) */
109     diff[0] = from[0] - center[0];
110     diff[1] = from[1] - center[1];
111
112     /* turn and translate back */
113     dest[0] = diff[0] * cosAngle - diff[1] * sinAngle + center[0];
114     dest[1] = diff[0] * sinAngle + diff[1] * cosAngle + center[1];
115 }
116 /*----------------------------------------------------------------------------*/
117 /* perform the translation of point from to point to with vector trans */
118 void translate2D(double from[2], double trans[2], double dest[2])
119 {
120     dest[0] = from[0] + trans[0];
121     dest[1] = from[1] + trans[1];
122 }
123 /*----------------------------------------------------------------------------*/
124 void iTranslate2D(int from[2], int trans[2], int dest[2])
125 {
126     dest[0] = from[0] + trans[0];
127     dest[1] = from[1] + trans[1];
128 }
129 /*----------------------------------------------------------------------------*/
130 void vectSubstract2D(const double vect1[2], const double vect2[], double res[2])
131 {
132     res[0] = vect1[0] - vect2[0];
133     res[1] = vect1[1] - vect2[1];
134 }
135 /*----------------------------------------------------------------------------*/
136 void vectAdd2D(const double v1[2], const double v2[2], double res[2])
137 {
138     res[0] = v1[0] + v2[0];
139     res[1] = v1[1] + v2[1];
140 }
141 /*----------------------------------------------------------------------------*/
142 void scalarMult2D(const double v[2], const double scalar, double res[2])
143 {
144     res[0] = scalar * v[0];
145     res[1] = scalar * v[1];
146 }
147 /*----------------------------------------------------------------------------*/
148 void normalize2d(double vect[2])
149 {
150     double norm = NORM_2D(vect);
151     vect[0] /= norm;
152     vect[1] /= norm;
153 }
154 /*----------------------------------------------------------------------------*/
155 void iNormalize2d(int vect[2])
156 {
157     double norm = NORM_2D(vect);
158     vect[0] = scilab_round(vect[0] / norm);
159     vect[1] = scilab_round(vect[1] / norm);
160 }
161 /*----------------------------------------------------------------------------*/
162 BOOL isPointInTriangle(const double point[2], const double a[2],
163                        const double b[2], const double c[2])
164 {
165     return (   areOnSameSideOfLine(point, a, b, c)
166                && areOnSameSideOfLine(point, b, a, c)
167                && areOnSameSideOfLine(point, c, a, b));
168 }
169 /*----------------------------------------------------------------------------*/
170 BOOL areOnSameSideOfLine(const double p1[2], const double p2[2],
171                          const double a[2], const double b[2])
172 {
173     // point are on the same if and only if (AB^AP1).(AB^AP2) >= 0
174     double ab[3];
175     double ap1[3];
176     double ap2[3];
177     double cp1[3];
178     double cp2[3];
179
180     ab[0] = b[0] - a[0];
181     ab[1] = b[1] - a[1];
182     ab[2] = 0.0;
183
184     ap1[0] = p1[0] - a[0];
185     ap1[1] = p1[1] - a[1];
186     ap1[2] = 0.0;
187
188     ap2[0] = p2[0] - a[0];
189     ap2[1] = p2[1] - a[1];
190     ap2[2] = 0.0;
191
192     crossProduct(ab, ap1, cp1);
193     crossProduct(ab, ap2, cp2);
194
195     return (DOT_PROD_3D(cp1, cp2) >= 0.0);
196
197 }
198 /*----------------------------------------------------------------------------*/
199 void crossProduct(const double v1[3], const double v2[3], double res[3])
200 {
201     /* save data to be able to use v1 o v2 as res */
202     double v10 = v1[0];
203     double v20 = v2[0];
204     double v11 = v1[1];
205     double v21 = v2[1];
206     res[0] = v11 * v2[2] - v1[2] * v21 ;
207     res[1] = v1[2] * v20 - v10 * v2[2] ;
208     res[2] = v10 * v21 - v11 * v20 ;
209 }
210 /*----------------------------------------------------------------------------*/
211 void vectSubstract3D(const double v1[3] , const double v2[3], double res[3])
212 {
213     res[0] = v1[0] - v2[0];
214     res[1] = v1[1] - v2[1];
215     res[2] = v1[2] - v2[2];
216 }
217 /*----------------------------------------------------------------------------*/
218 void vectAdd3D(const double v1[3], const double v2[3], double res[3])
219 {
220     res[0] = v1[0] + v2[0];
221     res[1] = v1[1] + v2[1];
222     res[2] = v1[2] + v2[2];
223 }
224 /*----------------------------------------------------------------------------*/
225 void scalarMult3D(const double v[3], double scalar, double res[3])
226 {
227     res[0] = scalar * v[0];
228     res[1] = scalar * v[1];
229     res[2] = scalar * v[2];
230 }
231 /*----------------------------------------------------------------------------*/
232 void normalize3D(double vect[3])
233 {
234     double norm = NORM_3D(vect) ;
235     vect[0] /= norm ;
236     vect[1] /= norm ;
237     vect[2] /= norm ;
238 }
239 /*----------------------------------------------------------------------------*/
240 void setToIdentity(double mat4D[4][4])
241 {
242     int i = 0;
243     int j = 0;
244     for (i = 0; i < 4; i++)
245     {
246         for (j = 0; j < 4; j++)
247         {
248             mat4D[i][j] = 0.0;
249         }
250         mat4D[i][i] = 1.0;
251     }
252 }
253 /*----------------------------------------------------------------------------*/
254 void mat4DMult(const double mat4D[4][4], const double vect3D[3], double res[3])
255 {
256     int i = 0;
257     double res4D[4];
258     // w coordinate of the vector is supposed to be 1;
259     for (i = 0; i < 4; i++)
260     {
261         res4D[i]  =  vect3D[0] * mat4D[i][0] + vect3D[1] * mat4D[i][1]
262                      + vect3D[2] * mat4D[i][2] + mat4D[i][3];
263     }
264     res[0] = res4D[0] / res4D[3];
265     res[1] = res4D[1] / res4D[3];
266     res[2] = res4D[2] / res4D[3];
267 }
268 /*----------------------------------------------------------------------------*/