#version 3.8; global_settings { assumed_gamma 1.0 } #include "math.inc" #declare E = 0.000001; #declare Camera_Orthographic = true; // <----################################################################################ #declare Camera_Position = <5, 5, -20>*10; // front view #declare Camera_Look_At = <5, 5, 0>*10 ; #declare Fraction = 6; // functions as a zoom for the orthographic view: 4 zooms in 4x, 8 zooms in 8x, etc. #declare Aspect = image_width/image_height; // ########################################### camera { #if (Camera_Orthographic = true) orthographic right x*image_width/(Fraction) up y*image_height/(Fraction) #else right x*image_width/image_height up y #end location Camera_Position look_at Camera_Look_At} // ########################################### sky_sphere {pigment {rgb 1}} light_source {<0, 0, -image_width/Fraction> color rgb 1} /*int *points_delaunay_naive_2d ( int n, double p[], int *ntri ) ****************************************************************************80 Purpose: POINTS_DELAUNAY_NAIVE_2D computes the Delaunay triangulation in 2D. Discussion: A naive and inefficient (but extremely simple) method is used. This routine is only suitable as a demonstration code for small problems. Its running time is of order N^4. Much faster algorithms are available. Given a set of nodes in the plane, a triangulation is a set of triples of distinct nodes, forming triangles, so that every point with the convex hull of the set of nodes is either one of the nodes, or lies on an edge of one or more triangles, or lies within exactly one triangle. Modified: 05 February 2005 Author: John Burkardt Reference: Joseph O'Rourke, Computational Geometry, Cambridge University Press, Second Edition, 1998, page 187. Parameters: Input, int N, the number of nodes. N must be at least 3. Input, double P[2*N], the coordinates of the nodes. Output, int *NTRI, the number of triangles. Output, int POINTS_DELAUNAY_NAIVE_2D[3*NTRI], the indices of the nodes making each triangle. { int count; int flag; int i; int j; int k; int m; int pass; int *tri; double xn; double yn; double zn; double *z; */ #macro Delaunay (n, p) #local Count = 0; #local Z = array [n]; #for (i, 0, n/2-1) //#debug concat ("i = ", str (i, 0, 0), "\n") #local Z[i] = p[0+i*2] * p[0+i*2] + p[1+i*2] * p[1+i*2]; #end // First pass counts triangles, // Second pass allocates triangles and sets them. #for (pass, 1, 2) #if ( pass = 2 ) #local tri = array [3*Count]; #end #local Count = 0; // For each triple (I,J,K): #for (i, 0, (n/2) - 2) //#for (i, 0, n - 2) #for (j, 0, n/2-1) #for (k, i+1, n/2-1) #if ( j != k ) #local xn = ( p[1+j*2] - p[1+i*2] ) * ( Z[k] - Z[i] ) - ( p[1+k*2] - p[1+i*2] ) * ( Z[j] - Z[i] ); #local yn = ( p[0+k*2] - p[0+i*2] ) * ( Z[j] - Z[i] ) - ( p[0+j*2] - p[0+i*2] ) * ( Z[k] - Z[i] ); #local zn = ( p[0+j*2] - p[0+i*2] ) * ( p[1+k*2] - p[1+i*2] ) - ( p[0+k*2] - p[0+i*2] ) * ( p[1+j*2] - p[1+i*2] ); #local flag = ( zn < 0 ); #if ( flag ) #for (m, 0, n/2-1) #local flag = (flag & ( ( p[0+m*2] - p[0+i*2] ) * xn + ( p[1+m*2] - p[1+i*2] ) * yn + ( Z[m] - Z[i] ) * zn <= 0 ) ); #end // end for #end // end if #if ( flag ) #if ( pass = 2 ) #local tri[0+Count*3] = i; #local tri[1+Count*3] = j; #local tri[2+Count*3] = k; #end // end if pass =2 #local Count = Count + 1; #end // end if flag #end // end if j != k #end // end for k #end // end for j #end // end for i #end // end for pass #declare NTriangles = Count; tri #end // end macro Delaunay //##################################################################################################### #declare Vertices = 40; #declare Seed = seed (123); #declare Points = array [Vertices] #declare NP = dimension_size (Points, 1); #for (i, 0, NP/2-1) #local Points[2*i] = floor (rand(Seed)*100); #local Points[2*i+1] = floor (rand(Seed)*100); #end #debug concat ("NP = ", str (NP, 0, 0), "\n") #declare Triangles = Delaunay (NP+1, Points); #declare Threshold = 0.01; //#macro VEq(V1, V2) (V1.x = V2.x & V1.y = V2.y & V1.z = V2.z) #macro VEq2 (V1, V2) (abs (V1.x - V2.x) < Threshold & abs (V1.y - V2.y) < Threshold & abs (V1.z - V2.z) < Threshold) #end #declare NumT = dimension_size (Triangles, 1); #debug concat ("NumT = ", str (NumT, 0, 0), "\n") #local NT = array; #local NT [0] = ; #local AT = 1; #while (AT <= NumT/3-1) #local Searching = 1; #local Match = 0; //#debug concat ("AT = ", str (AT, 0, 0), "\n") #local Tri = ; #local NTCount = dimension_size (NT, 1)-1; //#debug concat (" <", str (Triangles[0+AT*3], 0, 0), ", ", str (Triangles[1+AT*3], 0, 0), ", ", str (Triangles[2+AT*3], 0, 0), "> ", "\n") //#debug " ------- \n" #for (Test, 0, NTCount) //#debug concat (" <", str (NT[Test].x, 0, 0), ", ", str (NT[Test].y, 0, 0), ", ", str (NT[Test].z, 0, 0), "> ", "\n") //#debug concat (" Test = ", str (Test, 0, 0), "\n") // compare all permutations #local UTri1 = NT[Test]; #local UTri2 = ; #local UTri3 = ; #local UTri4 = ; #local UTri5 = ; #local UTri6 = ; #if (VEq2 (Tri, UTri1) | VEq (Tri, UTri2) | VEq (Tri, UTri3) | VEq (Tri, UTri4) | VEq (Tri, UTri5) | VEq (Tri, UTri6)) #local Match = 1; #end #end #if (Match) #else #local NT[NTCount+1] = Tri; //#debug " ^^^ *Unique triangle found* \n" #end #local AT = AT + 1; #end #local xmin = 0; #local xmax = 0; #local ymin = 0; #local ymax = 0; #for (i, 0, NP/2-1) #if (Points[2*i] > xmax) #local xmax = Points[2*i]; #elseif (Points[2*i] < xmin) #local xmin = Points[2*i]; #end #if (Points[2*i+1] > ymax) #local ymax = Points[2*i+1]; #elseif (Points[2*i+1] < ymin) #local ymin = Points[2*i+1]; #end #end // end for // draw box #local Line = 0.5; union { sphere {0, Line} cylinder {0, x*xmax, Line} sphere {x*xmax, Line} cylinder {x*xmax, x*xmax+y*ymax, Line} sphere {x*xmax+y*ymax, Line} cylinder {x*xmax+y*ymax, y*ymax, Line} sphere {y*ymax, Line} cylinder {y*ymax, 0, Line} } #for (P, 0, NP/2-1) sphere { Line*2 pigment {rgb 0}} #end #local Seed = seed (6745); #local E = 0.01; /* #for (T, 0, NumT/3-1) #local P1 = ; #local P2 = ; #local P3 = ; triangle {P1, P2, P3 pigment {rgb } translate -z*T*E} text { ttf "timrom.ttf", str (T*3, 0, 0), 0.02, 0.0 scale 5 translate (x+y)*(P1+P2+P3)/3 pigment {rgb 0} translate -z} // end of text object --------------------------------------------- #end */ #local NTCount = dimension_size (NT, 1)-1; #debug concat ("NTCount = ", str (NTCount, 0, 0), "\n") #for (T, 0, NTCount) #local P1 = ; #local P2 = ; #local P3 = ; triangle {P1, P2, P3 pigment {rgb }} text { ttf "timrom.ttf", str (T, 0, 0), 0.02, 0.0 scale 5 translate (x+y)*(P1+P2+P3)/3 pigment {rgb 0} translate -z*E} #debug "-------------------- \n" //#debug concat( "T", str (T, 0, 0), ": P1 = ", vstr(3, P1, ", ", 0, 5), " P2 = ", vstr(3, P2, ", ", 0, 5), " P3 = ", vstr(3, P3, ", ", 0, 5)," \n") #debug concat( "T", str (T, 0, 0), ",", vstr(3, P1, ", ", 0, 5), ",", vstr(3, P2, ", ", 0, 5), ",", vstr(3, P3, ", ", 0, 5)," \n") #end