//2D collage theorem demonstration, and linear IFS generator.
//By Dave Matthews
//This program does not yet bother to calculate the actual coefficients for the affine transformations,
//but it'd be easy enough to add that feature, and output to "debug" or a file. 
//Notes for changes necessary for extension to 3D IFS are included.   

#version 3.5;

#include "transforms.inc"


//These "switches" allow you to adjust the ball-and-stick triangles before generating the IFS

#declare Make_Triangles = 1;
#declare Make_IFS = 1;


///Here's the ball-and-stick parameters

#declare Ball_Rad = 0.1; //Destination sizes are by default 0.5*these
#declare Stick_Rad = 0.05;

///Here's the "dots" for the IFS

#declare Dot_Size = 0.05;
#declare Cursor = sphere {0, 1 pigment { color rgb x } 
                                scale Dot_Size no_shadow}; 
#declare N_Of_Dots = 20000;    


//Enter the source points.  This triangle gets remapped into the destination triangles
//Warning -- They must be non-collinear, to make a non-degenerate triangle.  This was not trapped for. 
// To convert to 3D IFS, you need a 4th (non-coplanar) source point, to make a non-degenrate tetrahedron.

#declare Source_Points = array[3];

#declare Source_Points[0] = <0, 0, 0>;
#declare Source_Points[1] = <16, 0, 0>;
#declare Source_Points[2] = <0, 16, 0>;   

//Enter Destination Points.  First decide on # of transforms

#declare Number_Of_Transforms = 8;

//Now enter each triad of destination points, followed by a relative weight, which should
//correspond roughly to the relative size of the destination triangle. These can add to 1, or if 
//they don't, they'll be normalized. 

//In order to get convergence, the destination points should define contraction mappings, 
//meaning that the destination triangles should be smaller than the source triangle.

// To convert to 3D IFS, each destination needs a 4th point.

#declare Destination_Points = array[Number_Of_Transforms][3]; 
#declare DP_Weight = array[Number_Of_Transforms];

#declare Destination_Points[0][0] = <0, 4.5, 0>;
#declare Destination_Points[0][1] = <0, 0, 0>;
#declare Destination_Points[0][2] = <0.5, 4.5, 0>;  

#declare DP_Weight[0] = 6;


#declare Destination_Points[1][0] = <1, 4, 0>;
#declare Destination_Points[1][1] = <3, 2, 0>;
#declare Destination_Points[1][2] = <1, 5, 0>;   

#declare DP_Weight[1] = 4;

#declare Destination_Points[2][0] = <1.5, 0, 0>;
#declare Destination_Points[2][1] = <3, 2, 0>;
#declare Destination_Points[2][2] = <1, 0.7, 0>;

#declare DP_Weight[2] = 4;  

#declare Destination_Points[3][0] = <4, 5, 0>;
#declare Destination_Points[3][1] = <4, 0, 0>;
#declare Destination_Points[3][2] = <5, 5, 0>;

#declare DP_Weight[3] = 6;

#declare Destination_Points[4][0] = <5, 4, 0>;
#declare Destination_Points[4][1] = <6, 2, 0>;
#declare Destination_Points[4][2] = <5, 5, 0>;

#declare DP_Weight[4] = 3;

#declare Destination_Points[5][0] = <7, 4, 0>;
#declare Destination_Points[5][1] = <6, 2, 0>;
#declare Destination_Points[5][2] = <7, 5, 0>;

#declare DP_Weight[5] = 3;

#declare Destination_Points[6][0] = <7, 5, 0>;
#declare Destination_Points[6][1] = <7, 0, 0>;
#declare Destination_Points[6][2] = <8, 5, 0>;

#declare DP_Weight[6] = 6; 

#declare Destination_Points[7][0] = <12, 4, 0>;
#declare Destination_Points[7][1] = <20, 12, 0>;
#declare Destination_Points[7][2] = <0, 16, 0>;

#declare DP_Weight[7] = 80;   

//Here's the normalization process:

#local TW = 0.0000000001;  //You can make this 0 if you make sure some DP's have non-zero weight

#local J = 0; 
#while (J < Number_Of_Transforms)
        #local TW = TW + DP_Weight[J];
        #local J = J + 1;
#end
#local J = 0;
#while (J < Number_Of_Transforms)
        #declare DP_Weight[J] = DP_Weight[J]/TW;
        #local J = J + 1;
#end

//////

// First, we make Ball-and-stick diagrams of the collage 
// In the 3D case, the extra vertex needs to bee connected, also
// (3 more sticks!) and you get a tetrahedron, rather than a triangle       

