1
2 /*
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) 1998-2001 - ENPC - Jean-Philippe Chancelier
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 #include "Scierror.h"
23 #include "sciprint.h"
24 #include "sci_malloc.h"
25 #include "Format.h"
26 #include "Contour.h"
27
28 #include "localization.h"
29
30 typedef void (level_f)(int ival, double Cont, double xncont, double yncont);
31 typedef void (*ptr_level_f)(int ival, double Cont, double xncont, double yncont);
32
33
34 static int contourI (ptr_level_f, double *, double *, double *, double *, int *, int *, int *);
35
36 static void
37 look(ptr_level_f, int i, int j, int ib, int jb, int qq, double Cont, int style);
38
39 static int ffnd (ptr_level_f, int, int, int, int, int,
40 int, int, int, int, int,
41 double, int *);
42
43 static int Gcont_size = 0;
44
45 static level_f GContStore2;
46 static void GContStore2Last(void);
47 static double x_cont(int i);
48 static double y_cont (int i);
49 static double phi_cont (int, int);
50 static double f_intercept (double, double, double, double, double);
51 static int not_same_sign (double val1, double val2);
52 static int get_itg_cont (int i, int j);
53 static void inc_itg_cont (int i, int j, int val);
54 static int oddp (int i);
55
56 /*-----------------------------------------------------------------------
57 * Level curves
58 * The computer journal vol 15 nul 4 p 382 (1972)
59 * from the Lisp Macsyma source (M.I.T)
60 * -------------------------------------------------------------------------*/
61
62 /*---------------------------------------------------------------------------
63 * General functions (could be changed in #define or
64 * inline functions to increase speed)
65 *---------------------------------------------------------------------------*/
66
67 static double *GX, *GY, *GZ;
68 static int Gn1, Gn2;
69
70 static double* Gxcont = NULL;
71 static double* Gycont = NULL;
72 static int currentContSize = 0;
73
InitValues(double * x,double * y,double * z,int n1,int n2)74 static void InitValues(double *x, double *y, double *z, int n1, int n2)
75 {
76 Gn1 = n1;
77 Gn2 = n2;
78 GX = x;
79 GY = y;
80 GZ = z;
81 }
82
83 /*--------return the value of f for a point on the grid-----*/
84
phi_cont(int i,int j)85 static double phi_cont(int i, int j)
86 {
87 return (GZ[i + Gn1 * j]);
88 }
89
90 /*---------return the coordinates between [xi,xj] along one axis
91 * for which the value of f is zCont */
92
f_intercept(double zCont,double fi,double xi,double fj,double xj)93 static double f_intercept(double zCont, double fi, double xi, double fj, double xj)
94 {
95 return (xi + (zCont - fi) * (xj - xi) / (fj - fi));
96 }
97
98 /* check for boundary points */
99
bdyp(int i,int j)100 static int bdyp(int i, int j)
101 {
102 return ( j == 0 || i == 0 || j == Gn2 - 1 || i == Gn1 - 1);
103 }
104
105 /* store or get flag values */
106
107 static int *itg_cont, *xbd_cont, *ybd_cont;
108
get_itg_cont(int i,int j)109 static int get_itg_cont(int i, int j)
110 {
111 return (itg_cont[i + Gn1 * j]);
112 }
113
inc_itg_cont(int i,int j,int val)114 static void inc_itg_cont(int i, int j, int val)
115 {
116 itg_cont[i + Gn1 * j] += val;
117 }
118
not_same_sign(double val1,double val2)119 static int not_same_sign(double val1, double val2)
120 {
121 if (ISNAN(val1) == 1 || ISNAN(val2) == 1)
122 {
123 return (0);
124 }
125 /** 0.0 est consid\'er\'e comme positif **/
126 if (val1 >= 0.0)
127 {
128 if (val2 < 0.0)
129 {
130 return (1) ;
131 }
132 else
133 {
134 return (0);
135 }
136 }
137 else
138 {
139 if (val2 >= 0.0)
140 {
141 return (1) ;
142 }
143 else
144 {
145 return (0);
146 }
147 }
148 }
149
oddp(int i)150 static int oddp(int i)
151 {
152 return (i == 1 || i == 3);
153 }
154
155 /*---------return the x-value of a grid point--------*/
156
x_cont(int i)157 static double x_cont(int i)
158 {
159 return GX[i] ;
160 }
161
162 /*---------return the y-value of a grid point --------*/
163
y_cont(int i)164 static double y_cont(int i)
165 {
166 return GY[i] ;
167 }
168
169
170 static char ContNumFormat[100];
171
172 /*--------------------------------------------------------------------
173 * the level curve is crossing the segment (i,j) (ib,jb)
174 * look store the level curve point and try to find the next segment to look at
175 * Cont: value of f along the contour
176 * ncont: number of contour
177 * c: indice of the contour Cont
178 *---------------------------------------------------------------------*/
179
look(ptr_level_f func,int i,int j,int ib,int jb,int qq,double Cont,int style)180 static void look(ptr_level_f func, int i, int j, int ib, int jb, int qq, double Cont, int style)
181 {
182 int ip = 0, jp = 0, im = 0, jm = 0, zds = 0, ent = 0, flag = 0, wflag = 0;
183 jp = j + 1;
184 ip = i + 1;
185 jm = j - 1;
186 im = i - 1;
187 /* on regarde comment est le segment de depart */
188 if (jb == jm)
189 {
190 flag = 1;
191 }
192 else
193 {
194 if (ib == im)
195 {
196 flag = 2;
197 }
198 else
199 {
200 if (jb == jp)
201 {
202 flag = 3;
203 }
204 else if (ib == ip)
205 {
206 flag = 4;
207 }
208 }
209 }
210 switch ( flag)
211 {
212 case 1 :
213 if (get_itg_cont(i, jm) > 1)
214 {
215 return;
216 }
217 ent = 1 ; /* le segment est vertical vers le bas */
218 /* Storing intersection point */
219 (*func)(0, Cont, x_cont(i),
220 f_intercept(Cont, phi_cont(i, jm),
221 y_cont(jm), phi_cont(i, j), y_cont(j)));
222 break;
223 case 2 :
224 if (get_itg_cont(im, j) == 1 || get_itg_cont(im, j) == 3)
225 {
226 return;
227 }
228 ent = 2 ; /* le segment est horizontal gauche */
229 /* Storing intersection point */
230 (*func)(0, Cont,
231 f_intercept(Cont, phi_cont(im, j),
232 x_cont(im), phi_cont(i, j), x_cont(i)), y_cont(j));
233 break ;
234 case 3 :
235 if (get_itg_cont(i, j) > 1)
236 {
237 return;
238 }
239 ent = 3 ; /* le segment est vertical haut */
240 /* Storing intersection point */
241 (*func)(0, Cont, x_cont(i), f_intercept(Cont, phi_cont(i, j),
242 y_cont(j), phi_cont(i, jp), y_cont(jp)));
243 break;
244 case 4 :
245 if (get_itg_cont(i, j) == 1 || get_itg_cont(i, j) == 3)
246 {
247 return;
248 }
249 ent = 4 ; /* le segment est horizontal droit */
250 /* Storing intersection point */
251 (*func)(0, Cont, f_intercept(Cont, phi_cont(i, j),
252 x_cont(i), phi_cont(ip, j), x_cont(ip)),
253 y_cont(j));
254 break;
255 default :
256 break;
257 }
258 wflag = 1;
259 while (wflag)
260 {
261 jp = j + 1;
262 ip = i + 1;
263 jm = j - 1;
264 im = i - 1;
265 switch (ent)
266 {
267 case 1 :
268 inc_itg_cont(i, jm, 2L);
269 ent = ffnd(func, i, ip, ip, i, j, j, jm, jm, ent, qq, Cont, &zds);
270 /* on calcule le nouveau point, ent donne la
271 direction du segment a explorer */
272 switch (ent)
273 {
274 case -1:
275 wflag = 0;
276 break;
277 case 1 :
278 i = ip ;
279 break;
280 case 2 :
281 i = ip;
282 j = jm;
283 break;
284 }
285 break;
286 case 2 :
287 inc_itg_cont(im, j, 1L);
288 ent = ffnd(func, i, i, im, im, j, jm, jm, j, ent, qq, Cont, &zds);
289 switch (ent)
290 {
291 case -1:
292 wflag = 0;
293 break;
294 case 2 :
295 j = jm ;
296 break;
297 case 3 :
298 i = im;
299 j = jm;
300 break;
301 }
302 break;
303 case 3 :
304 inc_itg_cont(i, j, 2L);
305 ent = ffnd(func, i, im, im, i, j, j, jp, jp, ent, qq, Cont, &zds);
306 switch (ent)
307 {
308 case -1:
309 wflag = 0;
310 break;
311 case 3 :
312 i = im;
313 break;
314 case 4 :
315 i = im;
316 j = jp;
317 break;
318 }
319 break;
320 case 4 :
321 inc_itg_cont(i, j, 1L);
322 ent = ffnd(func, i, i, ip, ip, j, jp, jp, j, ent, qq, Cont, &zds);
323 switch (ent)
324 {
325 case -1:
326 wflag = 0;
327 break;
328 case 4 :
329 j = jp;
330 break;
331 case 1 :
332 i = ip;
333 j = jp;
334 break;
335 }
336 break;
337 }
338
339 /** new segment is on the boundary **/
340 if (zds == 1)
341 {
342 switch (ent)
343 {
344 case 1 :
345 inc_itg_cont(i, (j - 1), 2L);
346 break ;
347 case 2 :
348 inc_itg_cont(i - 1, j, 1L);
349 break ;
350 case 3 :
351 inc_itg_cont(i, j, 2L);
352 break ;
353 case 4 :
354 inc_itg_cont(i, j, 1L);
355 break ;
356 }
357 /** we must quit the while loop **/
358 wflag = 0;
359 }
360 /** init point was inside the domain */
361 if (qq == 2)
362 {
363 switch (ent)
364 {
365 case 1 :
366 if (get_itg_cont (i, j - 1) > 1)
367 {
368 wflag = 0 ;
369 }
370 break ;
371 case 2 :
372 if (oddp(get_itg_cont(i - 1, j)))
373 {
374 wflag = 0 ;
375 }
376 break ;
377 case 3 :
378 if (get_itg_cont(i, j) > 1)
379 {
380 wflag = 0 ;
381 }
382 break ;
383 case 4 :
384 if (oddp(get_itg_cont(i, j)))
385 {
386 wflag = 0 ;
387 }
388 break ;
389 }
390 }
391 }
392 if (func == GContStore2)
393 {
394 GContStore2Last();
395 }
396 else
397 {
398 /* contour2di only computes level curves, not display them. */
399 sciprint(_("%s is only made to compute level curves and not display them.\n"), "Contourdi");
400 }
401 }
402
403
404 /*-------------------------------------------------------
405 * The function f is given on a grid and we want the level curves
406 * for the zCont[N[2]] values
407 * x : of size N[0] gives the x-values of the grid
408 * y : of size N[1] gives the y-values of the grid
409 * z : of size N[0]*N[1] gives the f-values on the grid
410 * style: size ncont (=N[2]) or empty int pointer
411 * gives the dash style for contour i
412 *-------------------------------------------------------*/
413
contourI(ptr_level_f func,double * x,double * y,double * z,double * zCont,int * N,int * style,int * err)414 static int contourI(ptr_level_f func, double *x, double *y, double *z, double *zCont, int *N, int *style, int *err)
415 {
416 int check = 1;
417 char *F = NULL;
418 int n1 = 0, n2 = 0, ncont = 0, i = 0, c = 0, j = 0, k = 0, n5 = 0;
419 int stylec = 0;
420 n1 = N[0];
421 n2 = N[1];
422 ncont = N[2];
423 F = getFPF();
424 if (F[0] == '\0')
425 {
426 ChoixFormatE1(ContNumFormat, zCont, N[2]);
427 }
428 InitValues(x, y, z, n1, n2);
429 n5 = 2 * (n1) + 2 * (n2) - 3;
430 /* Allocation */
431 Gcont_size = 0; /** initialize the array indices for storing contours **/
432 xbd_cont = (int*)MALLOC(n5 * sizeof(int));
433 ybd_cont = (int*)MALLOC(n5 * sizeof(int));
434 itg_cont = (int*)MALLOC(n1 * n2 * sizeof(int));
435 if ((xbd_cont == NULL) && n5 != 0)
436 {
437 check = 0;
438 }
439 if ((ybd_cont == NULL) && n5 != 0)
440 {
441 check = 0;
442 }
443 if ((itg_cont == NULL) && n1 * n2 != 0)
444 {
445 check = 0;
446 }
447 if (check == 0)
448 {
449 FREE(xbd_cont);
450 FREE(ybd_cont);
451 FREE(itg_cont);
452 Scierror(999, _("%s: No more memory.\n"), "contourI");
453 return -1;
454 }
455 /* just a parametrization of the boundary points */
456 for (i = 0 ; i < n2 ; i++)
457 {
458 ybd_cont[i] = i;
459 xbd_cont[i] = 0;
460 }
461 for (i = 1 ; i < n1 ; i++)
462 {
463 ybd_cont[n2 + i - 1] = n2 - 1;
464 xbd_cont[n2 + i - 1] = i ;
465 }
466 for (i = n2 - 2; i >= 0 ; i--)
467 {
468 ybd_cont[2 * n2 + n1 - 3 - i] = i;
469 xbd_cont[2 * n2 + n1 - 3 - i] = n1 - 1 ;
470 }
471 for (i = n1 - 2 ; i >= 0 ; i--)
472 {
473 ybd_cont[2 * n2 + 2 * n1 - 4 - i] = 0;
474 xbd_cont[2 * n2 + 2 * n1 - 4 - i] = i ;
475 }
476 for (c = 0 ; c < ncont ; c++)
477 {
478 stylec = (style != (int *) 0) ? style[c] : c;
479 /** itg-cont is a flag array to memorize checked parts of the grid **/
480 for (i = 0 ; i < n1; i++)
481 for (j = 0 ; j < n2 ; j++)
482 {
483 itg_cont[i + n1 * j] = 0;
484 }
485 /** all the boundary segments **/
486 for (k = 1 ; k < n5 ; k++)
487 {
488 int ib, jb;
489 i = xbd_cont[k] ;
490 j = ybd_cont[k];
491 ib = xbd_cont[k - 1] ;
492 jb = ybd_cont[k - 1];
493 if (not_same_sign (phi_cont(i, j) - zCont[c] ,
494 phi_cont(ib, jb) - zCont[c]))
495 {
496 look(func, i, j, ib, jb, 1L, zCont[c], stylec);
497 }
498 }
499 /** inside segments **/
500 for (i = 1 ; i < n1 - 1; i++)
501 for (j = 1 ; j < n2 - 1 ; j++)
502 if (not_same_sign (phi_cont(i, j) - zCont[c] ,
503 phi_cont(i, j - 1) - zCont[c]))
504 {
505 look(func, i, j, i, j - 1, 2L, zCont[c], stylec);
506 }
507 }
508 FREE(xbd_cont);
509 FREE(ybd_cont);
510 FREE(itg_cont);
511
512 return 0;
513 }
514
C2F(contourif)515 int C2F(contourif)(double *x, double *y, double *z, int *n1, int *n2, int *flagnz, int *nz, double *zz, int *style)
516 {
517 int err = 0;
518 static double *zconst = NULL;
519 double zmin = 0., zmax = 0.;
520 int N[3], i = 0;
521
522 zmin = (double) Mini(z, *n1 * (*n2));
523 zmax = (double) Maxi(z, *n1 * (*n2));
524
525 if (*flagnz == 0)
526 {
527 if ((zconst = (double*)MALLOC((*nz) * sizeof(double))) == 0)
528 {
529 Scierror(999, _("%s: No more memory.\n"), "contourif");
530 return -1;
531 }
532 for (i = 0 ; i < *nz ; i++)
533 {
534 zconst[i] = zmin + (i + 1) * (zmax - zmin) / (*nz + 1);
535 }
536 N[0] = *n1;
537 N[1] = *n2;
538 N[2] = *nz;
539 contourI(GContStore2, x, y, z, zconst, N, style, &err);
540 FREE(zconst);
541 zconst = NULL;
542 }
543 else
544 {
545 N[0] = *n1;
546 N[1] = *n2;
547 N[2] = *nz;
548 contourI(GContStore2, x, y, z, zz, N, style, &err);
549 }
550
551 return 0;
552 }
553
554
555
556
557
558 /*-----------------------------------------------------------------------
559 * ffnd : cette fonction recoit en entree quatre points
560 * on sait que la courbe de niveau passe entre le point 1 et le quatre
561 * on cherche a savoir ou elle resort,
562 * et on fixe une nouvelle valeur de ent qui indiquera le segment
563 * suivant a explorer
564 *-----------------------------------------------------------------------*/
565
ffnd(ptr_level_f func,int i1,int i2,int i3,int i4,int jj1,int jj2,int jj3,int jj4,int ent,int qq,double Cont,int * zds)566 static int ffnd (ptr_level_f func, int i1, int i2, int i3, int i4, int jj1, int jj2, int jj3, int jj4, int ent, int qq, double Cont, int *zds)
567 {
568 double phi1 = 0., phi2 = 0., phi3 = 0., phi4 = 0., xav = 0., yav = 0., phiav = 0.;
569 int revflag = 0, i = 0;
570 phi1 = phi_cont(i1, jj1) - Cont;
571 phi2 = phi_cont(i2, jj2) - Cont;
572 phi3 = phi_cont(i3, jj3) - Cont;
573 phi4 = phi_cont(i4, jj4) - Cont;
574 revflag = 0;
575 *zds = 0;
576 /* le point au centre du rectangle */
577 xav = (x_cont(i1) + x_cont(i3)) / 2.0 ;
578 yav = (y_cont(jj1) + y_cont(jj3)) / 2.0 ;
579 phiav = (phi1 + phi2 + phi3 + phi4) / 4.0;
580 if (ISNAN(phiav) == 1)
581 {
582 return -1;
583 }
584 if ( not_same_sign(phiav, phi4))
585 {
586 int l1, k1;
587 double phi;
588 revflag = 1 ;
589 l1 = i4;
590 k1 = jj4;
591 i4 = i1;
592 jj4 = jj1;
593 i1 = l1;
594 jj1 = k1;
595 l1 = i3;
596 k1 = jj3;
597 i3 = i2;
598 jj3 = jj2;
599 i2 = l1;
600 jj2 = k1;
601 phi = phi1;
602 phi1 = phi4;
603 phi4 = phi;
604 phi = phi3;
605 phi3 = phi2;
606 phi2 = phi;
607 }
608 /* on stocke un nouveau point */
609 (*func)(1, Cont, f_intercept(0.0, phi1, x_cont(i1), phiav, xav),
610 f_intercept(0.0, phi1, y_cont(jj1), phiav, yav));
611 /*
612 * on parcourt les segments du rectangle pour voir sur quelle face
613 * on sort
614 */
615 for (i = 0 ; ; i++)
616 {
617 int l1, k1;
618 double phi;
619 if (not_same_sign (phi1, phi2)) /** sortir du for **/
620 {
621 break ;
622 }
623 if (phiav != 0.0)
624 {
625 (*func)(1, Cont, f_intercept(0.0, phi2, x_cont(i2), phiav, xav),
626 f_intercept(0.0, phi2, y_cont(jj2), phiav, yav));
627 }
628 /** on permutte les points du rectangle **/
629 l1 = i1;
630 k1 = jj1;
631 i1 = i2;
632 jj1 = jj2;
633 i2 = i3;
634 jj2 = jj3;
635 i3 = i4;
636 jj3 = jj4;
637 i4 = l1;
638 jj4 = k1;
639 phi = phi1;
640 phi1 = phi2;
641 phi2 = phi3;
642 phi3 = phi4;
643 phi4 = phi;
644 }
645 (*func)(1, Cont, f_intercept(0.0, phi1, x_cont(i1), phi2, x_cont(i2)),
646 f_intercept(0.0, phi1, y_cont(jj1), phi2, y_cont(jj2)));
647 if (qq == 1 && bdyp(i1, jj1) && bdyp(i2, jj2))
648 {
649 *zds = 1;
650 }
651 if (revflag == 1 && ! oddp (i))
652 {
653 i = i + 2;
654 }
655 return (1 + ( (i + ent + 2) % 4));
656 }
657
658 /*--------------------------------------------------------------
659 * Following code is used to store the current level curves as
660 * double in order to access to the stored data at Scilab level
661 *----------------------------------------------------------------*/
662
663 static int last = -1;
664 static int count = 0;
665
666 /** used to bring back data to Scilab Stack **/
667
C2F(getconts)668 int C2F(getconts)(double **x, double **y, int *mm, int *n)
669 {
670 *x = Gxcont;
671 *y = Gycont;
672 *mm = 1;
673 *n = Gcont_size;
674 return 0;
675 }
676
GContStore2(int ival,double Cont,double xncont,double yncont)677 static void GContStore2(int ival, double Cont, double xncont, double yncont)
678 {
679 int n = 0;
680 if (ival == 0)
681 {
682 /* Here : ival == 0 means stop the current level curve and
683 * store data at the end but do reset Gcont_size to zero
684 */
685 n = Gcont_size + 2;
686 if (Gxcont == NULL)
687 {
688 Gxcont = (double*)MALLOC(n * sizeof(double));
689 Gycont = (double*)MALLOC(n * sizeof(double));
690 }
691 else
692 {
693 if (n > currentContSize)
694 {
695 currentContSize = (int)(n * 1.1);
696 Gxcont = (double*)REALLOC(Gxcont, currentContSize * sizeof(double));
697 Gycont = (double*)REALLOC(Gycont, currentContSize * sizeof(double));
698 }
699 }
700 if ((Gxcont == NULL) && n != 0)
701 {
702 return ;
703 }
704 if ((Gycont == NULL) && n != 0)
705 {
706 return;
707 }
708 Gxcont[Gcont_size] = Cont;
709 if (last != -1 && last < n)
710 {
711 Gycont[last] = count;
712 }
713 last = Gcont_size;
714 Gcont_size++;
715 count = 0;
716 }
717 else
718 {
719 n = Gcont_size + 1;
720 if (Gxcont == NULL)
721 {
722 Gxcont = (double*)MALLOC(n * sizeof(double));
723 Gycont = (double*)MALLOC(n * sizeof(double));
724 }
725 else
726 {
727 if (n > currentContSize)
728 {
729 currentContSize = (int)(n * 1.1);
730 Gxcont = (double*)REALLOC(Gxcont, currentContSize * sizeof(double));
731 Gycont = (double*)REALLOC(Gycont, currentContSize * sizeof(double));
732 }
733 }
734 if ((Gxcont == NULL) && n != 0)
735 {
736 return ;
737 }
738 if ((Gycont == NULL) && n != 0)
739 {
740 return;
741 }
742 }
743
744 if (Gxcont)
745 {
746 Gxcont[Gcont_size] = xncont;
747 }
748 if (Gycont)
749 {
750 Gycont[Gcont_size++] = yncont;
751 }
752 count++;
753 }
754
GContStore2Last(void)755 static void GContStore2Last(void)
756 {
757 if (last != -1)
758 {
759 Gycont[last] = count;
760 }
761 }
762
763
764
765
766
767
768
769
770
771
772
773