Remove graphics.h file.
[scilab.git] / scilab / modules / graphics / src / c / Contour.c
1 /*------------------------------------------------------------------------
2  *    Graphic library
3  *    Copyright (C) 1998-2001 Enpc/Jean-Philippe Chancelier
4  *    jpc@cermics.enpc.fr 
5  --------------------------------------------------------------------------*/
6
7 #include <string.h> /* in case of dbmalloc use */
8 #include <stdio.h>
9 #include <math.h>
10 #include "math_graphics.h"
11 #include "PloEch.h"
12 #include "Xcall1.h"
13 #include "sciprint.h"
14 #include "MALLOC.h"
15 #include "Format.h"
16 #include "Contour.h"
17
18 #include "GetProperty.h"
19
20
21 #ifdef _MSC_VER
22 extern void Scistring (char *str);
23 #endif
24
25 typedef void (level_f) __PARAMS((integer ival, double Cont, double xncont,
26                                double yncont));
27 typedef void (*ptr_level_f) __PARAMS((integer ival, double Cont, double xncont,
28                                double yncont));
29
30
31 static void 
32 contourI __PARAMS((ptr_level_f,double *, double *, double *,
33                   double *, integer *, integer *, integer *));
34 static void
35 look __PARAMS((ptr_level_f, integer i, integer j, integer ib,
36               integer jb, integer qq,double Cont, integer style));
37
38 static integer ffnd __PARAMS((ptr_level_f,integer,integer,integer,integer,integer,
39                              integer,integer,integer,integer,integer,
40                              double,integer *));
41
42 static int Gcont_size = 0;
43
44 static void ContourTrace __PARAMS((double Cont, integer style));
45 static level_f ContStore, ContStore1, ContStore2,GContStore2;
46 static void GContStore2Last __PARAMS((void));
47 static double x_cont __PARAMS((integer i));
48 static double y_cont __PARAMS((integer i));
49 static double phi_cont __PARAMS((integer, integer)); 
50 static double f_intercept  __PARAMS((double, double, double, double, double ));
51 static integer not_same_sign  __PARAMS((double val1, double val2)); 
52 static integer get_itg_cont  __PARAMS((integer i, integer j)); 
53 static void inc_itg_cont  __PARAMS((integer i, integer j, integer val)); 
54 static integer oddp  __PARAMS((integer 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 #define HIDDENFRAMECOLOR 2L /* default color for hidden frame */
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 integer Gn1,Gn2;
71
72 static double * Gxcont = NULL ;
73 static double * Gycont = NULL ;
74
75 static integer * xcont = NULL ;
76 static integer * ycont = NULL ;
77
78 static void InitValues(double *x, double *y, double *z, integer n1, integer n2)
79 {
80   Gn1=n1;  Gn2=n2;  GX = x;  GY = y;  GZ = z;
81 }
82
83 /*--------return the  value of f for a pointeger on the grid-----*/
84
85 static double phi_cont(integer i, integer 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
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
100 static integer  bdyp(integer i, integer 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  integer *itg_cont, *xbd_cont,*ybd_cont;
108
109 static integer get_itg_cont(integer i, integer j)
110 {
111   return( itg_cont[i+Gn1*j]);
112 }
113
114 static void inc_itg_cont(integer i, integer j, integer val)
115 {
116   itg_cont[i+Gn1*j] += val;
117 }
118
119 static integer not_same_sign(double val1, double val2)
120 {
121   if ( ISNAN(val1) ==1 || ISNAN(val2) == 1) return(0);
122   /** 0.0 est consid\'er\'e comme positif **/
123   if ( val1 >= 0.0) 
124     {
125       if (val2 < 0.0) return(1) ; else return(0);}
126   else 
127     {
128       if ( val2 >= 0.0) return(1) ; else return(0);}
129 }
130
131 static integer oddp(integer i) { return( i == 1 || i ==3 );}
132
133 /*---------return the x-value of a grid point--------*/
134
135 static double x_cont(integer i) { return GX[i] ;}
136
137 /*---------return the y-value of a grid point --------*/
138
139 static double y_cont(integer i) { return GY[i] ;}
140
141
142
143 static double ZC=0.0;
144 static char   ContNumFormat[100];
145
146 /*--------------------------------------------------------------------
147 *  the level curve is crossing the segment (i,j) (ib,jb)
148 *  look store the level curve point and try to find the next segment to look at
149 *  Cont: value of f along the contour 
150 *  ncont: number of contour 
151 *  c: indice of the contour Cont 
152 *---------------------------------------------------------------------*/
153
154 static void look(ptr_level_f func, integer i, integer j, integer ib, integer jb, integer qq, double Cont, integer style)
155 {
156   integer ip,jp,im,jm,zds,ent=0,flag=0,wflag;
157   jp= j+1; ip= i+1; jm=j-1;im=i-1;
158   /*  on regarde comment est le segment de depart */
159   if  ( jb == jm)  flag = 1; 
160   else  { 
161     if ( ib == im ) flag = 2 ;
162     else  {
163       if ( jb == jp ) flag = 3 ;
164       else  if ( ib == ip ) flag = 4;}}
165   switch  (  flag)
166   {
167   case  1 :
168     if  (get_itg_cont(i,jm) > 1) return;
169     ent=1 ; /* le segment est vertical vers le bas */
170     /* Storing intersection point */
171     (*func)(0,Cont, x_cont(i), 
172       f_intercept(Cont,phi_cont(i,jm),
173       y_cont(jm),phi_cont(i,j),y_cont(j)));
174     break;
175   case 2 : 
176     if  (get_itg_cont(im,j) == 1 || get_itg_cont(im,j)==3 ) return;
177     ent=2 ; /* le segment est horizontal gauche */
178     /* Storing intersection point */
179     (*func)( 0,Cont,
180       f_intercept(Cont,phi_cont(im,j),
181       x_cont(im),phi_cont(i,j),x_cont(i)), y_cont(j));
182     break ; 
183   case 3 :
184     if  (get_itg_cont(i,j) > 1 ) return;
185     ent=3 ; /* le segment est vertical haut */
186     /* Storing intersection point */
187     (*func)(0,Cont,x_cont(i), f_intercept(Cont,phi_cont(i,j),
188       y_cont(j),phi_cont(i,jp),y_cont(jp)));
189     break ;
190   case 4 :
191     if  (get_itg_cont(i,j) == 1 || get_itg_cont(i,j)==3 ) return;
192     ent=4 ; /* le segment est horizontal droit */
193     /* Storing intersection point */
194     (*func)(0,Cont,f_intercept(Cont,phi_cont(i,j),
195       x_cont(i),phi_cont(ip,j),x_cont(ip)),
196       y_cont(j));
197     break;
198   default :
199     sciprint(" Error in case wrong value ");
200     break;
201   }
202   wflag=1;
203   while ( wflag) 
204   { 
205     jp= j+1; ip= i+1; jm=j-1;im=i-1;
206     switch  ( ent) 
207     {case 1 :
208     inc_itg_cont(i,jm,2L);
209     ent = ffnd(func,i,ip,ip,i,j,j,jm,jm,ent,qq,Cont,&zds);
210     /* on calcule le nouveau point, ent donne la 
211     direction du segment a explorer */
212     switch ( ent)
213     {
214     case -1: wflag=0; break;
215     case 1 : i=ip ; break ;
216     case 2 : i=ip;j=jm; break ;
217     }
218     break ;
219     case 2  :
220       inc_itg_cont(im,j,1L);
221       ent = ffnd(func,i,i,im,im,j,jm,jm,j,ent,qq,Cont,&zds);
222       switch ( ent)
223       { 
224       case -1: wflag=0; break;
225       case 2 : j = jm ;break ;
226       case  3  : i=im;j=jm; break ;
227       }
228       break ;
229     case 3 :
230       inc_itg_cont(i,j,2L);
231       ent = ffnd(func,i,im,im,i,j,j,jp,jp,ent,qq,Cont,&zds);
232       switch ( ent)
233       { 
234       case -1: wflag=0; break;
235       case 3 : i=im; break ;
236       case 4 : i=im;j=jp; break ;
237       }
238       break ;
239     case 4 :
240       inc_itg_cont(i,j,1L);
241       ent = ffnd(func,i,i,ip,ip,j,jp,jp,j,ent,qq,Cont,&zds);
242       switch ( ent)
243       {
244       case -1: wflag=0; break;
245       case 4 :j=jp;break ;
246       case 1 :i=ip;j=jp;break ;
247       }
248       break ;
249     }
250
251     /** new segment is on the boundary **/
252     if ( zds == 1) 
253     {
254       switch ( ent) 
255       {
256       case 1 : inc_itg_cont(i,(j-1),2L); break ; 
257       case 2 : inc_itg_cont(i-1,j,1L);  break ; 
258       case 3 : inc_itg_cont(i,j,2L); break ; 
259       case 4 : inc_itg_cont(i,j,1L); break ; 
260       }
261       /** we must quit the while loop **/
262       wflag = 0 ;
263     }
264     /**  init point was inside the domain */
265     if ( qq == 2) 
266     {
267       switch ( ent) 
268       {
269       case 1 : if  ( get_itg_cont (i,j-1)  > 1) wflag = 0 ; break ; 
270       case 2 : if  ( oddp(get_itg_cont(i-1,j))) wflag = 0 ; break ; 
271       case 3 : if  ( get_itg_cont(i,j) > 1)     wflag = 0 ; break ; 
272       case 4 : if  ( oddp(get_itg_cont(i,j)))   wflag = 0 ; break ; 
273       }
274     }
275   }
276   if ( func == GContStore2 ) 
277     GContStore2Last();
278   else 
279     ContourTrace(Cont,style);
280 }
281
282
283 /*-------------------------------------------------------
284 *  The function f is given on a grid and we want the level curves 
285 *  for the zCont[N[2]] values 
286 *  x : of size N[0] gives the x-values of the grid 
287 *  y : of size N[1] gives the y-values of the grid 
288 *  z : of size N[0]*N[1]  gives the f-values on the grid 
289 *  style: size ncont (=N[2]) or empty integer pointer 
290 *  gives the dash style for contour i
291 *-------------------------------------------------------*/
292
293 static void contourI(ptr_level_f func, double *x, double *y, double *z, double *zCont, integer *N, integer *style, integer *err)
294 {
295   int check = 1;
296   char *F;
297   integer n1,n2,ncont,i,c,j,k,n5;
298   integer stylec;
299   n1=N[0];n2=N[1];ncont=N[2];
300   F=getFPF();
301   if ( F[0] == '\0') 
302     ChoixFormatE1(ContNumFormat,zCont,N[2]);
303   InitValues(x,y,z,n1,n2);
304   n5 =  2*(n1)+2*(n2)-3;
305   /* Allocation */
306   Gcont_size = 0; /** initialize the array indices for storing contours **/
307   xbd_cont = MALLOC( n5 * sizeof(int) ) ;
308   ybd_cont = MALLOC( n5 * sizeof(int) ) ;
309   itg_cont = MALLOC( n1*n2 * sizeof(int) ) ;
310   if ( (xbd_cont == NULL) && n5 != 0) check= 0;
311   if ( (ybd_cont == NULL) && n5 != 0) check= 0;
312   if ( (itg_cont == NULL) && n1*n2 != 0) check= 0;
313   if ( check == 0) 
314   {
315     FREE( xbd_cont ) ;
316     FREE( ybd_cont ) ;
317     FREE( itg_cont ) ;
318     sciprint("contourI_: Running out of memory\n");
319     return;
320   }
321   /* just a parametrization of the boundary points */
322   for ( i = 0 ; i <  n2 ; i++)
323   {
324     ybd_cont[i] = i ;
325     xbd_cont[i] = 0 ;
326   }
327   for ( i = 1 ; i <  n1 ; i++)
328   {
329     ybd_cont[n2+i-1] = n2-1 ;
330     xbd_cont[n2+i-1] = i  ;
331   }
332   for ( i = n2-2;  i >= 0  ; i--)
333   {
334     ybd_cont[2*n2 +n1-3-i] = i ;
335     xbd_cont[2*n2 +n1-3-i] = n1-1  ;
336   }
337   for ( i = n1-2 ; i >= 0 ; i--)
338   {
339     ybd_cont[2*n2 +2*n1-4-i] = 0 ;
340     xbd_cont[2*n2 +2*n1-4-i] = i  ;
341   }
342   for ( c= 0 ; c < ncont ; c++)
343   {
344     stylec = ( style != (integer *) 0) ? style[c] : c;
345     /** itg-cont is a flag array to memorize checked parts of the grid **/
346     for ( i = 0 ; i < n1; i++)
347       for ( j =0 ; j < n2 ; j++)
348         itg_cont[i+n1*j]=0 ;
349     /** all the boundary segments **/
350     for ( k = 1 ; k < n5 ; k++)
351     { integer ib,jb;
352     i = xbd_cont[k] ; j = ybd_cont[k];
353     ib = xbd_cont[k-1] ; jb= ybd_cont[k-1];
354     if  (not_same_sign (phi_cont(i,j)-zCont[c] , 
355       phi_cont(ib,jb)-zCont[c]))
356       look(func,i,j,ib,jb,1L,zCont[c],stylec);
357     }
358     /** inside segments **/
359     for ( i = 1 ; i < n1-1; i++)
360       for ( j = 1 ; j < n2-1 ; j++)
361         if  (not_same_sign ( phi_cont(i,j)-zCont[c] , 
362           phi_cont(i, j-1)-zCont[c]))
363           look(func,i,j,i,j-1,2L,zCont[c],stylec);
364   }
365   FREE( xbd_cont ) ;
366   FREE( ybd_cont ) ;
367   FREE( itg_cont ) ;
368
369 }
370
371 int C2F(contourif)(double *x, double *y, double *z, integer *n1, integer *n2, integer *flagnz, integer *nz, double *zz, integer *style)
372 {
373   integer err=0;
374   static double *zconst;
375   double zmin,zmax;
376   integer 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           sciprint("Running out of memory\r\n");
386           return 0;
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 integer ffnd (ptr_level_f func, integer i1, integer i2, integer i3, integer i4, integer jj1, integer jj2, integer jj3, integer jj4, integer ent, integer qq, double Cont, integer *zds)
417 {
418   double phi1,phi2,phi3,phi4,xav,yav,phiav;
419   integer 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       integer 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     { integer 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  * Storing and tracing level curves 
476  *----------------------------------------------------------------*/
477
478 static integer cont_size ;
479
480 /* 
481  * Explicit drawing of the current level curve with a dash style 
482  * The curve level is also drawn as a string according to current 
483  * floating point format 
484  */
485
486 static void ContourTrace(double Cont, integer style)
487
488   char *F;
489   integer verbose=0 ,Dnarg,Dvalue[10];
490   integer close=0, flag=0, uc;
491   double angle=0.0;
492   char str[100];
493
494   C2F(dr)("xget","use color",&verbose,&uc,&Dnarg,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
495   if (uc) {
496     C2F(dr)("xget","color",&verbose,Dvalue,&Dnarg,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
497     C2F(dr)("xset","color",&style,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
498     C2F(dr)("xlines","void",&cont_size,xcont,ycont,&close,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
499     C2F(dr)("xset","color",Dvalue,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
500   }
501   else {
502     C2F(dr)("xget","line style",&verbose,Dvalue,&Dnarg,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
503     C2F(dr)("xset","line style",&style,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
504     C2F(dr)("xlines","void",&cont_size,xcont,ycont,&close,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
505     C2F(dr)("xset","line style",Dvalue,PI0,PI0,PI0,PI0,PI0,PD0,PD0,PD0,PD0,0L,0L);
506   }
507
508   F=getFPF();
509   if ( F[0] == '\0') 
510     sprintf(str,ContNumFormat,Cont);
511   else 
512     sprintf(str,F,Cont);
513   C2F(dr)("xstring",str, &xcont[cont_size / 2],&ycont[cont_size /2],
514           PI0,&flag,PI0,PI0, &angle,PD0,PD0,PD0,0L,0L);
515 }
516
517 /*--------------------------------------------------------------
518  * Following code is used to store the current level curves as 
519  * double in order to access to the stored data at Scilab level 
520  *----------------------------------------------------------------*/
521
522 static int last=-1;
523 static int count=0; 
524  
525 /** used to bring back data to Scilab Stack **/
526
527 int C2F(getconts)(double **x, double **y, integer *mm, integer *n)
528 {
529   *x = Gxcont;
530   *y = Gycont;
531   *mm= 1;
532   *n= Gcont_size;
533   return 0;
534 }
535
536 static void GContStore2(integer ival, double Cont, double xncont, double yncont)
537 {
538   int n;
539   if ( ival == 0) 
540     {
541       /* Here : ival == 0 means stop the current level curve and 
542        * store data at the end but do reset Gcont_size to zero 
543        */
544       n= Gcont_size + 2;
545       if ( Gxcont == NULL )
546       {
547         Gxcont = MALLOC( n * sizeof(double) ) ;
548         Gycont = MALLOC( n * sizeof(double) ) ;
549       }
550       else
551       {
552         Gxcont = REALLOC( Gxcont, n * sizeof(double) ) ;
553         Gycont = REALLOC( Gycont, n * sizeof(double) ) ;
554       }
555       if ( (Gxcont == NULL) && n != 0) return ; 
556       if ( (Gycont == NULL) && n != 0) return ;
557       Gxcont[Gcont_size] = Cont;
558       if ( last != -1 && last < n ) Gycont[last]= count;
559       last = Gcont_size;
560       Gcont_size++;
561       count = 0;
562     }
563   else
564   {
565     n = Gcont_size + 1 ;
566     if ( Gxcont == NULL )
567     {
568       Gxcont = MALLOC( n * sizeof(double) ) ;
569       Gycont = MALLOC( n * sizeof(double) ) ;
570     }
571     else
572     {
573       Gxcont = REALLOC( Gxcont, n * sizeof(double) ) ;
574       Gycont = REALLOC( Gycont, n * sizeof(double) ) ;
575     }
576     if ( (Gxcont == NULL) && n != 0) return ; 
577     if ( (Gycont == NULL) && n != 0) return ;
578   }
579   Gxcont[Gcont_size]=xncont;
580   Gycont[Gcont_size++]=yncont;
581   count++;
582 }
583
584 static void GContStore2Last(void)
585 {
586   if ( last != -1 ) Gycont[last]= count;
587 }
588
589
590
591
592
593
594
595
596
597
598