1 /// Helper functions for working with geometrical shapes
2 module ggplotd.geometry;
3 
4 version( unittest )
5 {
6     import dunit.toolkit;
7 }
8 
9 struct Vertex3D
10 {
11     double x;
12     double y;
13     double z;
14 
15     this( in Vertex3D v )
16     {
17         x = v.x; y = v.y; z = v.z;
18     }
19 
20     this( in double _x, in double _y, in double _z )
21     {
22         x = _x; y = _y; z = _z;
23     }
24 
25     Vertex3D opBinary(string s)( in Vertex3D v2 ) const if (s == "-" || s == "+" )
26     {
27         mixin( "return Vertex3D( x " ~ s ~ " v2.x, y " ~ s ~ 
28             " v2.y, z " ~ s ~ " v2.z );");
29     }
30 }
31 unittest
32 {
33     auto v1 = Vertex3D(1,2,3);
34     auto v2 = Vertex3D(3,2,1);
35     auto v =  v1-v2;
36     assertEqual( v, Vertex3D(-2,0,2) );
37     v =  v1+v2;
38     assertEqual( v, Vertex3D(4,4,4) );
39 }
40 
41 Vertex3D crossProduct( in Vertex3D v1, in Vertex3D v2 )
42 {
43     return Vertex3D( v1.y*v2.z-v1.z*v2.y, 
44         v1.z*v2.x-v1.x*v2.z, 
45         v1.x*v2.y-v1.y*v2.x );
46 }
47 
48 unittest
49 {
50     auto v = Vertex3D(1,2,3).crossProduct( Vertex3D( 3,2,1 ) );
51     assertEqual( v, Vertex3D(-4,8,-4) );
52 }
53 
54 /// Calculate gradientVector based on a Triangle. The vertices of the triangle are assumed to be sorted by height.
55 auto gradientVector( T )( in T triangle )
56 {
57     import std.math : pow;
58     assert( triangle[0].z <= triangle[1].z && triangle[1].z <= triangle[2].z,
59         "gradientVector expects the triangle vertices to be sorted by height" );
60     auto gVector = [ Vertex3D(triangle[0]), Vertex3D(triangle[2]) ];
61 
62     if (triangle[0].z == triangle[2].z) {
63         return gVector;
64     }
65 
66     auto e1 = triangle[2]-triangle[0];
67     auto e2 = triangle[1]-triangle[0];
68     auto normal = e1.crossProduct( e2 );
69 
70     Vertex3D v;
71     if (normal.x == 0)
72     {
73         v = Vertex3D( normal.x/normal.y, 1, 
74             -(pow(normal.y,2)+pow(normal.x,2))/(normal.y*normal.z)
75             );
76     } else 
77     {
78         v = Vertex3D( 1, normal.y/normal.x, 
79             -(pow(normal.y,2)+pow(normal.x,2))/(normal.x*normal.z)
80             );
81     }
82 
83     auto scalar = (gVector[1].z-gVector[0].z)/v.z;
84 
85     gVector[1].x = gVector[0].x + scalar*v.x;
86     gVector[1].y = gVector[0].y + scalar*v.y;
87     gVector[1].z = gVector[0].z + scalar*v.z;
88 
89     return gVector; 
90 }
91 
92 unittest
93 {
94     import std.stdio;
95     auto v = gradientVector( [ Vertex3D( 0, 0, 1 ),
96         Vertex3D( 0, 1, 1 ),
97         Vertex3D( 1, 1, 1 )] );
98     assertEqual( v[0].z, v[1].z );
99 
100     v = gradientVector( [ Vertex3D( 0, 0, 0 ),
101         Vertex3D( 0, 1, 0.5 ),
102         Vertex3D( 1, 1, 1 )] );
103     assertEqual( v, [Vertex3D(0, 0, 0), Vertex3D(1, 1, 1)] );
104 
105     v = gradientVector( [ Vertex3D( 0, 0, 0 ),
106         Vertex3D( 0, 1, 0.1 ),
107         Vertex3D( 1, 1, 1 )] );
108     assertEqual( v[0], Vertex3D(0, 0, 0) );
109     assertApprox( v[1].x, 1.09756, 1e-5 );
110     assertApprox( v[1].y, 0.121951, 1e-5 );
111     assertApprox( v[1].z, 1.0 );
112 }
113