geom.c
Go to the documentation of this file.
1 #include "geom.h"
2 
3 #include "error.h"
4 #include "defines.h"
5 
6 #include <math.h>
7 
18 boolean pointLeftOfLine(double px,double py,
19  double x1,double y1,double x2,double y2)
20 {
21  // |x1 y1 1|
22  // det = |x2 y2 1| >= 0
23  // |px py 1|
24  double det;
25 
26  ROUNDDOWN;
27  det=(x2*py-px*y2)-(x1*py-px*y1)+(x1*y2-x2*y1);
28  ROUNDNEAR;
29 
30  return(det>=0);
31 }
32 
33 
34 boolean pointRightOfLine(double px,double py,
35  double x1,double y1,double x2,double y2)
36 {
37  // |x1 y1 1|
38  // det = |x2 y2 1| <= 0
39  // |px py 1|
40  double det;
41 
42  ROUNDUP;
43  det=(x2*py-px*y2)-(x1*py-px*y1)+(x1*y2-x2*y1);
44  ROUNDNEAR;
45 
46  return (det<=0);
47 }
48 
49 /*
50  * Computes two intervals that delimite the bounding box for a set of 'n' coordinates
51  * in the x,y plane
52  */
53 void ComputeBoundingBox(unsigned int n,double *x,double *y,
54  Tinterval *ix,Tinterval *iy)
55 {
56  unsigned int i;
57  double x_min,x_max;
58  double y_min,y_max;
59 
60  x_min=x_max=x[0];
61  y_min=y_max=y[0];
62  for(i=1;i<n;i++)
63  {
64  if (x[i]>x_max) x_max=x[i];
65  else {if (x[i]<x_min) x_min=x[i];}
66  if (y[i]>y_max) y_max=y[i];
67  else {if (y[i]<y_min) y_min=y[i];}
68  }
69  NewInterval(x_min,x_max,ix);
70  NewInterval(y_min,y_max,iy);
71 }
72 
73 void ComputeBoundingBox3d(unsigned int n,double *x,double *y,double *z,
74  Tinterval *ix,Tinterval *iy,Tinterval *iz)
75 {
76  unsigned int i;
77  double x_min,x_max;
78  double y_min,y_max;
79  double z_min,z_max;
80 
81  x_min=x_max=x[0];
82  y_min=y_max=y[0];
83  z_min=z_max=z[0];
84  for(i=1;i<n;i++)
85  {
86  if (x[i]>x_max) x_max=x[i];
87  else {if (x[i]<x_min) x_min=x[i];}
88  if (y[i]>y_max) y_max=y[i];
89  else {if (y[i]<y_min) y_min=y[i];}
90  if (z[i]>z_max) z_max=z[i];
91  else {if (z[i]<z_min) z_min=z[i];}
92  }
93  NewInterval(x_min,x_max,ix);
94  NewInterval(y_min,y_max,iy);
95  NewInterval(z_min,z_max,iz);
96 }
97 
99  double *x,double *y,
100  double tol,
101  unsigned int *n,
102  double *xi,double *yi)
103 {
104  Tinterval x0,y0,vx,vy,a,b,c,d,t,pos,sqrtD,e,f,l,vl,s;
105 
106  NewInterval(x[0],x[0],&x0);
107  NewInterval(y[0],y[0],&y0);
108 
109  IntervalOffset(&x0,-x[1],&vx);
110  IntervalInvert(&vx,&vx);
111  IntervalOffset(&y0,-y[1],&vy);
112  IntervalInvert(&vy,&vy);
113 
114  IntervalPow(&vx,2,&a);
115  IntervalPow(&vy,2,&t);
116  IntervalAdd(&a,&t,&a);
117 
118  IntervalScale(&a,2,&f);
119 
120  if (IsInside(0,0,&f))
121  Error("SegmentCircleIntersection with a 0-length segment");
122 
123  IntervalProduct(&x0,&vx,&b);
124  IntervalProduct(&y0,&vy,&t);
125  IntervalAdd(&b,&t,&b);
126  IntervalScale(&b,2,&b);
127 
128  IntervalPow(&x0,2,&c);
129  IntervalPow(&y0,2,&t);
130  IntervalAdd(&c,&t,&c);
131  IntervalOffset(&c,-r2,&c);
132 
133  IntervalPow(&b,2,&d);
134  IntervalProduct(&a,&c,&t);
135  IntervalScale(&t,4,&t);
136  IntervalSubstract(&d,&t,&d);
137 
138  NewInterval(0,INF,&pos);
139  Intersection(&d,&pos,&d);
140 
141  if (EmptyInterval(&d))
142  *n=0;
143  else
144  {
145  NewInterval(0-tol,1+tol,&vl);
146 
147  IntervalSqrt(&d,&sqrtD);
148 
149  *n=0;
150  IntervalSubstract(&sqrtD,&b,&e);
151  IntervalDivision(&e,&f,&l);
152  Intersection(&l,&vl,&l);
153  if (!EmptyInterval(&l))
154  {
155  IntervalProduct(&l,&vx,&s);
156  IntervalAdd(&x0,&s,&s);
157  xi[0]=IntervalCenter(&s);
158 
159  IntervalProduct(&l,&vy,&s);
160  IntervalAdd(&y0,&s,&s);
161  yi[0]=IntervalCenter(&s);
162 
163  (*n)++;
164  }
165 
166  IntervalInvert(&sqrtD,&sqrtD);
167  IntervalSubstract(&sqrtD,&b,&e);
168  IntervalDivision(&e,&f,&l);
169  Intersection(&l,&vl,&l);
170  if (!EmptyInterval(&l))
171  {
172  IntervalProduct(&l,&vx,&s);
173  IntervalAdd(&x0,&s,&s);
174  xi[*n]=IntervalCenter(&s);
175 
176  IntervalProduct(&l,&vy,&s);
177  IntervalAdd(&y0,&s,&s);
178  yi[*n]=IntervalCenter(&s);
179 
180  (*n)++;
181  }
182  }
183 }
184 
186  double *x,double *y,double *z,
187  double tol,
188  unsigned int *n,
189  double *xi,double *yi,double *zi)
190 {
191  Tinterval x0,y0,z0,vx,vy,vz,a,b,c,d,t,pos,sqrtD,e,f,l,vl,s;
192 
193  NewInterval(x[0],x[0],&x0);
194  NewInterval(y[0],y[0],&y0);
195  NewInterval(z[0],z[0],&z0);
196 
197  IntervalOffset(&x0,-x[1],&vx);
198  IntervalInvert(&vx,&vx);
199  IntervalOffset(&y0,-y[1],&vy);
200  IntervalInvert(&vy,&vy);
201  IntervalOffset(&z0,-z[1],&vz);
202  IntervalInvert(&vz,&vz);
203 
204  IntervalPow(&vx,2,&a);
205  IntervalPow(&vy,2,&t);
206  IntervalAdd(&a,&t,&a);
207  IntervalPow(&vz,2,&t);
208  IntervalAdd(&a,&t,&a);
209 
210  IntervalScale(&a,2,&f);
211 
212  if (IsInside(0,0,&f))
213  Error("SegmentCircleIntersection with a 0-length segment");
214 
215  IntervalProduct(&x0,&vx,&b);
216  IntervalProduct(&y0,&vy,&t);
217  IntervalAdd(&b,&t,&b);
218  IntervalProduct(&z0,&vz,&t);
219  IntervalAdd(&b,&t,&b);
220  IntervalScale(&b,2,&b);
221 
222  IntervalPow(&x0,2,&c);
223  IntervalPow(&y0,2,&t);
224  IntervalAdd(&c,&t,&c);
225  IntervalPow(&z0,2,&t);
226  IntervalAdd(&c,&t,&c);
227  IntervalOffset(&c,-r2,&c);
228 
229  IntervalPow(&b,2,&d);
230  IntervalProduct(&a,&c,&t);
231  IntervalScale(&t,4,&t);
232  IntervalSubstract(&d,&t,&d);
233 
234  NewInterval(0,INF,&pos);
235  Intersection(&d,&pos,&d);
236 
237  if (EmptyInterval(&d))
238  *n=0;
239  else
240  {
241  NewInterval(0-tol,1+tol,&vl);
242 
243  IntervalSqrt(&d,&sqrtD);
244 
245  *n=0;
246  IntervalSubstract(&sqrtD,&b,&e);
247  IntervalDivision(&e,&f,&l);
248  Intersection(&l,&vl,&l);
249  if (!EmptyInterval(&l))
250  {
251  IntervalProduct(&l,&vx,&s);
252  IntervalAdd(&x0,&s,&s);
253  xi[0]=IntervalCenter(&s);
254 
255  IntervalProduct(&l,&vy,&s);
256  IntervalAdd(&y0,&s,&s);
257  yi[0]=IntervalCenter(&s);
258 
259  IntervalProduct(&l,&vz,&s);
260  IntervalAdd(&z0,&s,&s);
261  zi[0]=IntervalCenter(&s);
262 
263  (*n)++;
264  }
265 
266  IntervalInvert(&sqrtD,&sqrtD);
267  IntervalSubstract(&sqrtD,&b,&e);
268  IntervalDivision(&e,&f,&l);
269  Intersection(&l,&vl,&l);
270  if (!EmptyInterval(&l))
271  {
272  IntervalProduct(&l,&vx,&s);
273  IntervalAdd(&x0,&s,&s);
274  xi[*n]=IntervalCenter(&s);
275 
276  IntervalProduct(&l,&vy,&s);
277  IntervalAdd(&y0,&s,&s);
278  yi[*n]=IntervalCenter(&s);
279 
280  IntervalProduct(&l,&vz,&s);
281  IntervalAdd(&z0,&s,&s);
282  zi[*n]=IntervalCenter(&s);
283 
284  (*n)++;
285  }
286  }
287 }
288 
289 void Line2Points(double *x,double *y,boolean normalized,double *c)
290 {
291  double vx,vy;
292  double n;
293 
294  vx=x[1]-x[0];
295  vy=y[1]-y[0];
296 
297  n=sqrt(vx*vx+vy*vy);
298  if (n<ZERO)
299  Error("Punctual segment in Line2Points");
300 
301  c[0]= vy;
302  c[1]=-vx;
303  c[2]=-x[0]*vy+y[0]*vx;
304 
305  if (normalized)
306  {
307  c[0]/=n;
308  c[1]/=n;
309  c[2]/=n;
310  }
311 }
312 
313 void Plane3Points(double *x,double *y,double *z,
314  boolean normalized,double *c)
315 {
316  double v1[3],v2[3];
317  double n;
318 
319  v1[0]=x[1]-x[0];
320  v1[1]=y[1]-y[0];
321  v1[2]=z[1]-z[0];
322 
323  v2[0]=x[2]-x[0];
324  v2[1]=y[2]-y[0];
325  v2[2]=z[2]-z[0];
326 
327  CrossProduct(v1,v2,c);
328 
329  n=sqrt(c[0]*c[0]+c[1]*c[1]+c[2]*c[2]);
330 
331  if (n<ZERO)
332  Error("Aligned points in Plane3Points");
333 
334  if (normalized)
335  {
336  c[0]/=n;
337  c[1]/=n;
338  c[2]/=n;
339  }
340 
341  c[3]=-(c[0]*x[0]+c[1]*y[0]+c[2]*z[0]);
342 }
343 
344 /*
345  * Performs the Clipping between a rectangle and a circle.
346  */
347 boolean RectangleCircleClipping(double r2,
348  Tinterval *x,Tinterval *y,
349  Tinterval *x_new,Tinterval *y_new)
350 {
351  Tinterval d,md,a,b,diameter,rad2;
352  boolean full;
353  double r;
354 
355  /* First, both x,y must be inside the square +/-r
356  This is basically to deal with infinite ranges */
357  ROUNDUP;
358  r=sqrt(r2);
359  ROUNDNEAR;
360  NewInterval(-r,+r,&diameter);
361  full=((Intersection(x,&diameter,x_new))&&
362  (Intersection(y,&diameter,y_new)));
363 
364  if (full)
365  {
366  NewInterval(r2-ZERO,r2+ZERO,&rad2);
367  /* y = sqrt(r2 - x^2) */
368  IntervalPow(x_new,2,&d);
369  IntervalInvert(&d,&d);
370  IntervalAdd(&d,&rad2,&d);
371  IntervalSqrt(&d,&d);
372  IntervalInvert(&d,&md);
373 
374  Intersection(y_new,&d,&a);
375  Intersection(y_new,&md,&b);
376 
377  full=Union(&a,&b,y_new);
378 
379  if (full)
380  {
381  /* x = sqrt(r2 - y_new^2) */
382  IntervalPow(y_new,2,&d);
383  IntervalInvert(&d,&d);
384  IntervalAdd(&d,&rad2,&d);
385  IntervalSqrt(&d,&d);
386  IntervalInvert(&d,&md);
387 
388  Intersection(x_new,&d,&a);
389  Intersection(x_new,&md,&b);
390 
391  full=Union(&a,&b,x_new);
392  }
393  }
394  return(full);
395 }
396 
397 
398 /*
399  * Performs the Clipping between a 3d box and a sphere.
400  */
401 boolean BoxSphereClipping(double r2,
402  Tinterval *x,Tinterval *y,Tinterval *z,
403  Tinterval *x_new,Tinterval *y_new,Tinterval *z_new)
404 {
405  Tinterval d,md,a,b,rad2,diameter;
406  boolean full;
407  double r;
408 
409  /* First, both x,y must be inside the square +/-r
410  This is basically to deal with infinite ranges */
411  ROUNDUP;
412  r=sqrt(r2);
413  ROUNDNEAR;
414  NewInterval(-r,+r,&diameter);
415  full=((Intersection(x,&diameter,x_new))&&
416  (Intersection(y,&diameter,y_new))&&
417  (Intersection(z,&diameter,z_new)));
418 
419  if (full)
420  {
421  /* z = sqrt(r2 - x^2 - y^2) */
422  NewInterval(r2-ZERO,r2+ZERO,&rad2);
423  IntervalPow(x_new,2,&a);
424  IntervalPow(y_new,2,&b);
425  IntervalAdd(&a,&b,&d);
426  IntervalInvert(&d,&d);
427  IntervalAdd(&d,&rad2,&d);
428  IntervalSqrt(&d,&d);
429  IntervalInvert(&d,&md);
430 
431  Intersection(z_new,&d,&a);
432  Intersection(z_new,&md,&b);
433  full=Union(&a,&b,z_new);
434 
435  if (full)
436  {
437  /* y = sqrt(r2 - x^2 - z_new^2) */
438  IntervalPow(x_new,2,&a);
439  IntervalPow(z_new,2,&b);
440  IntervalAdd(&a,&b,&d);
441  IntervalInvert(&d,&d);
442  IntervalAdd(&d,&rad2,&d);
443  IntervalSqrt(&d,&d);
444  IntervalInvert(&d,&md);
445 
446  Intersection(y_new,&d,&a);
447  Intersection(y_new,&md,&b);
448  full=Union(&a,&b,y_new);
449 
450  if (full)
451  {
452  /* x = sqrt(r2 - y_new^2 - z_new^2) */
453  IntervalPow(y_new,2,&a);
454  IntervalPow(z_new,2,&b);
455  IntervalAdd(&a,&b,&d);
456  IntervalInvert(&d,&d);
457  IntervalAdd(&d,&rad2,&d);
458  IntervalSqrt(&d,&d);
459  IntervalInvert(&d,&md);
460 
461  Intersection(x_new,&d,&a);
462  Intersection(x_new,&md,&b);
463  full=Union(&a,&b,x_new);
464  }
465  }
466  }
467  return(full);
468 }
469 
470 /*
471  * Clipping of a 2d box with the function x^2 + \alpha y= \beta
472 */
473 boolean RectangleParabolaClipping(Tinterval *x,double alpha,double beta,Tinterval *y,
474  Tinterval *x_new,Tinterval *y_new)
475 {
476  Tinterval d,md,a,b,s;
477  boolean full;
478  double s1,s2;
479 
480  /* y = (1/alpha) *(beta - x^2) */
481  ROUNDDOWN;
482  s1=1/alpha;
483  ROUNDUP;
484  s2=1/alpha;
485  ROUNDNEAR;
486  NewInterval(s1,s2,&s);
487 
488  IntervalPow(x,2,&d);
489  IntervalInvert(&d,&d);
490  IntervalOffset(&d,beta,&d);
491  IntervalProduct(&d,&s,&d);
492 
493  full=Intersection(y,&d,y_new);
494 
495  if (full)
496  {
497  /* x = +/- sqrt(beta-alpha * y_new)*/
498  IntervalScale(y_new,alpha,&d);
499  IntervalInvert(&d,&d);
500  IntervalOffset(&d,beta,&d);
501 
502  IntervalSqrt(&d,&d);
503  IntervalInvert(&d,&md);
504 
505  Intersection(x,&d,&a);
506  Intersection(x,&md,&b);
507  full=Union(&a,&b,x_new);
508  }
509 
510  return(full);
511 }
512 
513 /*
514  A saddle is a function of the form
515  x * y + \alpha z = \beta
516  This implements the clipping between a saddle and a 3d box
517 */
518 boolean BoxSaddleClipping(Tinterval *x,Tinterval *y,double alpha,double beta,Tinterval *z,
519  Tinterval *x_new,Tinterval *y_new,Tinterval *z_new)
520 
521 {
522  Tinterval d,s;
523  boolean full;
524  double s1,s2;
525 
526  /* z=(1/alpha)*(beta-x*y) */
527  ROUNDDOWN;
528  s1=1/alpha;
529  ROUNDUP;
530  s2=1/alpha;
531  ROUNDNEAR;
532  NewInterval(s1,s2,&s);
533 
534  IntervalProduct(x,y,&d);
535  IntervalInvert(&d,&d);
536  IntervalOffset(&d,beta,&d);
537  IntervalProduct(&s,&d,&d);
538 
539  full=Intersection(z,&d,z_new);
540 
541  if (full)
542  {
543  /* x=(beta-alpha*z_new)/y */
544  if (IsInside(0.0,0.0,y))
545  CopyInterval(x_new,x);
546  else
547  {
548  IntervalScale(z_new,-alpha,&d);
549  IntervalOffset(&d,beta,&d);
550  IntervalDivision(&d,y,&d);
551  full=Intersection(x,&d,x_new);
552  }
553 
554  if (full)
555  {
556  /* y=(beta-alpha*z_new)/x_new */
557  if (IsInside(0.0,0.0,x_new))
558  CopyInterval(y_new,y);
559  else
560  {
561  IntervalScale(z_new,-alpha,&d);
562  IntervalOffset(&d,beta,&d);
563  IntervalDivision(&d,x_new,&d);
564  full=Intersection(y,&d,y_new);
565  }
566  }
567  }
568  return(full);
569 }
570 
571 void DefineNormalVector(double *v,double *n)
572 {
573  unsigned int i;
574  double norm;
575 
576  /* First check if we have a trivial case (X,Y,Z vector) */
577  if ((n[1]==0.0)&&(n[2]==0.0)) // X-vector
578  {
579  v[0]=0.0;v[1]=1.0;v[2]=0.0; // Y
580  }
581  else
582  {
583  if ((n[0]==0.0)&&(n[2]==0.0)) // Y-vector
584  {
585  v[0]=-1.0;v[1]=0.0;v[2]=0.0; // -X (Stay in the XY plane just in case)
586  }
587  else
588  {
589  if ((n[0]==0.0)&&(n[1]==0.0)) // Z-vector (already out of the XY plane)
590  {
591  v[0]=1.0;v[1]=0.0;v[2]=0.0; // X
592  }
593  else
594  {
595  if (n[2]==0.0) /* A vector in the X-Y plane */
596  {
597  v[0]=-n[1]; /* x=-y */
598  v[1]= n[0]; /* y= x */
599  v[2]=0.0;
600  }
601  else
602  {
603  /* The general case */
604  if (n[0]!=0.0)
605  {
606  v[1]=v[2]=1.0;
607  v[0]=-(n[1]+n[2])/n[0];
608  }
609  else
610  {
611  if (n[1]!=0.0)
612  {
613  v[0]=v[2]=1.0;
614  v[1]=-(n[0]+n[2])/n[1];
615  }
616  else
617  {
618  if (n[2]!=0.0)
619  {
620  v[0]=v[1]=1.0;
621  v[2]=-(n[0]+n[1])/n[2];
622  }
623  else
624  Error("Null input vector in DefineNormalVector");
625  }
626  }
627  }
628  norm=0.0;
629  for(i=0;i<3;i++)
630  norm=norm+v[i]*v[i];
631  norm=sqrt(norm);
632  for(i=0;i<3;i++)
633  v[i]/=norm;
634  }
635  }
636  }
637 }
638 
639 void CrossProduct(double *v1,double *v2,double *v3)
640 {
641  double a,b,c;
642 
643  a= (v1[1]*v2[2])-(v1[2]*v2[1]);
644  b=-(v1[0]*v2[2])+(v1[2]*v2[0]);
645  c= (v1[0]*v2[1])-(v1[1]*v2[0]);
646 
647  v3[0]=a;
648  v3[1]=b;
649  v3[2]=c;
650 }
651 
652 double DotProduct(double *v1,double *v2)
653 {
654  return(v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]);
655 }
656 
657 boolean NeighbouringTriangles(unsigned int *v1,unsigned int *v2)
658 {
659  unsigned int i;
660  unsigned int ns;
661 
662  ns=0;
663  for(i=0;i<3;i++)
664  {
665  if ((v1[i]==v2[0])||(v1[i]==v2[1])||(v1[i]==v2[2]))
666  ns++;
667  }
668  return(ns==2);
669 }
670 
671 boolean SameTriangle(unsigned int *v1,unsigned int *v2)
672 {
673  unsigned int i;
674  boolean same;
675 
676  same=TRUE;
677  for(i=0;((same)&&(i<3));i++)
678  same=((v1[i]==v2[0])||(v1[i]==v2[1])||(v1[i]==v2[2]));
679  return(same);
680 }
681 
682 double BallVolume(unsigned int n,double r)
683 {
684  unsigned int i,s;
685  double v;
686 
687  if ((n%2)==1)
688  {
689  // 1,3,5,...
690  s=3;
691  v=2*r;
692  }
693  else
694  {
695  // 0,2,4,....
696  s=2;
697  v=1;
698  }
699  /* V_n = (2*pi*R^2/n) * V_{n-2} */
700  for(i=s;i<=n;i+=2)
701  v*=(2.0*M_PI*r*r/(double)i);
702 
703  return(v);
704 }
705 
706 void FirstCombination(unsigned int n,unsigned int m,unsigned int *cb)
707 {
708  unsigned int i;
709 
710  for(i=0;i<n;i++)
711  cb[i]=0;
712 }
713 
714 boolean NextCombination(unsigned int n,unsigned int m,unsigned int *cb)
715 {
716  unsigned int i;
717  boolean end;
718 
719  end=FALSE;
720  for(i=0;((!end)&&(i<n));i++)
721  {
722  if (cb[i]<m)
723  {
724  cb[i]++;
725  end=TRUE;
726  }
727  else
728  cb[i]=0;
729  }
730 
731  return(end);
732 }
733 
734 double Det2x2(double *c1,double *c2)
735 {
736  return((c1[0]*c2[1])-(c1[1]*c2[0]));
737 }
738 
739 double Det3x3(double *c1,double *c2,double *c3)
740 {
741  return((c1[0]*c2[1]*c3[2]+c1[2]*c2[3]*c3[1]+c2[1]*c3[2]*c1[3])-
742  (c1[2]*c2[1]*c3[0]+c1[1]*c2[0]*c3[3]+c1[0]*c2[3]*c3[1]));
743 }
744 
745 void PrintReal(FILE *f,double r)
746 {
747  if (fabs(r)<ZERO)
748  fprintf(f,"0");
749  else
750  {
751  if (fabs(1-r)<ZERO)
752  fprintf(f,"1");
753  else
754  {
755  if (fabs(-1-r)<ZERO)
756  fprintf(f,"-1");
757  else
758  fprintf(f,"%.16f",r);
759  }
760  }
761 }
762 
763 void Print3Reals(FILE *f,double r1,double r2,double r3)
764 {
765  PrintReal(f,r1);
766  fprintf(f,",");
767  PrintReal(f,r2);
768  fprintf(f,",");
769  PrintReal(f,r3);
770 }
771 
772 
773 void PrintAngle(FILE *f,double a)
774 {
775  if (fabs(a)<ZERO)
776  fprintf(f,"0");
777  else
778  {
779  if (fabs(M_PI-a)<ZERO)
780  fprintf(f,"pi");
781  else
782  {
783  if (fabs(-M_PI-a)<ZERO)
784  fprintf(f,"-pi");
785  else
786  {
787  if (fabs(M_PI-a)<ZERO)
788  fprintf(f,"pi/2");
789  else
790  {
791  if (fabs(-M_PI-a)<ZERO)
792  fprintf(f,"-pi/2");
793  else
794  fprintf(f,"%.16f",a);
795  }
796  }
797  }
798  }
799 }
boolean Intersection(Tinterval *i1, Tinterval *i2, Tinterval *i_out)
Computes the intersection of two intervals.
Definition: interval.c:285
Definition of basic functions.
#define FALSE
FALSE.
Definition: boolean.h:30
boolean pointLeftOfLine(double px, double py, double x1, double y1, double x2, double y2)
Determines if a point is at the left of a vector.
Definition: geom.c:18
#define ROUNDDOWN
Sets the floating point operations in rounding down mode.
Definition: defines.h:219
void FirstCombination(unsigned int n, unsigned int m, unsigned int *cb)
Initializes a combination.
Definition: geom.c:706
double Det3x3(double *c1, double *c2, double *c3)
Determinant of a 3x3 matrix.
Definition: geom.c:739
void IntervalAdd(Tinterval *i1, Tinterval *i2, Tinterval *i_out)
Addition of two intervals.
Definition: interval.c:423
void IntervalInvert(Tinterval *i, Tinterval *i_out)
Change of sign of a given interval.
Definition: interval.c:461
#define TRUE
TRUE.
Definition: boolean.h:21
double Det2x2(double *c1, double *c2)
Determinant of a 2x2 matrix.
Definition: geom.c:734
void Error(const char *s)
General error function.
Definition: error.c:80
void SegmentSphereIntersection(double r2, double *x, double *y, double *z, double tol, unsigned int *n, double *xi, double *yi, double *zi)
Intersects a segment and a sphere.
Definition: geom.c:185
boolean RectangleParabolaClipping(Tinterval *x, double alpha, double beta, Tinterval *y, Tinterval *x_new, Tinterval *y_new)
Clips a 2D box and a scaled parabola.
Definition: geom.c:473
void ComputeBoundingBox(unsigned int n, double *x, double *y, Tinterval *ix, Tinterval *iy)
Determines the axis aligned bounding box of a set of points in R^2.
Definition: geom.c:53
#define ZERO
Floating point operations giving a value below this constant (in absolute value) are considered 0...
Definition: defines.h:37
boolean NeighbouringTriangles(unsigned int *v1, unsigned int *v2)
Determines if two sets of vertices define a neighbouring triangle.
Definition: geom.c:657
void CopyInterval(Tinterval *i_dst, Tinterval *i_org)
Copy constructor.
Definition: interval.c:59
boolean BoxSaddleClipping(Tinterval *x, Tinterval *y, double alpha, double beta, Tinterval *z, Tinterval *x_new, Tinterval *y_new, Tinterval *z_new)
Clips a 3D box and a scaled saddle.
Definition: geom.c:518
boolean NextCombination(unsigned int n, unsigned int m, unsigned int *cb)
Moves to the next combination.
Definition: geom.c:714
Error and warning functions.
boolean IntervalSqrt(Tinterval *i, Tinterval *i_out)
Interval square root.
Definition: interval.c:530
double DotProduct(double *v1, double *v2)
Computes the dot product of two 3d vectors.
Definition: geom.c:652
Definitions of constants and macros used in several parts of the cuik library.
void IntervalSubstract(Tinterval *i1, Tinterval *i2, Tinterval *i_out)
Substraction of two intervals.
Definition: interval.c:442
void SegmentCircleIntersection(double r2, double *x, double *y, double tol, unsigned int *n, double *xi, double *yi)
Intersects a segment and a circle.
Definition: geom.c:98
void PrintAngle(FILE *f, double a)
Pretty print an angle.
Definition: geom.c:773
void ComputeBoundingBox3d(unsigned int n, double *x, double *y, double *z, Tinterval *ix, Tinterval *iy, Tinterval *iz)
Determines the axis aligned bounding box of a set of points in R^3.
Definition: geom.c:73
#define ROUNDNEAR
Sets the floating point operations in round near mode.
Definition: defines.h:225
void PrintReal(FILE *f, double r)
Pretty print a real number.
Definition: geom.c:745
void Plane3Points(double *x, double *y, double *z, boolean normalized, double *c)
Defines a plane given 3 points.
Definition: geom.c:313
boolean pointRightOfLine(double px, double py, double x1, double y1, double x2, double y2)
Determines if a point is at the right of a vector.
Definition: geom.c:34
#define M_PI
Pi.
Definition: defines.h:83
void IntervalDivision(Tinterval *num, Tinterval *den, Tinterval *i_out)
Interval division.
Definition: interval.c:556
boolean Union(Tinterval *i1, Tinterval *i2, Tinterval *i_out)
Computes the union of two intervals.
Definition: interval.c:297
void IntervalScale(Tinterval *i1, double e, Tinterval *i_out)
Scales an interval.
Definition: interval.c:360
double IntervalCenter(Tinterval *i)
Gets the interval center.
Definition: interval.c:129
boolean RectangleCircleClipping(double r2, Tinterval *x, Tinterval *y, Tinterval *x_new, Tinterval *y_new)
Clips a rectangle and a circle.
Definition: geom.c:347
#define ROUNDUP
Sets the floating point operations in rounding up mode.
Definition: defines.h:213
boolean EmptyInterval(Tinterval *i)
Checks if a given interval is empty.
Definition: interval.c:335
boolean IsInside(double p, double tol, Tinterval *i)
Checks if a value is inside an interval.
Definition: interval.c:348
boolean SameTriangle(unsigned int *v1, unsigned int *v2)
Identifies triangles formed by the same set of vertices.
Definition: geom.c:671
void IntervalProduct(Tinterval *i1, Tinterval *i2, Tinterval *i_out)
Product of two intervals.
Definition: interval.c:389
void IntervalPow(Tinterval *i, unsigned int p, Tinterval *i_out)
Power of a given interval by a integer power factor.
Definition: interval.c:494
void NewInterval(double lower, double upper, Tinterval *i)
Constructor.
Definition: interval.c:47
#define INF
Infinite.
Definition: defines.h:70
void DefineNormalVector(double *v, double *n)
Defines a unitary vector normal to a given vector.
Definition: geom.c:571
boolean BoxSphereClipping(double r2, Tinterval *x, Tinterval *y, Tinterval *z, Tinterval *x_new, Tinterval *y_new, Tinterval *z_new)
Clips a 3D box and a sphere.
Definition: geom.c:401
Defines a interval.
Definition: interval.h:33
void CrossProduct(double *v1, double *v2, double *v3)
Computes the cross product of two 3d vectors.
Definition: geom.c:639
double BallVolume(unsigned int n, double r)
Volume of of a n-ball.
Definition: geom.c:682
void Line2Points(double *x, double *y, boolean normalized, double *c)
Defines a line given two points.
Definition: geom.c:289
void Print3Reals(FILE *f, double r1, double r2, double r3)
Pretty print three real number.
Definition: geom.c:763
void IntervalOffset(Tinterval *i, double offset, Tinterval *i_out)
Interval offset.
Definition: interval.c:627