Finish to replace sciprint by Scierror when needed in graphics
[scilab.git] / scilab / modules / graphics / src / c / Contour.c
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  * 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 <string.h> /* in case of dbmalloc use */
19 #include <stdio.h>
20 #include <math.h>
21 #include "math_graphics.h"
22 #include "PloEch.h"
23 #include "Scierror.h"
24 #include "sciprint.h"
25 #include "MALLOC.h"
26 #include "Format.h"
27 #include "Contour.h"
28
29 #include "GetProperty.h"
30 #include "localization.h"
31
32 typedef void (level_f)(int ival, double Cont, double xncont,double yncont);
33 typedef void (*ptr_level_f)(int ival, double Cont, double xncont,double yncont);
34
35
36 static int contourI (ptr_level_f,double *, double *, double *, double *, int *, int *, int *);
37
38 static void
39 look(ptr_level_f, int i, int j, int ib,int jb, int qq,double Cont, int style);
40
41 static int ffnd (ptr_level_f,int,int,int,int,int,
42                              int,int,int,int,int,
43                              double,int *);
44
45 static int Gcont_size = 0;
46
47 static level_f GContStore2;
48 static void GContStore2Last(void);
49 static double x_cont(int i);
50 static double y_cont (int i);
51 static double phi_cont (int, int); 
52 static double f_intercept  (double, double, double, double, double );
53 static int not_same_sign  (double val1, double val2); 
54 static int get_itg_cont  (int i, int j); 
55 static void inc_itg_cont  (int i, int j, int val); 
56 static int oddp  (int i);
57
58 /*-----------------------------------------------------------------------
59  *  Level curves 
60  *  The computer journal vol 15 nul 4 p 382 (1972)
61  *  from the Lisp Macsyma source (M.I.T)
62  * -------------------------------------------------------------------------*/
63
64 /*---------------------------------------------------------------------------
65  * General functions (could be changed in #define or
66  *   inline functions to increase speed)
67  *---------------------------------------------------------------------------*/
68
69 static double *GX,*GY,*GZ;
70 static int Gn1,Gn2;
71
72 static double * Gxcont = NULL ;
73 static double * Gycont = NULL ;
74
75 static void InitValues(double *x, double *y, double *z, int n1, int n2)
76 {
77   Gn1=n1;  Gn2=n2;  GX = x;  GY = y;  GZ = z;
78 }
79
80 /*--------return the  value of f for a point on the grid-----*/
81
82 static double phi_cont(int i, int j)
83 {
84   return(GZ[i+Gn1*j]);
85 }
86
87 /*---------return the coordinates between  [xi,xj] along one axis 
88  *  for which the value of f is zCont */ 
89
90 static double f_intercept(double zCont, double fi, double xi, double fj, double xj)
91 {
92   return( xi+ (zCont-fi)*(xj-xi)/ (fj-fi));
93 }
94
95 /* check for boundary points */
96
97 static int  bdyp(int i, int j)
98 {
99   return (  j == 0 || i == 0 || j == Gn2-1 || i == Gn1-1 );
100 }
101
102 /* store or get flag values */
103
104 static  int *itg_cont, *xbd_cont,*ybd_cont;
105
106 static int get_itg_cont(int i, int j)
107 {
108   return( itg_cont[i+Gn1*j]);
109 }
110
111 static void inc_itg_cont(int i, int j, int val)
112 {
113   itg_cont[i+Gn1*j] += val;
114 }
115
116 static int not_same_sign(double val1, double val2)
117 {
118   if ( ISNAN(val1) ==1 || ISNAN(val2) == 1) return(0);
119   /** 0.0 est consid\'er\'e comme positif **/
120   if ( val1 >= 0.0) 
121     {
122       if (val2 < 0.0) return(1) ; else return(0);}
123   else 
124     {
125       if ( val2 >= 0.0) return(1) ; else return(0);}
126 }
127
128 static int oddp(int i) { return( i == 1 || i ==3 );}
129
130 /*---------return the x-value of a grid point--------*/
131
132 static double x_cont(int i) { return GX[i] ;}
133
134 /*---------return the y-value of a grid point --------*/
135
136 static double y_cont(int i) { return GY[i] ;}
137
138
139 static char   ContNumFormat[100];
140
141 /*--------------------------------------------------------------------
142 *  the level curve is crossing the segment (i,j) (ib,jb)
143 *  look store the level curve point and try to find the next segment to look at
144 *  Cont: value of f along the contour 
145 *  ncont: number of contour 
146 *  c: indice of the contour Cont 
147 *---------------------------------------------------------------------*/
148
149 static void look(ptr_level_f func, int i, int j, int ib, int jb, int qq, double Cont, int style)
150 {
151   int ip,jp,im,jm,zds,ent=0,flag=0,wflag;
152   jp= j+1; ip= i+1; jm=j-1;im=i-1;
153   /*  on regarde comment est le segment de depart */
154   if  ( jb == jm)  flag = 1; 
155   else  { 
156     if ( ib == im ) flag = 2 ;
157     else  {
158       if ( jb == jp ) flag = 3 ;
159       else  if ( ib == ip ) flag = 4;}}
160   switch  (  flag)
161   {
162   case  1 :
163     if  (get_itg_cont(i,jm) > 1) return;
164     ent=1 ; /* le segment est vertical vers le bas */
165     /* Storing intersection point */
166     (*func)(0,Cont, x_cont(i), 
167       f_intercept(Cont,phi_cont(i,jm),
168       y_cont(jm),phi_cont(i,j),y_cont(j)));
169     break;
170   case 2 : 
171     if  (get_itg_cont(im,j) == 1 || get_itg_cont(im,j)==3 ) return;
172     ent=2 ; /* le segment est horizontal gauche */
173     /* Storing intersection point */
174     (*func)( 0,Cont,
175       f_intercept(Cont,phi_cont(im,j),
176       x_cont(im),phi_cont(i,j),x_cont(i)), y_cont(j));
177     break ; 
178   case 3 :
179     if  (get_itg_cont(i,j) > 1 ) return;
180     ent=3 ; /* le segment est vertical haut */
181     /* Storing intersection point */
182     (*func)(0,Cont,x_cont(i), f_intercept(Cont,phi_cont(i,j),
183       y_cont(j),phi_cont(i,jp),y_cont(jp)));
184     break ;
185   case 4 :
186     if  (get_itg_cont(i,j) == 1 || get_itg_cont(i,j)==3 ) return;
187     ent=4 ; /* le segment est horizontal droit */
188     /* Storing intersection point */
189     (*func)(0,Cont,f_intercept(Cont,phi_cont(i,j),
190       x_cont(i),phi_cont(ip,j),x_cont(ip)),
191       y_cont(j));
192     break;
193   default :
194     break;
195   }
196   wflag=1;
197   while ( wflag) 
198   { 
199     jp= j+1; ip= i+1; jm=j-1;im=i-1;
200     switch  ( ent) 
201     {case 1 :
202     inc_itg_cont(i,jm,2L);
203     ent = ffnd(func,i,ip,ip,i,j,j,jm,jm,ent,qq,Cont,&zds);
204     /* on calcule le nouveau point, ent donne la 
205     direction du segment a explorer */
206     switch ( ent)
207     {
208     case -1: wflag=0; break;
209     case 1 : i=ip ; break ;
210     case 2 : i=ip;j=jm; break ;
211     }
212     break ;
213     case 2  :
214       inc_itg_cont(im,j,1L);
215       ent = ffnd(func,i,i,im,im,j,jm,jm,j,ent,qq,Cont,&zds);
216       switch ( ent)
217       { 
218       case -1: wflag=0; break;
219       case 2 : j = jm ;break ;
220       case  3  : i=im;j=jm; break ;
221       }
222       break ;
223     case 3 :
224       inc_itg_cont(i,j,2L);
225       ent = ffnd(func,i,im,im,i,j,j,jp,jp,ent,qq,Cont,&zds);
226       switch ( ent)
227       { 
228       case -1: wflag=0; break;
229       case 3 : i=im; break ;
230       case 4 : i=im;j=jp; break ;
231       }
232       break ;
233     case 4 :
234       inc_itg_cont(i,j,1L);
235       ent = ffnd(func,i,i,ip,ip,j,jp,jp,j,ent,qq,Cont,&zds);
236       switch ( ent)
237       {
238       case -1: wflag=0; break;
239       case 4 :j=jp;break ;
240       case 1 :i=ip;j=jp;break ;
241       }
242       break ;
243     }
244
245     /** new segment is on the boundary **/
246     if ( zds == 1) 
247     {
248       switch ( ent) 
249       {
250       case 1 : inc_itg_cont(i,(j-1),2L); break ; 
251       case 2 : inc_itg_cont(i-1,j,1L);  break ; 
252       case 3 : inc_itg_cont(i,j,2L); break ; 
253       case 4 : inc_itg_cont(i,j,1L); break ; 
254       }
255       /** we must quit the while loop **/
256       wflag = 0 ;
257     }
258     /**  init point was inside the domain */
259     if ( qq == 2) 
260     {
261       switch ( ent) 
262       {
263       case 1 : if  ( get_itg_cont (i,j-1)  > 1) wflag = 0 ; break ; 
264       case 2 : if  ( oddp(get_itg_cont(i-1,j))) wflag = 0 ; break ; 
265       case 3 : if  ( get_itg_cont(i,j) > 1)     wflag = 0 ; break ; 
266       case 4 : if  ( oddp(get_itg_cont(i,j)))   wflag = 0 ; break ; 
267       }
268     }
269   }
270   if ( func == GContStore2 )
271   {
272     GContStore2Last();
273   }
274   else
275   {
276     /* contour2di only computes level curves, not display them. */
277     sciprint(_("%s is only made to compute level curves and not display them.\n"),"Contourdi");
278   }
279 }
280
281
282 /*-------------------------------------------------------
283 *  The function f is given on a grid and we want the level curves 
284 *  for the zCont[N[2]] values 
285 *  x : of size N[0] gives the x-values of the grid 
286 *  y : of size N[1] gives the y-values of the grid 
287 *  z : of size N[0]*N[1]  gives the f-values on the grid 
288 *  style: size ncont (=N[2]) or empty int pointer 
289 *  gives the dash style for contour i
290 *-------------------------------------------------------*/
291
292 static int contourI(ptr_level_f func, double *x, double *y, double *z, double *zCont, int *N, int *style, int *err)
293 {
294   int check = 1;
295   char *F;
296   int n1,n2,ncont,i,c,j,k,n5;
297   int stylec;
298   n1=N[0];n2=N[1];ncont=N[2];
299   F=getFPF();
300   if ( F[0] == '\0') 
301     ChoixFormatE1(ContNumFormat,zCont,N[2]);
302   InitValues(x,y,z,n1,n2);
303   n5 =  2*(n1)+2*(n2)-3;
304   /* Allocation */
305   Gcont_size = 0; /** initialize the array indices for storing contours **/
306   xbd_cont = MALLOC( n5 * sizeof(int) ) ;
307   ybd_cont = MALLOC( n5 * sizeof(int) ) ;
308   itg_cont = MALLOC( n1*n2 * sizeof(int) ) ;
309   if ( (xbd_cont == NULL) && n5 != 0) check= 0;
310   if ( (ybd_cont == NULL) && n5 != 0) check= 0;
311   if ( (itg_cont == NULL) && n1*n2 != 0) check= 0;
312   if ( check == 0) 
313   {
314     FREE( xbd_cont ) ;
315     FREE( ybd_cont ) ;
316     FREE( itg_cont ) ;
317     Scierror(999, _("%s: No more memory.\n"),"contourI");
318     return -1;
319   }
320   /* just a parametrization of the boundary points */
321   for ( i = 0 ; i <  n2 ; i++)
322   {
323     ybd_cont[i] = i ;
324     xbd_cont[i] = 0 ;
325   }
326   for ( i = 1 ; i <  n1 ; i++)
327   {
328     ybd_cont[n2+i-1] = n2-1 ;
329     xbd_cont[n2+i-1] = i  ;
330   }
331   for ( i = n2-2;  i >= 0  ; i--)
332   {
333     ybd_cont[2*n2 +n1-3-i] = i ;
334     xbd_cont[2*n2 +n1-3-i] = n1-1  ;
335   }
336   for ( i = n1-2 ; i >= 0 ; i--)
337   {
338     ybd_cont[2*n2 +2*n1-4-i] = 0 ;
339     xbd_cont[2*n2 +2*n1-4-i] = i  ;
340   }
341   for ( c= 0 ; c < ncont ; c++)
342   {
343     stylec = ( style != (int *) 0) ? style[c] : c;
344     /** itg-cont is a flag array to memorize checked parts of the grid **/
345     for ( i = 0 ; i < n1; i++)
346       for ( j =0 ; j < n2 ; j++)
347         itg_cont[i+n1*j]=0 ;
348     /** all the boundary segments **/
349     for ( k = 1 ; k < n5 ; k++)
350     { int ib,jb;
351     i = xbd_cont[k] ; j = ybd_cont[k];
352     ib = xbd_cont[k-1] ; jb= ybd_cont[k-1];
353     if  (not_same_sign (phi_cont(i,j)-zCont[c] , 
354       phi_cont(ib,jb)-zCont[c]))
355       look(func,i,j,ib,jb,1L,zCont[c],stylec);
356     }
357     /** inside segments **/
358     for ( i = 1 ; i < n1-1; i++)
359       for ( j = 1 ; j < n2-1 ; j++)
360         if  (not_same_sign ( phi_cont(i,j)-zCont[c] , 
361           phi_cont(i, j-1)-zCont[c]))
362           look(func,i,j,i,j-1,2L,zCont[c],stylec);
363   }
364   FREE( xbd_cont ) ;
365   FREE( ybd_cont ) ;
366   FREE( itg_cont ) ;
367         
368         return 0;
369 }
370
371 int C2F(contourif)(double *x, double *y, double *z, int *n1, int *n2, int *flagnz, int *nz, double *zz, int *style)
372 {
373   int err=0;
374   static double *zconst;
375   double zmin,zmax;
376   int N[3],i;
377
378   zmin=(double) Mini(z,*n1*(*n2)); 
379   zmax=(double) Maxi(z,*n1*(*n2));
380
381   if (*flagnz==0)
382     {
383       if ( ( zconst = MALLOC( (*nz) * sizeof(double) ) ) == 0 ) 
384         {
385           Scierror(999, _("%s: No more memory.\n"),"contourif");
386           return -1;
387         }
388       for ( i =0 ; i < *nz ; i++) 
389         zconst[i]=zmin + (i+1)*(zmax-zmin)/(*nz+1);
390       N[0]= *n1;N[1]= *n2;N[2]= *nz;
391       contourI(GContStore2,x,y,z,zconst,N,style,&err);
392       FREE(zconst) ;
393       zconst = NULL ;
394     }
395   else
396     {
397       N[0]= *n1;N[1]= *n2;N[2]= *nz;
398       contourI(GContStore2,x,y,z,zz,N,style,&err);
399     }
400
401   return 0 ;
402 }
403
404
405
406
407
408 /*-----------------------------------------------------------------------
409  *   ffnd : cette fonction  recoit en entree quatre points 
410  *       on sait que la courbe de niveau passe entre le point 1 et le quatre 
411  *       on cherche a savoir ou elle resort, 
412  *       et on fixe une nouvelle valeur de ent qui indiquera le segment 
413  *       suivant a explorer 
414  *-----------------------------------------------------------------------*/
415
416 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)
417 {
418   double phi1,phi2,phi3,phi4,xav,yav,phiav;
419   int revflag,i;
420   phi1=phi_cont(i1,jj1)-Cont;
421   phi2=phi_cont(i2,jj2)-Cont;
422   phi3=phi_cont(i3,jj3)-Cont;
423   phi4=phi_cont(i4,jj4)-Cont;
424   revflag = 0;
425   *zds = 0;
426   /* le point au centre du rectangle */
427   xav = ( x_cont(i1)+ x_cont(i3))/2.0 ; 
428   yav = ( y_cont(jj1)+ y_cont(jj3))/2.0 ; 
429   phiav = ( phi1+phi2+phi3+phi4) / 4.0;
430   if (ISNAN(phiav)==1) 
431     {
432       return -1;
433     }
434   if (  not_same_sign( phiav,phi4)) 
435     {
436       int l1, k1; 
437       double phi;
438       revflag = 1 ; 
439       l1= i4; k1= jj4;
440       i4=i1; jj4 = jj1; i1= l1; jj1= k1;
441       l1= i3; k1= jj3;
442       i3=i2; jj3= jj2; i2=l1; jj2= k1;
443       phi = phi1; phi1 = phi4; phi4= phi;
444       phi = phi3; phi3 = phi2; phi2= phi;
445     }
446   /* on stocke un nouveau point  */
447   (*func)(1,Cont,f_intercept(0.0,phi1,x_cont(i1),phiav,xav),
448             f_intercept(0.0,phi1,y_cont(jj1),phiav,yav));
449   /*
450    * on parcourt les segments du rectangle pour voir sur quelle face
451    * on sort 
452    */
453   for  ( i = 0 ;  ; i++)
454     { int l1,k1;
455       double phi;
456       if ( not_same_sign ( phi1,phi2))   /** sortir du for **/ break ; 
457       if  ( phiav != 0.0 ) 
458         {
459           (*func)(1,Cont,f_intercept(0.0,phi2,x_cont(i2),phiav,xav),
460                     f_intercept(0.0,phi2,y_cont(jj2),phiav,yav));
461         } 
462       /** on permutte les points du rectangle **/
463       l1=i1; k1= jj1;
464       i1=i2;jj1=jj2;i2=i3;jj2=jj3;i3=i4;jj3=jj4;i4=l1;jj4=k1;
465       phi=phi1; phi1=phi2;phi2=phi3;phi3=phi4;phi4=phi;
466     }
467   (*func)(1,Cont,f_intercept(0.0,phi1,x_cont(i1),phi2,x_cont(i2)),
468             f_intercept(0.0,phi1,y_cont(jj1),phi2,y_cont(jj2)));
469   if ( qq==1 && bdyp(i1,jj1) && bdyp(i2,jj2)) *zds = 1 ;
470   if ( revflag == 1  &&  ! oddp (i) )  i = i+2;
471   return ( 1 + (  ( i + ent + 2) % 4));
472 }
473
474 /*--------------------------------------------------------------
475  * Following code is used to store the current level curves as 
476  * double in order to access to the stored data at Scilab level 
477  *----------------------------------------------------------------*/
478
479 static int last=-1;
480 static int count=0; 
481  
482 /** used to bring back data to Scilab Stack **/
483
484 int C2F(getconts)(double **x, double **y, int *mm, int *n)
485 {
486   *x = Gxcont;
487   *y = Gycont;
488   *mm= 1;
489   *n= Gcont_size;
490   return 0;
491 }
492
493 static void GContStore2(int ival, double Cont, double xncont, double yncont)
494 {
495   int n;
496   if ( ival == 0) 
497     {
498       /* Here : ival == 0 means stop the current level curve and 
499        * store data at the end but do reset Gcont_size to zero 
500        */
501       n= Gcont_size + 2;
502       if ( Gxcont == NULL )
503       {
504         Gxcont = MALLOC( n * sizeof(double) ) ;
505         Gycont = MALLOC( n * sizeof(double) ) ;
506       }
507       else
508       {
509         Gxcont = REALLOC( Gxcont, n * sizeof(double) ) ;
510         Gycont = REALLOC( Gycont, n * sizeof(double) ) ;
511       }
512       if ( (Gxcont == NULL) && n != 0) return ; 
513       if ( (Gycont == NULL) && n != 0) return ;
514       Gxcont[Gcont_size] = Cont;
515       if ( last != -1 && last < n ) Gycont[last]= count;
516       last = Gcont_size;
517       Gcont_size++;
518       count = 0;
519     }
520   else
521   {
522     n = Gcont_size + 1 ;
523     if ( Gxcont == NULL )
524     {
525       Gxcont = MALLOC( n * sizeof(double) ) ;
526       Gycont = MALLOC( n * sizeof(double) ) ;
527     }
528     else
529     {
530       Gxcont = REALLOC( Gxcont, n * sizeof(double) ) ;
531       Gycont = REALLOC( Gycont, n * sizeof(double) ) ;
532     }
533     if ( (Gxcont == NULL) && n != 0) return ; 
534     if ( (Gycont == NULL) && n != 0) return ;
535   }
536   Gxcont[Gcont_size]=xncont;
537   Gycont[Gcont_size++]=yncont;
538   count++;
539 }
540
541 static void GContStore2Last(void)
542 {
543   if ( last != -1 ) Gycont[last]= count;
544 }
545
546
547
548
549
550
551
552
553
554
555