#macro DoTriangles()

#macro B_And_S_Triangle(V_Array,BR,SR)
        union {
        sphere { 0, 1  scale BR pigment { color rgb x } translate V_Array[0]}
        sphere { 0, 1  scale BR pigment { color rgb y } translate V_Array[1]}
        sphere { 0, 1  scale BR pigment { color rgb z } translate V_Array[2]}
        cylinder { V_Array[0], V_Array[1], SR pigment { color rgb x+y } } 
        cylinder { V_Array[0], V_Array[2], SR pigment { color rgb x+y } }
        cylinder { V_Array[1], V_Array[2], SR pigment { color rgb x+y } }  
        }  
#end



#macro Ball_And_Stick(N_Of_T,SArray,DArray,B_Rad,S_Rad)

#declare Source_Triangle = B_And_S_Triangle(SArray,B_Rad,S_Rad)
        
#declare Destination_Triangles = union {
        #local T_Count = 0; 
        #while (T_Count < N_Of_T)
                #local Temp_Array = array[3];
                #local J = 0;
                #while (J < 3)
                        #local Temp_Array[J] = DArray[T_Count][J];
                        #local J = J + 1;
                #end
        B_And_S_Triangle(Temp_Array,B_Rad/2,S_Rad/2)
        #local T_Count = T_Count + 1;
        #end
        }  
union {        
object { Source_Triangle }
object { Destination_Triangles translate <0, 0, -B_Rad> }  
        }

#end



object { Ball_And_Stick(Number_Of_Transforms,Source_Points,Destination_Points,Ball_Rad, Stick_Rad) }


#end

// Here, we determine the IFS matrices To do 3DIFS, add in the corresponding ".z" components in TransArray,
// as well as the K[3] vertex
     

#macro DoIFS()

#declare R1 = seed(12);


#macro TransArray(K)
  
#local A = K[1].x - K[0].x;
#local B = K[2].x - K[0].x;
#local C = K[1].y - K[0].y;
#local D = K[2].y - K[0].y;
#local E = K[0].x;
#local F = K[0].y;

transform { matrix 
            < A, C, 0,
              B, D, 0,
              0, 0, 1,
              E, F, 0 >
           } 
#end

#declare SMat = TransArray(Source_Points);
 
//Note that DMat is actually an array of matrices!

#declare DMat = array[Number_Of_Transforms];                                

#local J = 0;   
#declare Temp = array[3];
#while (J < Number_Of_Transforms) 
        #local I = 0;
        #while (I < 3)
          #declare Temp[I] = Destination_Points[J][I];
          #local I = I + 1;
        #end
        #declare DMat[J] = TransArray(Temp);
        #local J = J + 1;
#end


/////////////////////  This is the actual IFS machinery.  It should not need to change for 3D

                            

#declare IFS = union {
#local Old_V = <0,0,0>;
#local J = 0;
#while (J < N_Of_Dots )
        #local TD = rand(R1); 
        #local New_V = vinv_transform(Old_V,SMat);
        #local I = 0;
        #local P = 0;
        #local Flag = 0;
               #while ((Flag = 0) & (I < Number_Of_Transforms - 1))
                  #local P = P + DP_Weight[I];  
                  #if (TD < P)
                   #local New_V = vtransform(New_V,DMat[I]);
                   #local Flag = 1;
                  #else 
                   #local I = I + 1;
                  #end 
                #end
                #if (Flag = 0)
                  #local New_V = vtransform(New_V,DMat[Number_Of_Transforms-1]);
                #end 
        #if (J > 20)
        object {Cursor translate New_V}
        #end
        #local J = J + 1;
        #local Old_V = New_V;
#end
}


object {IFS }  

#end


// Check the switches

#if (Make_Triangles = 1)
   DoTriangles()
#end


#if (Make_IFS = 1)
  DoIFS()
#end



//Adjust the camera to get the whole IFS in the scene
//You probably won't want an orthographic view for 3D IFS

camera { orthographic
                             right     <4/3,0,0>*35
                             up        <  0,1,0>*35
                             location  < 0, 15,-40.0>
                             look_at   < 0, 15, 0.00>
                           }

light_source {<4, 3, -100> color rgb 1 }

// Change This If You Dislike The Checker Colors (Most People Do)
plane { z, 1 pigment { checker 
                        pigment {color rgb <0.95, 0.95, 1> } 
                        pigment { color rgb <1, 0.95, 0.95> } 
                      }
                      }