home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 07 - 1991 / 07.05 May 91 / Fractal Code / MountainsAlgorithm.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-14  |  13.1 KB  |  507 lines  |  [TEXT/KAHL]
  1. /*
  2.     MountainsAlgorithm.c
  3.     A fractal mountain generating algorithm
  4.     By Ben Haller 1/1/90
  5.     Technical Assistance by Eli Meir
  6.     ©1990 Ben Haller
  7. */
  8.  
  9. #include <math.h>
  10. #include "Mountains.h"
  11.  
  12. /* Imported declarations from MountainsToolbox.c    */
  13. extern CursHandle watch;
  14. extern int wx,wy,wd;
  15. extern char colorF;
  16. extern char menuCk[5];
  17.  
  18. /* Algorithmic declarations */
  19. #define maxIter    9L
  20.  
  21. long bwPat[17][2]={        /* Black and white patterns for shading                */
  22.     {0xFFFFFFFF, 0xFFFFFFFF},
  23.     {0xFFFBFFBF, 0xFFFBFFBF},
  24.     {0xFFEEFFBB, 0xFFEEFFBB},
  25.     {0xFFEEFFAA, 0xFFEEFFAA},
  26.     {0xBBEEBBEE, 0xBBEEBBEE},
  27.     {0xAA77AAFF, 0xAA77AAFF},
  28.     {0xAA77AADD, 0xAA77AADD},
  29.     {0xAA77AA55, 0xAA77AA55},
  30.     {0xAA55AA55, 0xAA55AA55},
  31.     {0xAA55AA11, 0xAA55AA11},
  32.     {0xAA44AA11, 0xAA44AA11},
  33.     {0xAA44AA00, 0xAA44AA00},
  34.     {0x22882288, 0x22882288},
  35.     {0x2200AA00, 0x2200AA00},
  36.     {0x22008800, 0x22008800},
  37.     {0x80000800, 0x80000800},
  38.     {0x00000000, 0x00000000}
  39. };
  40.  
  41. typedef struct {        /* 3D point - double FP precision                    */
  42.     double x,y,z;
  43. } point3;
  44.  
  45. double xTh=-1.396;        /* Angle of tilt (in radians)                        */
  46. double tm[3][3];        /* Transformation matrix                            */
  47. long *map;                /* A pointer to the depth values of the map            */
  48.  
  49. /* Variables dependant on menu selections */
  50. long nIter;                /* Number of iterations of divide-and-conquer        */
  51. long mapS;                /* Length of one side of the map                    */
  52. long mapD;                /* Map dimension - equals mapS-1                    */
  53.  
  54. double contour;            /* Power to raise heights to, to make contour        */
  55. double roughness;        /* Base to raise to get scale variance                */
  56. long normalHeight;        /* Normalized height                                */
  57. char profile;            /* Profile number                                    */
  58.  
  59. /* Drawing variables - precalculated coefficients */
  60. double xc,sc;
  61.  
  62. /* The Map #define lets us access map easily, although
  63.     it is a variable-dimension array.                        */
  64. #define Map(x,y)    map[(x)+mapS*(y)]
  65.  
  66. /********************************************************************/
  67. /*                                                                    */
  68. /*        Calculation Routines                                        */
  69. /*                                                                    */
  70. /********************************************************************/
  71.  
  72. /* Make the mountains a standard height, and simultaneously apply the contour transform */
  73. void NormalizeMap()
  74. {
  75.     register long n,m,k;
  76.     register long *a;
  77.     register double i,z;
  78.     
  79.     a=map;
  80.     m=0x80000000;
  81.     for (n=mapS*mapS-1; n>=0; n--)            /* Find highest point in map                */
  82.         if (*(a+n)>m)
  83.             m=*(a+n);
  84.     z=(double)m;
  85.     z=pow(z,contour);                        /* Apply contour transform to it            */
  86.     z=(double)normalHeight/z;                /* Find normalizing coefficient such that    */
  87.                                             /* the range will be the chosen height        */
  88.     for (n=mapS*mapS-1; n>=0; n--) {        /* For each point in the map                */
  89.         k=*(a+n);
  90.         if (k>0) {                            /* Special case for positive & neg. heights    */
  91.             i=(double)k;
  92.             i=pow(i,contour)*z;                /* Apply contour transform, normalize        */
  93.             *(a+n)=(long)i;                    /* Poke back in map                            */
  94.         } else {
  95.             i=(double)(-k);
  96.             i=pow(i,contour)*z;                /* Apply contour transform, normalize        */
  97.             *(a+n)=-((long)i);                /* Poke back in map, preserving sign        */
  98.         }
  99.     }
  100. }
  101.  
  102. /* Returns a random long integer.  This routine just generates two integers using    */
  103. /* QuickDraw's random number generator and tacks them together.  I have no idea how    */
  104. /* good the numbers thus generated are, but I suspect not very good at all.  But,    */
  105. /* they're good enough.                                                                */
  106. long Rand(first, last)
  107. long first, last;
  108. {
  109.     long r;
  110.     
  111.     do {
  112.         r=(((long)Random())<<16)+Random();
  113.     } while (r<0);
  114.     return(r % (last - first + 1) + first);
  115. }
  116.  
  117. /* Returns the maximum deviation a point could attain given an iteration depth.    */
  118. /* The exact number returned depends on roughness, but this function is always    */
  119. /* strictly decreasing monotonic for given depth ic.                            */
  120. long MaxDeviation(ic)
  121. long ic;
  122. {
  123.     if (roughness>0)
  124.         return ((long)(pow(roughness,(double)(ic-1L))*8.0));
  125.     else return 100000L;
  126. }
  127.  
  128. /* This routine deviates a point by a random amount from -m to m, m being the    */
  129. /* maximum deviation for the given iteration.  If roughness==0 (symbolic of        */
  130. /* infinite smoothness), it is not deviated at all.                                */
  131. void DeviatePoint(long,long);
  132. void DeviatePoint(o,ic)
  133. long o,ic;
  134. {
  135.     long v;
  136.     
  137.     if (roughness<0) return;
  138.     v=MaxDeviation(ic);
  139.     map[o]+=Rand(-v,v);
  140. }
  141.  
  142. /* Passed a triangle with "side" (i.e. horizontal side) (sx1,sy1)-(sx2,sy2)    */
  143. /* and "apex" (ax,ay).  It calculates midpoints and recurses.  Parameter    */
  144. /* 'c' gives a standard iteration count                                        */
  145. void IterCalc(long,long,long,long);
  146. void IterCalc(s1,s2,a,c)
  147. long s1,s2,a,c;
  148. {
  149.     register long ns1,ns2,na;
  150.     
  151.     /* Decrement iteration count */
  152.     c--;
  153.     
  154.     /* Find midpoints */
  155.     ns1=(s1+a)>>1;
  156.     ns2=(s2+a)>>1;
  157.     na=(s1+s2)>>1;
  158.     
  159.     /* For each midpoint, if it hasn't already been set, set it to the    */
  160.     /* average of it's endpoints, and deviate by a random amount        */
  161.     if (map[ns1]==0x80000000L) {
  162.         map[ns1]=(map[s1]+map[a])>>1;
  163.         DeviatePoint(ns1,c);
  164.     }
  165.     if (map[ns2]==0x80000000L) {
  166.         map[ns2]=(map[s2]+map[a])>>1;
  167.         DeviatePoint(ns2,c);
  168.     }
  169.     if (map[na]==0x80000000L) {
  170.         map[na]=(map[s1]+map[s2])>>1;
  171.         DeviatePoint(na,c);
  172.     }
  173.     
  174.     /* Iterate calculations on sub-triangles if we haven't yet reached    */
  175.     /* the maximum resolution of the map                                */
  176.     if (ns1+1!=ns2) {
  177.         IterCalc(s1,na,ns1,c);
  178.         IterCalc(na,s2,ns2,c);
  179.         IterCalc(ns1,ns2,na,c);
  180.         IterCalc(ns1,ns2,a,c);
  181.     }
  182. }
  183.  
  184. /* Initialize the entire map to 0x80000000 (<not set> value) */
  185. void InitMap()
  186. {
  187.     register long n;
  188.     register long *a;
  189.     
  190.     a=map;
  191.     for (n=mapS*mapS-1; n>=0; n--)
  192.         *(a+n)=0x80000000L;
  193. }
  194.  
  195. /* Initialize values from menu choices, allocate map, etc. */
  196. void InitMountains()
  197. {
  198.     nIter=menuCk[0]+2;
  199.     map=0L;
  200.     mapD=1L<<nIter;
  201.     mapS=mapD+1;
  202.     
  203.     contour=0.25*menuCk[1];
  204.     if (menuCk[1]==9) contour=3.0;
  205.     else if (menuCk[1]==10) contour=5.0;
  206.     
  207.     roughness=0.75+0.25*menuCk[2];
  208.     if (menuCk[2]==7) roughness=2.75;
  209.     else if (menuCk[2]==8) roughness=3.5;
  210.     else if (menuCk[2]==9) roughness=5.0;
  211.     else if (menuCk[2]==10) roughness=-1.0;
  212.     
  213.     normalHeight=10000L*menuCk[3];
  214.     
  215.     profile=menuCk[4];
  216.     
  217.     if (map==0L) map=NewPtr(mapS*mapS*4L);
  218.     if (map==0L) ExitToShell();    /* Out of memory */
  219.     InitMap();
  220. }
  221.  
  222. /* Calculate mountain range using current parameters, profile, etc. */
  223. void CalcMountains()
  224. {
  225.     register long q;
  226.     
  227.     SetCursor(*watch);
  228.     InitMountains();
  229.     
  230.     /* Generate starting profile to build on */
  231.     q=MaxDeviation(maxIter+1)/2;
  232.     switch (profile) {
  233.         case 1:
  234.             Map(0,0)=q;
  235.             Map(mapD,0)=0;
  236.             Map(0,mapD)=0;
  237.             Map(mapD,mapD)=-q;
  238.             break;
  239.         case 2:
  240.             Map(0,0)=q;
  241.             Map(mapD,0)=0;
  242.             Map(0,mapD)=0;
  243.             Map(mapD,mapD)=0;
  244.             Map(mapD>>1,mapD>>1)=0;
  245.             Map(mapD>>1,mapD)=0;
  246.             Map(mapD,mapD>>1)=0;
  247.             break;
  248.         case 3:
  249.             Map(0,0)=0;
  250.             Map(mapD,0)=0;
  251.             Map(0,mapD)=0;
  252.             Map(mapD,mapD)=-q;
  253.             break;
  254.         case 4:
  255.             Map(0,0)=0;
  256.             Map(mapD,0)=0;
  257.             Map(0,mapD)=0;
  258.             Map(mapD,mapD)=0;
  259.             Map(mapD>>1,mapD>>1)=q/2;
  260.             Map(mapD>>1,0)=q;
  261.             Map(0,mapD>>1)=q;
  262.             break;
  263.     }
  264.     
  265.     /* Generate each main triangle recursively */
  266.     IterCalc(0L,mapS-1,mapS*mapS-1,maxIter+1);
  267.     IterCalc(mapS*(mapS-1),mapS*mapS-1,0L,maxIter+1);
  268.     
  269.     /* Normalize to standard height, apply contour transform */
  270.     NormalizeMap();
  271.     
  272.     SetCursor(&arrow);
  273. }
  274.  
  275. /********************************************************************/
  276. /*                                                                    */
  277. /*        Drawing Routines                                            */
  278. /*                                                                    */
  279. /********************************************************************/
  280.  
  281. /* Transform from map coordinates to screen coordinates*/
  282. void CalcPoint3(p)
  283. point3 *p;
  284. {
  285.     register double xp,yp,zp;
  286.     
  287.     xp=(p->x*2-p->y+mapD)*xc;
  288.     yp=(p->y*2)*xc;
  289.     zp=p->z*0.00217;
  290.     
  291.     p->x=xp*sc;
  292.     p->y=(yp*tm[1][1]+zp*tm[2][1]+230)*sc;
  293.     p->z=(yp*tm[1][2]+zp*tm[2][2])*sc;
  294. }
  295.  
  296. /* This routine actually draws a triangle, given by a set of three (x,y,z) triplets.    */
  297. /* It determines the color and shading according to altitude and lighting, constructs    */
  298. /* a QuickDraw polygon for the triangle, and draws it.                                    */
  299. void _DrawTriangle(p)
  300. point3 *p;
  301. {
  302.     double nx,ny,nz,lx,ly,lz,i;
  303.     RGBColor c;
  304.     PolyHandle ph;
  305.     long color,slope;
  306.     
  307.     /* Figure out what base color our triangle will be : blue, green, or gray */
  308.     if ((p[0].z==0.0) && (p[1].z==0.0) && (p[2].z==0.0))
  309.         color=-1;                            /* blue        */
  310.     else {
  311.         double az,treeline;
  312.         
  313.         treeline=.4*normalHeight+10000;                /* Calculate treeline    */
  314.         az=(p[0].z+p[1].z+p[2].z)/3.0;                /* Find avg. height        */
  315.         if (az>treeline+9000) color=150;            /* gray                    */
  316.         else if (az<treeline) color=0;                /* green                */
  317.         else color=(az-treeline)/60;                /* intermediate         */
  318.     }
  319.     
  320.     /* Transform into viewing space */
  321.     CalcPoint3(p+0);
  322.     CalcPoint3(p+1);
  323.     CalcPoint3(p+2);
  324.     
  325.     /* Create a Quickdraw polygon to be the triangle in question */
  326.     ph=OpenPoly();
  327.     MoveTo((int)p[0].x,(int)p[0].y);
  328.     LineTo((int)p[1].x,(int)p[1].y);
  329.     LineTo((int)p[2].x,(int)p[2].y);
  330.     LineTo((int)p[0].x,(int)p[0].y);
  331.     ClosePoly();
  332.     
  333.     /* Calculate the normal vector to our triangle, by means of a cross product */
  334.     {
  335.         double v1x,v1y,v1z,v2x,v2y,v2z;
  336.     
  337.         v1x=p[1].x-p[0].x;
  338.         v1y=p[1].y-p[0].y;
  339.         v1z=p[1].z-p[0].z;
  340.         
  341.         v2x=p[2].x-p[0].x;
  342.         v2y=p[2].y-p[0].y;
  343.         v2z=p[2].z-p[0].z;
  344.         
  345.         nx=v1y*v2z-v1z*v2y;
  346.         ny=v1z*v2x-v1x*v2z;
  347.         nz=v1y*v2x-v1x*v2y;
  348.     }
  349.     
  350.     /* Initialize our light source */
  351.     lx=-1.0/sqrt(3.0);
  352.     ly=-lx; lz=-lx;
  353.     
  354.     /* Figure out what color we want */
  355.     if (color==-1) {
  356.         c.red=0x9FFF; c.green=0x9FFF; c.blue=0xFFFF;        /* Water */
  357.         if (colorF) RGBForeColor(&c);
  358.         else PenPat(bwPat[8]);
  359.     } else {
  360.         unsigned long cv[3];
  361.         
  362.         /* Start with green */
  363.         cv[0]=0x00005555L; cv[1]=0x0000FFFFL; cv[2]=0x0000AAAAL;
  364.         
  365.         /* Blend to brown and white over treeline */
  366.         cv[0]=((150L-color)*(cv[0]-0x00002000L))/150L+0x00002000L;
  367.         cv[1]=(cv[1]*(150L-color))/150L;
  368.         cv[2]=0x0000EFFFL-((150L-color)*(0x0000EFFFL-cv[2]))/150L;
  369.         
  370.         /* Move into an RGBColor structure */
  371.         c.red=cv[0]; c.green=cv[1]; c.blue=cv[2];
  372.         
  373.         /* Scale brightness according to the incidence of light */
  374.         i=((lx*nx+ly*ny+lz*nz)/sqrt(nx*nx+ny*ny+nz*nz))*0.5+0.5;
  375.         c.blue*=i;
  376.         
  377.         if (colorF) {        /* Pick our color */
  378.             HSV2RGB(&c,&c);
  379.             RGBForeColor(&c);
  380.         } else {            /* Pick a pattern based on brightness */
  381.             int pat;
  382.             
  383.             pat=(((unsigned long)c.blue)*36L)/0x00010000L;
  384.             pat-=10;
  385.             if (pat<0) pat=0;
  386.             if (color==150) {
  387.                 if (pat>16) pat=16;
  388.             } else if (pat>15) pat=15;
  389.             PenPat(bwPat[pat]);
  390.         }
  391.     }
  392.     
  393.     PaintPoly(ph);            /* Draw our triangle    */
  394.     KillPoly(ph);            /* Get rid of it        */
  395. }
  396.  
  397. /* Draw a given triangle.  This routine is mainly concerned with the possibility that    */
  398. /* a triangle could span the waterline.  If this occurs, this procedure breaks it up    */
  399. /* into three smaller triangles, each of which is either above or below water.  All        */
  400. /* actual drawing or coloration is delegated to _DrawTriangle, above.                    */
  401. void DrawTriangle(x1,y1,x2,y2,x3,y3)
  402. long x1,y1,x2,y2,x3,y3;
  403. {
  404.     long z[3];
  405.     point3 p[3];
  406.     
  407.     /* Get depths of the three points */
  408.     z[0]=Map(x1,y1);
  409.     z[1]=Map(x2,y2);
  410.     z[2]=Map(x3,y3);
  411.     if ((z[0]<=0) && (z[1]<=0) && (z[2]<=0)) {        /* All underwater */
  412.         p[0].x=x1; p[0].y=y1; p[0].z=0;
  413.         p[1].x=x2; p[1].y=y2; p[1].z=0;
  414.         p[2].x=x3; p[2].y=y3; p[2].z=0;
  415.         _DrawTriangle(p);
  416.     } else if ((z[0]<0) || (z[1]<0) || (z[2]<0)) {    /* Some underwater - split at waterline */
  417.         point3 ap,s[2],m[2];
  418.         char w[3];
  419.         int n;
  420.         double f;
  421.         
  422.         p[0].x=x1; p[0].y=y1; p[0].z=z[0];
  423.         p[1].x=x2; p[1].y=y2; p[1].z=z[1];
  424.         p[2].x=x3; p[2].y=y3; p[2].z=z[2];
  425.         for (n=0; n<3; n++)
  426.             w[n]=(z[n]<0);
  427.         if ((w[0]!=w[1]) && (w[0]!=w[2])) {
  428.             ap=p[0]; s[0]=p[1]; s[1]=p[2];
  429.         } else {
  430.             if (w[1]!=w[0]) {
  431.                 s[1]=p[0]; ap=p[1]; s[0]=p[2];
  432.             } else {
  433.                 s[0]=p[0]; s[1]=p[1]; ap=p[2];
  434.             }
  435.         }
  436.         
  437.         /* At this point, ap is the "odd man out" - either it is above water and the other    */
  438.         /* two are below, or it is below and the other two are above.  Which corner s[0]    */
  439.         /* is and which s[1] is *is* important - if we get the wrong order, the normal        */
  440.         /* vector used to find the shading coefficient is the wrong sign.  This is true        */
  441.         /* whenever we are manipulating corners - the ordering is always important.            */
  442.         
  443.         /* Find the "midpoints" between ap and s[0]&s[1] - this is where we split our big    */
  444.         /* triangle into smaller triangles.  Actually it is not a normal midpoint, but a    */
  445.         /* weighted midpoint, such that the z component is 0 - waterline.                    */
  446.         for (n=0; n<2; n++) {
  447.             f=-((double)ap.z)/(double)(s[n].z-ap.z);
  448.             
  449.             m[n].x=ap.x-(ap.x-s[n].x)*f;
  450.             m[n].y=ap.y-(ap.y-s[n].y)*f;
  451.             m[n].z=0;
  452.         }
  453.         
  454.         /* Set whichever triangles are below water to 0 altitude */
  455.         if (ap.z<0) ap.z=0;
  456.         else {
  457.             s[0].z=0; s[1].z=0;
  458.         }
  459.         
  460.         /* Draw our three triangles */
  461.         p[0]=ap; p[1]=m[0]; p[2]=m[1];
  462.         _DrawTriangle(p);
  463.         p[0]=m[0]; p[1]=s[0]; p[2]=s[1];
  464.         _DrawTriangle(p);
  465.         p[0]=m[0]; p[1]=s[1]; p[2]=m[1];
  466.         _DrawTriangle(p);
  467.     } else {                                        /* All above water */
  468.         p[0].x=x1; p[0].y=y1; p[0].z=z[0];
  469.         p[1].x=x2; p[1].y=y2; p[1].z=z[1];
  470.         p[2].x=x3; p[2].y=y3; p[2].z=z[2];
  471.         _DrawTriangle(p);
  472.     }
  473. }
  474.  
  475. /* Draw the entire mountain range.  Non-recursive. */
  476. void DrawMountains()
  477. {
  478.     long x,y;
  479.     
  480.     if (map==0L) return;
  481.     SetCursor(*watch);
  482.     
  483.     /* Set drawing coefficients (used in CalcPoint3) */
  484.     xc=0.4073*(double)(1<<(maxIter-nIter));
  485.     sc=((double)wd)/630.0;
  486.     
  487.     /* Make transformation matrix */
  488.     tm[0][0]=1;
  489.     tm[1][0]=0;
  490.     tm[2][0]=0;
  491.     tm[0][1]=0;
  492.     tm[1][1]=cos(xTh);
  493.     tm[2][1]=sin(xTh);
  494.     tm[0][2]=0;
  495.     tm[1][2]=-sin(xTh);
  496.     tm[2][2]=cos(xTh);
  497.     
  498.     /* Go back to front, left to right, and draw each triangle */
  499.     for (y=0; y<mapS-1; y++) {
  500.         for (x=0; x<mapS-1; x++) {
  501.             DrawTriangle(x,y,x,y+1,x+1,y+1);
  502.             DrawTriangle(x,y,x+1,y+1,x+1,y);
  503.         }
  504.     }
  505.     
  506.     SetCursor(&arrow);
  507. }