home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / ifs.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  34.4 KB  |  1,104 lines
  1. /*
  2.  * ifs.c
  3.  *
  4.  * Philippe Bekaert - bekaert@cs.kuleuven.ac.be - 19 may 93
  5.  */
  6. #include <stdlib.h>
  7. #include <math.h>
  8. #include "geom.h"
  9. #include "ifs.h"
  10. #include "bounds.h"
  11. #include "libcommon/sampling.h"
  12. #include "libshade/viewing.h"
  13. #include "libcommon/common.h"
  14.  
  15. static Methods *iIfsMethods = NULL, *iCIfsMethods = NULL;
  16. static char ifsName[] = "ifs", CifsName[] = "condensed ifs";
  17. unsigned long IfsTests, IfsHits, CIfsTests, CIfsHits;;
  18.  
  19. #define OPTIMAL_BV  /* for better bounding spheres */
  20.  
  21. /* return less than zero if first map has higher rang number than the second,
  22.  *      zero if two maps have equal rang numbers, 
  23.  *        greater than zero if first map has lower rang number than the second */
  24. static int rangcmp(map1, map2)
  25. IfsMap *map1, *map2;
  26. {
  27.       return map2->rang - map1->rang;
  28. }
  29.  
  30. /* computes radius of smallest sphere with centre (x,y,z) which encloses
  31.  * the n spheres with centre c[i] and radius r[i]: to be minimized */
  32. static Float REval(n, c, r, x, y, z)
  33. int n;
  34. Vector *c; Float *r;
  35. Float x, y, z;
  36. {
  37.     Vector C; 
  38.     Float R, R1;
  39.     int i;
  40.  
  41.     C.x = x; C.y = y; C.z = z;
  42.     R = VecDist(&C, &c[0]) + r[0];
  43.     for (i=1; i<n; i++)
  44.         if ((R1 = VecDist(&C, &c[i]) + r[i]) > R) R = R1;
  45.     return R;
  46. }
  47.  
  48. static Float SmallestBounding(n, c, r, C)
  49. int n;
  50. Vector *c; Float *r;
  51. Vector *C;
  52. {
  53.     Float R, R1, step;
  54.     int i, dir, go_on;
  55.  
  56. /* initial guess for center point: gravitation centre of centre points */
  57.     *C = c[0];
  58.     for (i=1; i<n; i++)
  59.         VecAdd(*C, c[i], C);
  60.     VecScale(1./(float)n, *C, C);
  61.  
  62. /* calculate radius and minimize */
  63.     R = REval(n, c, r, C->x, C->y, C->z);
  64.     step = R/8.;
  65. #ifdef OPTIMAL_BV
  66.     do {
  67.         dir = 0;
  68.         do {
  69.             go_on = FALSE;
  70.             if (dir != -1 && (R1 = REval(n, c, r, C->x + step, C->y, C->z)) < R) 
  71.                       { R = R1; C->x += step; dir = 1; go_on = TRUE;}
  72.             if (dir !=  1 && (R1 = REval(n, c, r, C->x - step, C->y, C->z)) < R) 
  73.                       { R = R1; C->x -= step; dir = -1; go_on = TRUE;}
  74.         } while (go_on); 
  75.  
  76.         dir = 0;
  77.         do {
  78.             go_on = FALSE;
  79.             if (dir != -1 && (R1 = REval(n, c, r, C->x, C->y + step, C->z)) < R) 
  80.                       { R = R1; C->y += step; dir = 1; go_on = TRUE;}
  81.             if (dir !=  1 && (R1 = REval(n, c, r, C->x, C->y - step, C->z)) < R) 
  82.                       { R = R1; C->y -= step; dir = -1; go_on = TRUE;}
  83.         } while (go_on); 
  84.  
  85.         dir = 0;
  86.         do {
  87.             go_on = FALSE;
  88.             if (dir != -1 && (R1 = REval(n, c, r, C->x, C->y, C->z + step)) < R) 
  89.                       { R = R1; C->z += step; dir = 1; go_on = TRUE;}
  90.             if (dir !=  1 && (R1 = REval(n, c, r, C->x, C->y, C->z - step)) < R) 
  91.                       { R = R1; C->z -= step; dir = -1; go_on = TRUE;}
  92.         } while (go_on); 
  93.  
  94.         step /= 2.;
  95.     } while (step > EPSILON);
  96. #endif /* OPTIMAL_BV */
  97.  
  98.     return R;
  99. }
  100.  
  101. static void IfsPartialBounding(nrmaps, maps, condensed, C, R, minsize, maxdepth)
  102. int nrmaps, condensed, maxdepth;
  103. IfsMap *maps;
  104. Vector *C;
  105. Float *R, minsize;
  106. {
  107.     int i, N, d;
  108.     Float *r, R1, delta, maxcontractivity;
  109.     Vector *c, C1;
  110.  
  111. /* largest Lipschitz constant of the maps */
  112.     for (i = 0, maxcontractivity = 0.; i < nrmaps; i++)
  113.         if (maps[i].contractivity > maxcontractivity) 
  114.             maxcontractivity = maps[i].contractivity;
  115.  
  116. /* for IFS without condensation: make new sphere which contains the images
  117.  * of this sphere and repeat as long the new sphere gets smaller. For IFS
  118.  * with condensation set: make new sphere which contains the old and it's
  119.  * images and repeat as long as it gets larger */
  120.     N = condensed ? nrmaps+1 : nrmaps;
  121.         c = (Vector *)Malloc(N * sizeof(Vector)); 
  122.         r =  (Float *)Malloc(N * sizeof( Float)); 
  123.  
  124.         for (delta = R1 = *R, C1 = *C, d = 0;
  125.          delta > EPSILON && delta >= minsize * maxcontractivity && d <= maxdepth; 
  126.          delta *= maxcontractivity, d++) {
  127.                 if (!condensed)
  128.                         if (R1 > *R) break; 
  129.                 *C = C1; *R = R1;
  130.  
  131.                 for (i=0; i<nrmaps; i++) {
  132.             c[i] = *C;
  133.                         PointTransform(&c[i], &maps[i].matrix);
  134.             r[i] = *R * maps[i].contractivity;
  135.         }
  136.                 if (condensed) {
  137.                         c[nrmaps] = *C;
  138.             r[nrmaps] = *R;
  139.         }
  140.  
  141.         R1 = SmallestBounding(N, c, r, &C1);
  142.         }
  143.     if (condensed) { *C = C1; *R = R1; }
  144.         free(c); free(r);
  145.  
  146.         *R += EPSILON;          /* to be on the safe side */
  147. }
  148.  
  149. /* Computes a sphere which is supposed to fit the attractor of the IFS quite tightly */
  150. static void IfsBoundingSphere(ifs, R, C)
  151. Ifs *ifs;
  152. Float *R;
  153. Vector *C;
  154. {
  155.         Float R1;
  156.         int i, n, condensed;
  157.     IfsMap *m;
  158.  
  159.     condensed = ifs->condensation ? TRUE : FALSE;
  160.         if (!condensed) {
  161. /* first guess: a sphere with centre in the origin which contains all it's
  162.  * images: The radius must be equal or larger 
  163.  *
  164.  *      |B| / (1-s) 
  165.  *
  166.  * for each IFS map x -> Ax+B with contractivity s.
  167.  * hierarchical IFS: use only the maps with highest rang number */
  168.                 C->x = C->y = C->z = 0.;
  169.                 for (i=0, *R=0.; i < ifs->nrmaps && ifs->maps[i].rang == ifs->maps->rang; i++) 
  170.                         if ((R1 = VecNorm(ifs->maps[i].matrix.translate) /
  171.                                 (1. - ifs->maps[i].contractivity)) > *R) *R = R1;
  172.         } else {
  173. /*
  174.  * construct a sphere around the condensation bounding box 
  175.  */
  176.                 C->x = (ifs->condensation->bounds[HIGH][X] + ifs->condensation->bounds[LOW][X])/2.;
  177.                 C->y = (ifs->condensation->bounds[HIGH][Y] + ifs->condensation->bounds[LOW][Y])/2.;
  178.                 C->z = (ifs->condensation->bounds[HIGH][Z] + ifs->condensation->bounds[LOW][Z])/2.;
  179.         *R = ifs->condensize/2.;
  180.         }
  181.  
  182. /* process maps in groups of same rang */
  183.     for (i = 0; i < ifs->nrmaps; i += n, condensed = TRUE) {
  184. /* count the number of maps with same rang number as ifs->maps[i] */
  185.         for (m = &ifs->maps[i], n = 0; i+n < ifs->nrmaps && m[n].rang == m->rang; n++) {}
  186. /* construct a bounding sphere for the partial IFS defined by these maps of same rang */
  187.         IfsPartialBounding(n, m, condensed, C, R, ifs->options.size, ifs->options.depth);
  188.     }
  189. }
  190.  
  191. /* global variables for depth first processing of the tree of images */
  192. static IfsMap *maps;
  193. static int nrmaps, maxdepth, condensation, oppenheimer_tree;
  194. static Float minsize, isize;
  195. static Float ibounds[2][3], cbounds[2][3], lbounds[2][3];
  196. static Vector ipoint;
  197.  
  198. static IfsBB(map, depth, size, rang)
  199. RSMatrix *map;
  200. int depth;
  201. Float size;
  202. int rang;
  203. {
  204.     int i;
  205.  
  206.     if (depth >= maxdepth || (minsize > EPSILON && size <= minsize)) {
  207.         Vector p;
  208.         Float b[2][3];
  209.         if (oppenheimer_tree) {    /* enlarge with bounds of image of leaves */
  210.             BoundsCopy(lbounds, b);
  211.             BoundsTransform(map, b);
  212.             BoundsEnlarge(ibounds, b);    
  213.         } else if (condensation) {    /* enlarge with bounds of condensation */
  214.             BoundsCopy(cbounds, b);
  215.             BoundsTransform(map, b);
  216.             BoundsEnlarge(ibounds, b);    
  217.         } else {        /* enlarge with sphere with centre at image of ipoint */
  218.             p = ipoint;
  219.             PointTransform(&p, map);
  220.             BoundsEnlargeSphere(ibounds, &p, size/2);
  221.         }
  222.         return;
  223.     }
  224.  
  225.     if (condensation) {        /* enlarge with bounds of image of condensation set */
  226.         Float b[2][3];
  227.         BoundsCopy(cbounds, b);
  228.         BoundsTransform(map, b);
  229.         BoundsEnlarge(ibounds, b);    
  230.     }
  231.  
  232.     for (i=0; i<nrmaps; i++) {
  233.         RSMatrix m;
  234.         Float lambda[3], siz;
  235.  
  236.         if (maps[i].rang < rang) break;        /* for hierarchical IFS */
  237.         MatrixMult(&maps[i].matrix, map, &m);    
  238.             if (maps[i].uniform)            
  239.                 siz = size * maps[i].contractivity;
  240.             else {
  241.                     MatrixScaling(&m, lambda);
  242.                     siz = isize * Lipschitz(lambda);
  243.             }
  244.         IfsBB(&m, depth+1, siz, maps[i].rang);    
  245.     }
  246. }
  247.  
  248. static IfsBoundingBox(ifs, bounds)
  249. Ifs *ifs;
  250. Float bounds[2][3];
  251. {
  252.     RSMatrix unity_transform;
  253.  
  254.     MatrixInit(&unity_transform);
  255.     maps = ifs->maps;
  256.     nrmaps = ifs->nrmaps;
  257.     condensation = (ifs->condensation != NULL);
  258.     oppenheimer_tree = (ifs->leaves != NULL);
  259.     maxdepth = ifs->options.depth;
  260.     minsize = ifs->options.size;
  261.     if (minsize < EPSILON)
  262.         minsize = ifs->extent.r/20;
  263.     isize = ifs->condensation ? ifs->condensize : ifs->size ;
  264.     ipoint = ifs->extent.c;
  265.  
  266.     BoundsInit(ibounds);
  267.     if (ifs->condensation) 
  268.         BoundsCopy(ifs->condensation->bounds, cbounds);
  269.     if (ifs->leaves)
  270.         BoundsCopy(ifs->leaves->bounds, lbounds);
  271.  
  272.     IfsBB(&unity_transform, 0, isize, -99999999);
  273.     BoundsCopy(ibounds, bounds);
  274. }
  275.  
  276. /* compute a bunding box and bounding sphere for the IFS attractor */
  277. static void IfsComputeExtent(ifs)
  278. Ifs *ifs;
  279. {
  280.         Float R, sphavpa, boxavpa;
  281.         Vector C, d;
  282.  
  283. /* compute bounds of condensation set, if any */
  284.     if (ifs->condensation) {
  285.                 GeomComputeBounds(ifs->condensation);
  286.         if (UNBOUNDED(ifs->condensation)) 
  287.                     RLerror(RL_ABORT, "Can't compute bounds of IFS with unbounded condensation.\n");
  288.  
  289.                 d.x = ifs->condensation->bounds[HIGH][X] - ifs->condensation->bounds[LOW][X];
  290.                 d.y = ifs->condensation->bounds[HIGH][Y] - ifs->condensation->bounds[LOW][Y];
  291.                 d.z = ifs->condensation->bounds[HIGH][Z] - ifs->condensation->bounds[LOW][Z];
  292.         ifs->condensize = VecNorm(d);
  293.         }
  294.  
  295. /* compute bounds of leaves, if any */
  296.         if (ifs->leaves) {
  297.                 GeomComputeBounds(ifs->leaves);
  298.                 if (UNBOUNDED(ifs->leaves))
  299.                         RLerror(RL_ABORT, "Oppenheimer tree with unbounded leaves.\n");
  300.  
  301.                 d.x = ifs->leaves->bounds[HIGH][X] - ifs->leaves->bounds[LOW][X];
  302.                 d.y = ifs->leaves->bounds[HIGH][Y] - ifs->leaves->bounds[LOW][Y];
  303.                 d.z = ifs->leaves->bounds[HIGH][Z] - ifs->leaves->bounds[LOW][Z];
  304.         ifs->leavesize = VecNorm(d);
  305.         }
  306.  
  307. /* construct a sphere fitting tightly the attractor of the IFS */
  308.         IfsBoundingSphere(ifs, &R, &C);
  309.  
  310.         ifs->extent.r = R;
  311.         ifs->extent.rsq = R*R;
  312.         ifs->extent.c = C;
  313.         ifs->size = 2*R;
  314.  
  315. /* construct a box fitting tightly the attractor of the IFS */
  316.     IfsBoundingBox(ifs, ifs->bounds);
  317.  
  318. /* compare average projected area of the box and the sphere. The avaerage projected
  319.  * area for a convex volume is one fourth of the surface area of the volume [Glassner,
  320.  * An Introduction to Ray Tracing, Academic Press, 1989 p212]. The volume
  321.  * with the smallest average projected area will be used for building a bounding
  322.  * volume hierarchy for the attractor of the IFS */
  323.     sphavpa = 3.14 * ifs->extent.rsq;
  324.     d.x = ifs->bounds[HIGH][X] - ifs->bounds[LOW][X];
  325.     d.y = ifs->bounds[HIGH][Y] - ifs->bounds[LOW][Y];
  326.     d.z = ifs->bounds[HIGH][Z] - ifs->bounds[LOW][Z];
  327.     boxavpa = (d.x * d.y + d.y * d.z + d.z * d.x) / 2;
  328.     
  329.     if (sphavpa < boxavpa)    ifs->boxbounding = FALSE;
  330.     else            ifs->boxbounding = TRUE;
  331.     if (ifs->options.boxbounding) ifs->boxbounding = TRUE;
  332.     if (ifs->options.spherebounding) ifs->boxbounding = FALSE;
  333. }
  334.  
  335. /*
  336.  * for debugging or listing the volumes that will be rendered
  337.  */
  338. static TransPrint(trans, out)
  339. RSMatrix *trans;
  340. FILE *out;
  341. {
  342.     int i;
  343.  
  344.     for (i=0; i<3; i++) 
  345.         fprintf(out, "%lf %lf %lf   ", 
  346.             (double)trans->matrix[i][0],
  347.             (double)trans->matrix[i][1],
  348.             (double)trans->matrix[i][2]  );
  349.     fprintf(out, "%lf %lf %lf",
  350.             (double)trans->translate.x,
  351.             (double)trans->translate.y,
  352.             (double)trans->translate.z  );
  353. }
  354.  
  355. /* produce a list of volumes that will be rendered */
  356. static FILE *listfp;    
  357. static long volcount;
  358.  
  359. static ListIt(map, depth, size, rang)
  360. RSMatrix *map;
  361. int depth;
  362. Float size;
  363. int rang;
  364. {
  365.     int i;
  366.     RSMatrix m;
  367.     Float lambda[3];
  368.  
  369.     if (depth >= maxdepth || (minsize > EPSILON && size <= minsize)) {
  370.         if (oppenheimer_tree)
  371.             fprintf(listfp, "object leaf         ");
  372.         else if (condensation)
  373.             fprintf(listfp, "object condensation ");
  374.         else
  375.             fprintf(listfp, "object extent ");
  376.         fprintf(listfp, "transform   ");
  377.         TransPrint(map, listfp);
  378.         fprintf(listfp, "\n");
  379.         volcount++;
  380.         return;
  381.     }
  382.  
  383.     if (condensation) {
  384.         fprintf(listfp, "object condensation ");
  385.         fprintf(listfp, "transform   ");
  386.         TransPrint(map, listfp);
  387.         fprintf(listfp, "\n");
  388.         volcount++;
  389.     }
  390.  
  391.     for (i=0; i<nrmaps; i++) {
  392.             Float siz;
  393.         if (maps[i].rang < rang) break;
  394.         MatrixMult(&maps[i].matrix, map, &m);    
  395.             if (maps[i].uniform)            
  396.                 siz = size * maps[i].contractivity;
  397.             else {
  398.                     MatrixScaling(&m, lambda);
  399.                     siz = isize * Lipschitz(lambda);
  400.             }
  401.         ListIt(&m, depth+1, siz, maps[i].rang);
  402.     }
  403. }
  404.  
  405. static IfsList(ifs, filename)
  406. Ifs *ifs;
  407. char *filename;
  408. {
  409.     RSMatrix unity_transform;
  410.  
  411. /* initialise some global variables */
  412.     MatrixInit(&unity_transform);
  413.     maps = ifs->maps;
  414.     nrmaps = ifs->nrmaps;
  415.     condensation = (ifs->condensation != NULL);
  416.     oppenheimer_tree = (ifs->leaves != NULL);
  417.     maxdepth = ifs->options.depth;
  418.     minsize = ifs->options.size;
  419.     isize = ifs->condensation ? ifs->condensize : ifs->size ;
  420.     if (maxdepth > 999999 && minsize < EPSILON) {
  421.         minsize = isize/100;
  422.         fprintf(stderr, "Warning: no stopcriterium for IFS list of volumes. Using 'minsize %lf'\n", (double)minsize);
  423.         return;
  424.     }
  425.  
  426. /* open output file */
  427.     if ((listfp = fopen(filename, "w")) == NULL) {
  428.         fprintf(stderr, "IfsList: Cannot open %s for writing\n", filename);
  429.         perror(ifs->options.listfilename);
  430.         return;
  431.     }
  432.  
  433.     fprintf(listfp, "/* maxdepth=%d, initial_size=%lf, minsize=%lf */\n\n",
  434.              maxdepth, (double)isize, (double)minsize);
  435.  
  436.     if (!ifs->boxbounding) fprintf(listfp, "/* ");
  437.     fprintf(listfp, "name extent box  %lf %lf %lf  %lf %lf %lf", 
  438.         (double)ifs->bounds[LOW][X],
  439.         (double)ifs->bounds[LOW][Y],
  440.         (double)ifs->bounds[LOW][Z],
  441.         (double)ifs->bounds[HIGH][X],
  442.         (double)ifs->bounds[HIGH][Y],
  443.         (double)ifs->bounds[HIGH][Z]  );
  444.     if (!ifs->boxbounding) fprintf(listfp, "*/");
  445.     fprintf(listfp, "\n");
  446.  
  447.     if (ifs->boxbounding) fprintf(listfp, "/* ");
  448.     fprintf(listfp, "name extent sphere  %lf  %lf %lf %lf", 
  449.         (double)ifs->extent.r,
  450.         (double)ifs->extent.c.x, 
  451.         (double)ifs->extent.c.y, 
  452.         (double)ifs->extent.c.z        );
  453.     if (ifs->boxbounding) fprintf(listfp, "*/");
  454.     fprintf(listfp, "\n");
  455.  
  456.     if (condensation)
  457.         fprintf(listfp, "name condensation /* fill in the condensation-set object */\n");
  458.     if (oppenheimer_tree)
  459.         fprintf(listfp, "name leaf /* fill in the leaf object */\n");
  460.  
  461.     volcount = 0;
  462.     fprintf(listfp, "\ngrid 1 1 1 /* change to your personal taste */\n");
  463.     ListIt(&unity_transform, 0, isize, -99999999);
  464.     fprintf(listfp, "end /* %ld objects */\n", volcount);
  465.  
  466.     fclose(listfp);
  467. }
  468.  
  469. /*
  470.  * normal weighting functions 
  471.  */
  472. static Float HighPassNormalWeighting(Size, size)
  473. Float Size, size;
  474. {
  475.         return (Size-size)/Size;
  476. }
  477.  
  478. static Float LowPassNormalWeighting(Size, size)
  479. Float Size, size;
  480. {
  481.         return size/Size;
  482. }
  483.  
  484. static Float ConstantNormalWeighting(Size, size)
  485. Float Size, size;
  486. {
  487.         return 1.;
  488. }
  489.  
  490.  
  491. /*
  492.  * Create & return reference to an ifs.
  493.  */
  494. Ifs *
  495. IfsCreate(maps, condensation, leaves, options)
  496. IfsTransList *maps;
  497. Geom *condensation, *leaves;
  498. IfsOptions *options;
  499. {
  500.         Ifs *ifs;
  501.         IfsTransList *map; 
  502.         Float contractivity, maxcontractivity;
  503.         Float lambda[3];
  504.         int n;
  505.  
  506. /* first: count the number of maps */
  507.         for (n=0, map=maps; map; n++, map=map->next) {}
  508.         if (n<=0) {
  509.                 RLerror(RL_WARN, "I don't create IFS with no maps.\n");
  510.                 return (Ifs *)NULL;
  511.         }
  512.  
  513. /* create Ifs structure */
  514.         ifs = (Ifs *)share_malloc(sizeof(Ifs));
  515.         ifs->condensation = condensation;       /* NULL = no cendensation */
  516.         ifs->leaves = leaves;   /* if not NULL -> Oppenheimer tree */
  517.         ifs->options = *options;
  518.  
  519. /* allocate space for keeping the maps */
  520.         ifs->nrmaps = n;
  521.         ifs->maps = (IfsMap *)share_malloc(n*sizeof(IfsMap));
  522.  
  523. /* the maps are in reverse order */
  524.         ifs->animated = FALSE;
  525.         ifs->uniform_maps = TRUE;
  526.         maxcontractivity = 0.;
  527.         for (n=ifs->nrmaps-1, map=maps; map; n--, map=map->next) {
  528.                 ifs->maps[n].trans = map->trans;
  529.                 ifs->maps[n].transtail = map->transtail;
  530.                 ifs->maps[n].rang = map->rang;
  531.  
  532. /* is the map animated ? */
  533.                 ifs->maps[n].animated = (map->trans->animated || map->trans->next!=(Trans *)NULL);
  534. /* if so, the IFS is animated */
  535.                 if (ifs->maps[n].animated)
  536.                         ifs->animated = TRUE;
  537. /* if not, we can do some computations already */
  538.                 else {
  539. /* "compose" the transformations - they are already composed during the 
  540.  * parsing, so, just copy them. They are also garanteed to be non-singular. */
  541.                         MatrixCopy(&map->trans->trans,  &ifs->maps[n].matrix);
  542.                         MatrixCopy(&map->trans->itrans, &ifs->maps[n].imatrix);
  543.  
  544. /* compute the contractivity of the map */
  545.                         MatrixScaling(&ifs->maps[n].matrix, lambda);
  546.                         contractivity = Lipschitz(lambda);
  547.                         if (contractivity >= 1.) 
  548.                                 RLerror(RL_WARN, "Ifs map %d is not contractive (contractivity = %lf)\n", n+1, contractivity);
  549.                         if (contractivity > maxcontractivity)
  550.                                 maxcontractivity = contractivity;
  551.                         ifs->maps[n].contractivity = contractivity;
  552.                         ifs->maps[n].uniform = Uniform(lambda);
  553.                         if (!ifs->maps[n].uniform) ifs->uniform_maps = FALSE;
  554.                 }
  555.         }
  556.  
  557. /* some procedures may never end for non-contractive IFS's */
  558. /*        if (maxcontractivity >= 1.0) {
  559.                 free(ifs->maps);
  560.                 free(ifs);
  561.                 return (Ifs *)NULL;
  562.         }
  563.  * just warn the user. we'll see ... */
  564.  
  565. /* for hierarchical IFS: sort maps to descending rang number (those maps 
  566.  * that can follow any other map will be processed first) */
  567.     qsort(ifs->maps, ifs->nrmaps, sizeof(IfsMap), rangcmp);
  568.  
  569.         if (!ifs->animated) ifs->contractivity = maxcontractivity;
  570.  
  571. /* the IFS is also animated if its condensation is -- no code to support this yet */
  572. /*        ifs->animated |= ifs->condensation->animtrans; */
  573.  
  574. /* if the IFS is not animated, compute its extent here */
  575.         if (!ifs->animated && !ifs->condensation)
  576.                 IfsComputeExtent(ifs);
  577.  
  578. /* normal weighting function for IFS without condensation */
  579.         if (!ifs->condensation) {
  580.                 switch (ifs->options.normalweighting) {
  581.                         case HIGHPASS_NORMAL_WEIGHTING:
  582.                                 ifs->normalweight = HighPassNormalWeighting; break;
  583.                         case LOWPASS_NORMAL_WEIGHTING:
  584.                                 ifs->normalweight = LowPassNormalWeighting; break;
  585.                         case CONSTANT_NORMAL_WEIGHTING:
  586.                                 ifs->normalweight = ConstantNormalWeighting; break;
  587.                         default: {
  588.                                 RLerror(RL_WARN, "invalid normal weighting method.\n");
  589.                                 free(ifs->maps);
  590.                                 free(ifs);
  591.                                 return (Ifs *)NULL;
  592.                         } 
  593.                 }
  594.         } 
  595.  
  596.         ifs->frame = -1;        /* impossible value */
  597.  
  598.         return ifs;
  599. }
  600.  
  601. /* resolve geometric associates -- untested */
  602. static void IfsResolveAssoc(ifs)
  603. Ifs *ifs;
  604. {
  605.         Trans tr;
  606.         Float contractivity, maxcontractivity;
  607.         Float lambda[3];
  608.         int i;  
  609.  
  610.         maxcontractivity = 0.;
  611.         ifs->uniform_maps = TRUE;
  612.         for (i = 0; i < ifs->nrmaps; i++) {
  613.                 if (ifs->maps[i].animated) {
  614.                         TransResolveAssoc(ifs->maps[i].trans);
  615.                         TransComposeList(ifs->maps[i].trans, &tr);
  616.             ifs->maps[i].matrix = tr.trans;
  617.             ifs->maps[i].imatrix = tr.itrans;
  618.  
  619. /* compute the contractivity of the map */
  620.                         MatrixScaling(&ifs->maps[i].matrix, lambda);
  621.                         contractivity = Lipschitz(lambda);
  622.                         if (contractivity >= 1.) 
  623.                                 RLerror(RL_WARN, "Ifs map %d is not contractive (contractivity = %lf)\n", i+1, contractivity);
  624.                         ifs->maps[i].contractivity = contractivity;
  625.                         ifs->maps[i].uniform = Uniform(lambda);
  626.                 } else 
  627.                         contractivity = ifs->maps[i].contractivity;
  628.  
  629.                 if (contractivity > maxcontractivity)
  630.                         maxcontractivity = ifs->maps[i].contractivity;
  631.                 if (!ifs->maps[i].uniform) ifs->uniform_maps = FALSE;
  632.         }
  633.  
  634.         ifs->contractivity = maxcontractivity;
  635. /*        if (ifs->contractivity >= 1.)
  636.                 RLerror(RL_ABORT, "Can't go on with non-contractive IFS.\n");
  637. */
  638.         ifs->frame = Sampling.framenum;
  639. }
  640.  
  641. Methods *
  642. IfsMethods()
  643. {
  644.     if (iIfsMethods == (Methods *)NULL) {
  645.         iIfsMethods = MethodsCreate();
  646.         iIfsMethods->create = (GeomCreateFunc *)IfsCreate;
  647.         iIfsMethods->methods = IfsMethods;
  648.         iIfsMethods->name = IfsName;
  649.         iIfsMethods->intersect = IfsIntersect;
  650.         iIfsMethods->normal = IfsNormal;
  651.         iIfsMethods->uv = NULL;    /* IfsUV; */
  652.         iIfsMethods->enter = IfsEnter;
  653.         iIfsMethods->bounds = IfsBounds;
  654.         iIfsMethods->stats = IfsStats;
  655.         iIfsMethods->checkbounds = TRUE;
  656.         iIfsMethods->closed = FALSE;
  657.     }
  658.  
  659.     return iIfsMethods;
  660. }
  661.  
  662. Methods *
  663. CIfsMethods()
  664. {
  665.     if (iCIfsMethods == (Methods *)NULL) {
  666.         iCIfsMethods = MethodsCreate();
  667.         iCIfsMethods->create = (GeomCreateFunc *)IfsCreate;
  668.         iCIfsMethods->methods = CIfsMethods;
  669.         iCIfsMethods->name = CIfsName;
  670.         iCIfsMethods->intersect = CIfsIntersect;
  671.         iCIfsMethods->bounds = IfsBounds;
  672.         iCIfsMethods->stats = CIfsStats;
  673.         iCIfsMethods->checkbounds = TRUE;
  674.         iCIfsMethods->closed = TRUE;    /* depends on its condensation */
  675.     }
  676.  
  677.     return iCIfsMethods;
  678. }
  679.  
  680. /* Transformed Extends Lists */
  681. static TEList *TENew()
  682. {
  683.         return (TEList *)Malloc(sizeof(TEList));
  684. }
  685.  
  686. static void TEDispose(p)
  687. TEList *p;
  688. {
  689.         free(p);
  690. }
  691.  
  692. static void TEDisposeList(list)
  693. TEList *list;
  694. {
  695.         TEList *p;
  696.  
  697. /* de-allocate nodes in reverse order of allocation */
  698.         for (p=list; p->next; p=p->next) {}
  699.         list=p;
  700.         while (list) {
  701.                 p = list->prev;
  702.                 TEDispose(list);
  703.                 list = p;
  704.         }
  705. }
  706.  
  707. /* ray - sphere intersection */
  708. static int IfsIntersectSphere(sph, ray, mindist, maxdist)
  709. struct IfsExtent *sph;
  710. Ray *ray;
  711. Float mindist, *maxdist;
  712. {
  713.     Float xadj, yadj, zadj;
  714.     Float b, d, t, s;
  715.  
  716.     xadj = sph->c.x - ray->pos.x;
  717.     yadj = sph->c.y - ray->pos.y;
  718.     zadj = sph->c.z - ray->pos.z;
  719.         d = xadj * xadj + yadj * yadj + zadj * zadj - sph->rsq;
  720.  
  721.     b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
  722.     t = b * b - d;
  723.     if (t < 0.)
  724.         return FALSE;
  725.     t = (Float)sqrt((double)t);
  726.     s = b - t;
  727.     if (s > mindist) {
  728.         if (s < *maxdist) {
  729.             *maxdist = s;
  730.             return TRUE;
  731.         }
  732.         return FALSE;
  733.     }
  734.     s = b + t;
  735.     if (s > mindist && s < *maxdist) {
  736.         *maxdist = s;
  737.         return TRUE;
  738.     }
  739.     return FALSE;
  740. }
  741.  
  742. static int IfsExtentIntersect(ifs, ray, mindist, maxdist)
  743. Ifs *ifs;
  744. Ray *ray;
  745. Float mindist, *maxdist;
  746. {
  747.         Vector vtmp;
  748.  
  749.         VecAddScaled(ray->pos, mindist, ray->dir, &vtmp);
  750.         if (ifs->boxbounding) {       /*  test against solid box */
  751.                 if (OutOfBounds(&vtmp, ifs->bounds))
  752.                         return BoundsIntersect(ray, ifs->bounds, mindist, maxdist);
  753.         } else {                       /* test against solid sphere */
  754.                 VecSub(vtmp, ifs->extent.c, &vtmp);
  755.                 if (vtmp.x*vtmp.x + vtmp.y*vtmp.y + vtmp.z*vtmp.z > ifs->extent.rsq)
  756.                         return IfsIntersectSphere(&ifs->extent, ray, mindist, maxdist);
  757.         } 
  758.         *maxdist = mindist; 
  759.         return TRUE;
  760. }
  761.  
  762. /* for IFS without condensation */
  763. static void IfsExtentNormal(ifs, ray, dist, normal)
  764. Ifs *ifs;
  765. Ray *ray;
  766. Float dist;
  767. Vector *normal;    
  768. {
  769.     Vector pos;
  770.  
  771.         pos.x = ray->pos.x + dist * ray->dir.x;
  772.     pos.y = ray->pos.y + dist * ray->dir.y;
  773.         pos.z = ray->pos.z + dist * ray->dir.z;
  774.  
  775.     if (ifs->boxbounding) {     /* this produces nicer normals, so ... */
  776.         normal->x = pos.x - (ifs->bounds[HIGH][X] + ifs->bounds[LOW][X])/2.;
  777.         normal->y = pos.y - (ifs->bounds[HIGH][Y] + ifs->bounds[LOW][Y])/2.;
  778.         normal->z = pos.z - (ifs->bounds[HIGH][Z] + ifs->bounds[LOW][Z])/2.;
  779.         VecNormalize(normal); 
  780.     } else {
  781.             normal->x = (pos.x - ifs->extent.c.x) / ifs->extent.r;
  782.             normal->y = (pos.y - ifs->extent.c.y) / ifs->extent.r;
  783.             normal->z = (pos.z - ifs->extent.c.z) / ifs->extent.r;
  784.         }
  785. }
  786.  
  787. static int DoIfsIntersect();
  788.  
  789. /* an IFS without condensation set is a primitive object, so no hitlist here */
  790. int 
  791. IfsIntersect(ifs, ray, mindist, maxdist)
  792. Ifs *ifs;
  793. Ray *ray;
  794. Float mindist, *maxdist;
  795. {
  796.         IfsTests++;
  797.         if (DoIfsIntersect(ifs, ray, (HitList *)NULL, mindist, maxdist)) {
  798.                 IfsHits++;
  799.                 return TRUE;
  800.         }
  801.         return FALSE;
  802. }
  803.  
  804. /* an IFS with condensation is like an aggregate */
  805. int 
  806. CIfsIntersect(ifs, ray, hitlist, mindist, maxdist)
  807. Ifs *ifs;
  808. Ray *ray;
  809. HitList *hitlist;
  810. Float mindist, *maxdist;
  811. {
  812.         CIfsTests++;
  813.         if (DoIfsIntersect(ifs, ray, hitlist, mindist, maxdist)) {
  814.                 CIfsHits++;
  815.                 return TRUE;
  816.         }
  817.         return FALSE;
  818. }
  819.  
  820. /*
  821.  * Ray/ifs intersection test.
  822.  */
  823. static int DoIfsIntersect(ifs, ray, hitlist, mindist, maxdist)
  824. Ifs *ifs;
  825. Ray *ray;
  826. HitList *hitlist;
  827. Float mindist, *maxdist;
  828. {
  829.         TEList *list, *closest, *p, *q, *new;
  830.         Float nmindist, nmaxdist, distfac, pixelsize;
  831.         Ray newray;
  832.         Vector normal;
  833.         int n, intersection_found, construct_children;
  834.         Float lambda[3];
  835.         Geom *obj;
  836.     HitNode *np;
  837.     IfsMap *map;
  838.  
  839. /*        if (ifs->animated && !equal(ifs->time, ray->time)) {
  840.                 IfsResolveAssoc(ifs);
  841.                 IfsComputeExtent(ifs);
  842.                 ifs->time = ray->time;
  843.         }
  844. */
  845. /* does the ray intersect the IFS's extent ? */
  846.         distfac = 1.;
  847.         newray = *ray;
  848.         nmaxdist = *maxdist * distfac;
  849.         if (!IfsExtentIntersect(ifs, &newray, EPSILON, &nmaxdist))
  850.                 return FALSE;
  851.  
  852. /* put the extent on the list */
  853.         new = TENew();
  854.         MatrixInit(&new->matrix);
  855.         MatrixInit(&new->imatrix);
  856.         new->ray = newray;
  857.         new->distfac = distfac;
  858.         new->dist = nmaxdist / distfac;
  859.         new->size = ifs->size;
  860.         new->depth = 0;
  861.         new->rang = -99999999;
  862.         if (!ifs->condensation) {
  863.                 IfsExtentNormal(ifs, &newray, nmaxdist, &normal);
  864.                 NormalTransform(&normal, &new->imatrix);
  865.                 VecScale((*ifs->normalweight)(ifs->size, new->size), normal, &new->normal);
  866.         }
  867.         new->next = new->prev = (TEList *)NULL;
  868.         list = new;
  869.         intersection_found = FALSE;
  870.  
  871. /* while list not empty */
  872.         while (list) {
  873. /* find closest entry in list */
  874.                 closest = list;
  875.                 for (p = list->next; p; p = p->next) 
  876.                         if (p->dist < closest->dist) closest = p;
  877.  
  878.                 construct_children = TRUE;
  879.                 pixelsize = PixelSize(ray, closest->dist);
  880.  
  881. /* does the (inverse transformed) ray intersect the condensation set (if any) ? */
  882.                 if (ifs->condensation) {
  883. /* Oppenheimer trees: if TE small enough, than test against leave, not against condensation */
  884.                         obj = ifs->condensation;
  885.             if (ifs->condensize * closest->size / ifs->size < ifs->options.size ||
  886.                 closest->depth >= ifs->options.depth) {
  887.                 if (ifs->leaves) obj = ifs->leaves;
  888.                 construct_children = FALSE;
  889.             }
  890.  
  891. /* instantiate the object */
  892.                         nmindist = mindist * closest->distfac;
  893.                         nmaxdist = *maxdist * closest->distfac;
  894.                         if (intersect(obj, &closest->ray, hitlist, nmindist, &nmaxdist)) {
  895.                                 *maxdist = nmaxdist / closest->distfac;
  896.                                 intersection_found = TRUE;
  897.  
  898. if (*maxdist <= closest->dist) fprintf(stderr, "impossible situation: closest = %p, closest->dist = %lf\n", closest, closest->dist);
  899. /*                if (closest->depth > 0) 
  900.                       MatrixInvert(&closest->matrix, &closest->imatrix); 
  901. */
  902.                 np = &hitlist->data[hitlist->nodes-1];
  903.                 if (np->dotrans) {
  904.                     MatrixMult(&closest->matrix, &np->trans.trans, &np->trans.trans);
  905.                     MatrixMult(&np->trans.itrans, &closest->imatrix, &np->trans.itrans);
  906.                 } else {
  907.                                 MatrixCopy(&closest->matrix, &np->trans.trans);
  908.                                 MatrixCopy(&closest->imatrix, &np->trans.itrans);
  909.                 }
  910.                 np->dotrans = TRUE;
  911.  
  912. /* unless we don't need the *closest* intersection (shadow-rays!!) 
  913.  * we cannot stop here because the generated hierarchy of extents is
  914.  * *not* well-behaved: it is possible that we find a closer intersection
  915.  * later.  */
  916.                                 if (ray->type == SHADOW_RAY) {
  917.                                         TEDisposeList(list);
  918.                                         return TRUE;
  919.                                 }
  920.  
  921. /* But we can remove all the transformed extents which are behind this one. Note
  922.  * that at least one TE (closest) stays in the list. */
  923.                                  for (p = list; p;) {
  924.                                         if (p->dist > *maxdist + EPSILON) {
  925.                                                 if (p->next) p->next->prev = p->prev;
  926.                                                 if (p->prev) p->prev->next = p->next;
  927.                                                 else list = p->next;      
  928.                                                 q = p->next;
  929.                                                 TEDispose(p);
  930.                                                 p = q;
  931.                                         } else
  932.                                                 p = p->next;
  933.                                 }
  934.                         }
  935.                 }
  936.  
  937.  
  938. /* or is the (transformed) extent smaller than one pixel ? */
  939.                 if (closest->size < pixelsize || closest->size < ifs->options.size || closest->depth >= ifs->options.depth) {
  940.                         if (!ifs->condensation &&
  941.                             closest->dist >= mindist && closest->dist <= *maxdist) {
  942.                                 intersection_found = TRUE;
  943.                                 *maxdist = closest->dist;
  944. /* remember hierarchically computed normal and intersection position */
  945.                                 ifs->normal = closest->normal;
  946.                                 VecNormalize(&ifs->normal);
  947.                                 VecAddScaled(ray->pos, *maxdist, ray->dir, &ifs->intersectpos);
  948. /* and we can stop here */
  949.                                 TEDisposeList(list);
  950.                                 return TRUE;
  951.                         }
  952.  
  953. /* anyway, we don't construct the children of the transformed extent */
  954.                         construct_children = FALSE;
  955.                 }
  956.  
  957.  
  958. /* construct the (transformed) extent's "children" and add any that intersect to the list */
  959.                if (construct_children) for (n=0, map=ifs->maps; n<ifs->nrmaps; n++, map++) {
  960. /* for hierarchical IFS: use only maps of equal or higher rang. The maps are sorted to
  961.  * descending rang number */
  962.                         if (map->rang < closest->rang) break;
  963.                         distfac = closest->distfac;
  964.                         newray = closest->ray;
  965.                         distfac *= RayTransform(&newray, &map->imatrix);
  966.                         nmaxdist = *maxdist * distfac;
  967.  
  968.                         if (IfsExtentIntersect(ifs, &newray, EPSILON*distfac, &nmaxdist)) {
  969. /* add the "transformed extent" to the list */
  970.                                 new = TENew();
  971.                                 new->ray = newray;
  972.                                 new->distfac = distfac;
  973.                                 new->dist = nmaxdist / distfac;
  974.                                 new->depth = closest->depth+1;
  975.                                 new->rang = map->rang;
  976. /* pre-multiplication !! */
  977.                 if (ifs->condensation || !ifs->uniform_maps)
  978.                                     MatrixMult(&map->matrix, &closest->matrix, &new->matrix);
  979.                                 MatrixMult(&closest->imatrix, &map->imatrix, &new->imatrix);
  980.  
  981. /* compute size of transformed extent */
  982.                                 if (map->uniform)
  983.                                         new->size = closest->size * map->contractivity;
  984.                                 else {
  985.                                         MatrixScaling(&new->matrix, lambda);
  986.                                         new->size = ifs->size * Lipschitz(lambda);
  987.                                 }
  988.  
  989. /* if IFS without condensation, compute hierarchically composed normal */
  990.                                 if (!ifs->condensation) {
  991.                                         IfsExtentNormal(ifs, &newray, nmaxdist, &normal);
  992.                                         NormalTransform(&normal, &new->imatrix);
  993.                                         VecAddScaled(closest->normal,
  994.                                                      (*ifs->normalweight)(ifs->size, new->size),
  995.                                                      normal,
  996.                                                      &new->normal);
  997.                                 }
  998. /* add to list */
  999.  
  1000.                                 new->prev = list->prev; list->prev = new;
  1001.                                 new->next = list; list = new;
  1002.                         }
  1003.                 }
  1004.  
  1005. /* remove closest from list */
  1006.                 if (closest->next) closest->next->prev = closest->prev;
  1007.                 if (closest->prev) closest->prev->next = closest->next;
  1008.                 else list = closest->next;      /* closest was the first item */
  1009.                 TEDispose(closest);
  1010.         }
  1011.  
  1012.         return intersection_found;
  1013. }
  1014.  
  1015. #define VecEqual(a, b)  ( (equal((a).x, (b).x)) && (equal((a).y, (b).y)) && (equal((a).z, (b).z)) )
  1016.  
  1017. /*
  1018.  * Compute normal to ifs at position of last intersection found.
  1019.  */
  1020. int
  1021. IfsNormal(ifs, pos, nrm, gnrm)
  1022. Ifs *ifs;
  1023. Vector *pos, *nrm, *gnrm;
  1024. {
  1025. /* this situation does not happen in rayshade.4.0.6 */
  1026.         if (!VecEqual(*pos, ifs->intersectpos))
  1027.                 RLerror(RL_ABORT, "Can't compute ifs normal in other point than last intersection.\n");
  1028.  
  1029. /* normal is the hierarchically computed normal */
  1030.         if (ifs->normal.x == 0. && ifs->normal.y == 0. && ifs->normal.z == 0.)
  1031.                 ifs->normal.z = 1.;     /* may happen when the whole IFS is smaller than one pixel */
  1032.         *nrm = *gnrm = ifs->normal;
  1033.         return FALSE;
  1034. }
  1035.  
  1036.  
  1037. /*
  1038.  * Determine if ray enters (TRUE) or leaves (FALSE) ifs at position of
  1039.  * last intersection.
  1040.  */
  1041. int
  1042. IfsEnter(ifs, ray, mind, hitd)
  1043. Ifs *ifs;
  1044. Ray *ray;
  1045. Float mind, hitd;
  1046. {
  1047.         return TRUE;
  1048. }
  1049.  
  1050. void
  1051. IfsBounds(ifs, bounds)
  1052. Ifs *ifs;
  1053. Float bounds[2][3];
  1054. {
  1055.         if (ifs->animated) 
  1056.                 IfsResolveAssoc(ifs);
  1057.  
  1058.         if (ifs->animated || ifs->condensation)
  1059.                 IfsComputeExtent(ifs);
  1060.  
  1061.         BoundsCopy(ifs->bounds, bounds);
  1062.  
  1063. /* produce a list of volumes that will be rendered */
  1064.     if (ifs->options.listfilename) 
  1065.         IfsList(ifs, ifs->options.listfilename);
  1066. }
  1067.  
  1068. char *
  1069. IfsName()
  1070. {
  1071.     return ifsName;
  1072. }
  1073.  
  1074. char *
  1075. CIfsName()
  1076. {
  1077.     return CifsName;
  1078. }
  1079.  
  1080. void
  1081. IfsStats(tests, hits)
  1082. unsigned long *tests, *hits;
  1083. {
  1084.     *tests = IfsTests;
  1085.     *hits = IfsHits;
  1086. }
  1087.  
  1088. void
  1089. CIfsStats(tests, hits)
  1090. unsigned long *tests, *hits;
  1091. {
  1092.     *tests = CIfsTests;
  1093.     *hits = CIfsHits;
  1094. }
  1095.  
  1096. void
  1097. IfsMethodRegister(meth)
  1098. UserMethodType meth;
  1099. {
  1100.     if (iIfsMethods)
  1101.         iIfsMethods->user = meth;
  1102. }
  1103.  
  1104.