1 /*
2 Copyright (c) 2013-2017 Timur Gafarov
3 
4 Boost Software License - Version 1.0 - August 17th, 2003
5 
6 Permission is hereby granted, free of charge, to any person or organization
7 obtaining a copy of the software and accompanying documentation covered by
8 this license (the "Software") to use, reproduce, display, distribute,
9 execute, and transmit the Software, and to prepare derivative works of the
10 Software, and to permit third-parties to whom the Software is furnished to
11 do so, all subject to the following:
12 
13 The copyright notices in the Software and this entire statement, including
14 the above license grant, this restriction and the following disclaimer,
15 must be included in all copies of the Software, in whole or in part, and
16 all derivative works of the Software, unless such copies or derivative
17 works are solely in the form of machine-executable object code generated by
18 a source language processor.
19 
20 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
21 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
23 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
24 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
25 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 DEALINGS IN THE SOFTWARE.
27 */
28 
29 module dmech.mpr;
30 
31 import dlib.math.vector;
32 import dlib.math.matrix;
33 import dlib.math.transformation;
34 import dlib.math.utils;
35 
36 import dmech.shape;
37 import dmech.contact;
38 
39 /*
40  * Implementation of the Minkowski Portal Refinement algorithm
41  */
42 
43 // TODO: use ShapeComponent.supportPointGlobal instead of this
44 void supportTransformed(ShapeComponent s, Vector3f dir, out Vector3f result)
45 {
46     Matrix4x4f m = s.transformation;
47 
48     result.x = ((dir.x * m.a11) + (dir.y * m.a21)) + (dir.z * m.a31);
49     result.y = ((dir.x * m.a12) + (dir.y * m.a22)) + (dir.z * m.a32);
50     result.z = ((dir.x * m.a13) + (dir.y * m.a23)) + (dir.z * m.a33);
51 
52     result = s.geometry.supportPoint(result);
53 
54     float x = ((result.x * m.a11) + (result.y * m.a12)) + (result.z * m.a13);
55     float y = ((result.x * m.a21) + (result.y * m.a22)) + (result.z * m.a23);
56     float z = ((result.x * m.a31) + (result.y * m.a32)) + (result.z * m.a33);
57 
58     result.x = m.a14 + x;
59     result.y = m.a24 + y;
60     result.z = m.a34 + z;
61 }
62 
63 /*
64  * TODO:
65  * - write c.fact here
66  */
67 
68 bool MPRCollisionTest(
69     ShapeComponent s1,
70     ShapeComponent s2,
71     ref Contact c)
72 {
73     enum float collideEpsilon = 1e-4f;
74     enum maxIterations = 10;
75 
76     // Used variables
77     Vector3f temp1;
78     Vector3f v01, v02, v0;
79     Vector3f v11, v12, v1;
80     Vector3f v21, v22, v2;
81     Vector3f v31, v32, v3;
82     Vector3f v41, v42, v4;
83 
84     Matrix4x4f transform1 = s1.transformation;
85     Matrix4x4f transform2 = s2.transformation;
86 
87     // Initialization of the output
88     c.point = Vector3f(0.0f, 0.0f, 0.0f);
89     c.normal = Vector3f(0.0f, 0.0f, 0.0f);
90     c.penetration = 0.0f;
91 
92     // Get the center of shape1 in world coordinates
93     v01 = transform1.translation;
94 
95     // Get the center of shape2 in world coordinates
96     v02 = transform2.translation;
97 
98     // v0 is the center of the Minkowski difference
99     v0 = v02 - v01;
100 
101     // Avoid case where centers overlap - any direction is fine in this case
102     if (v0.isAlmostZero)
103         v0 = Vector3f(0.00001f, 0.0f, 0.0f);
104 
105     // v1 = support in direction of origin
106     c.normal = -v0;
107 
108     supportTransformed(s1, v0, v11);
109     supportTransformed(s2, c.normal, v12);
110     v1 = v12 - v11;
111 
112     if (dot(v1, c.normal) <= 0.0f)
113         return false;
114 
115     // v2 = support perpendicular to v1,v0
116     c.normal = cross(v1, v0);
117 
118     if (c.normal.isAlmostZero)
119     {
120         c.normal = v1 - v0;
121         c.normal.normalize();
122 
123         c.point = v11;
124         c.point += v12;
125         c.point *= 0.5f;
126 
127         c.penetration = dot(v12 - v11, c.normal);
128 
129         return true;
130     }
131 
132     supportTransformed(s1, -c.normal, v21);
133     supportTransformed(s2,  c.normal, v22);
134     v2 = v22 - v21;
135 
136     if (dot(v2, c.normal) <= 0.0f)
137         return false;
138 
139     // Determine whether origin is on + or - side of plane (v1,v0,v2)
140     c.normal = cross(v1 - v0, v2 - v0);
141 
142     float dist = dot(c.normal, v0);
143 
144     // If the origin is on the - side of the plane, reverse the direction of the plane
145     if (dist > 0.0f)
146     {
147         swap(&v1, &v2);
148         swap(&v11, &v21);
149         swap(&v12, &v22);
150         c.normal = -c.normal;
151     }
152 
153     int phase2 = 0;
154     int phase1 = 0;
155     bool hit = false;
156 
157     // Phase One: Identify a portal
158     while (true)
159     {
160         if (phase1 > maxIterations)
161             return false;
162 
163         phase1++;
164 
165         // Obtain the support point in a direction perpendicular to the existing plane
166         // Note: This point is guaranteed to lie off the plane
167         supportTransformed(s1, -c.normal, v31);
168         supportTransformed(s2,  c.normal, v32);
169         v3 = v32 - v31;
170 
171         if (dot(v3, c.normal) <= 0.0f)
172             return false;
173 
174         // If origin is outside (v1,v0,v3), then eliminate v2 and loop
175         temp1 = cross(v1, v3);
176         if (dot(temp1, v0) < 0.0f)
177         {
178             v2 = v3;
179             v21 = v31;
180             v22 = v32;
181             c.normal = cross(v1 - v0, v3 - v0);
182             continue;
183         }
184 
185         // If origin is outside (v3,v0,v2), then eliminate v1 and loop
186         temp1 = cross(v3, v2);
187         if (dot(temp1, v0) < 0.0f)
188         {
189             v1 = v3;
190             v11 = v31;
191             v12 = v32;
192             c.normal = cross(v3 - v0, v2 - v0);
193             continue;
194         }
195 
196         // Phase Two: Refine the portal
197         // We are now inside of a wedge...
198         while (true)
199         {
200             phase2++;
201 
202             // Compute normal of the wedge face
203             c.normal = cross(v2 - v1, v3 - v1);
204 
205             // Can this happen? Can it be handled more cleanly?
206             //if (c.normal.isAlmostZero)
207             //    return true;
208 
209             c.normal.normalize();
210 
211             // Compute distance from origin to wedge face
212             float d = dot(c.normal, v1);
213 
214             // If the origin is inside the wedge, we have a hit
215             if (d >= 0 && !hit)
216                 hit = true;
217 
218             // Find the support point in the direction of the wedge face
219             supportTransformed(s1, -c.normal, v41);
220             supportTransformed(s2,  c.normal, v42);
221             v4 = v42 - v41;
222 
223             float delta = dot(v4 - v3, c.normal);
224             c.penetration = dot(v4, c.normal);
225 
226             // If the boundary is thin enough or the origin is outside
227             // the support plane for the newly discovered vertex, then we can terminate
228             if (delta <= collideEpsilon || c.penetration <= 0.0f || phase2 > maxIterations)
229             {
230                 if (hit)
231                 {
232                     float b0 = dot(cross(v1, v2), v3);
233                     float b1 = dot(cross(v3, v2), v0);
234                     float b2 = dot(cross(v0, v1), v3);
235                     float b3 = dot(cross(v2, v1), v0);
236 
237                     float sum = b0 + b1 + b2 + b3;
238 
239                     if (sum <= 0)
240                     {
241                         b0 = 0;
242                         b1 = dot(cross(v2, v3), c.normal);
243                         b2 = dot(cross(v3, v1), c.normal);
244                         b3 = dot(cross(v1, v2), c.normal);
245 
246                         sum = b1 + b2 + b3;
247                     }
248 
249                     float inv = 1.0f / sum;
250 
251                     c.point = v01 * b0;
252                     c.point += v11 * b1;
253                     c.point += v21 * b2;
254                     c.point += v31 * b3;
255 
256                     c.point += v02 * b0;
257                     c.point += v12 * b1;
258                     c.point += v22 * b2;
259                     c.point += v32 * b3;
260 
261                     c.point *= inv * 0.5f;
262                 }
263 
264                 return hit;
265             }
266 
267             // Compute the tetrahedron dividing face (v4,v0,v3)
268             temp1 = cross(v4, v0);
269             float d2 = dot(temp1, v1);
270 
271             if (d2 >= 0.0f)
272             {
273                 d2 = dot(temp1, v2);
274                 if (d2 >= 0.0f)
275                 {
276                     // Inside d1 & inside d2 ==> eliminate v1
277                     v1 = v4;
278                     v11 = v41;
279                     v12 = v42;
280                 }
281                 else
282                 {
283                     // Inside d1 & outside d2 ==> eliminate v3
284                     v3 = v4;
285                     v31 = v41;
286                     v32 = v42;
287                 }
288             }
289             else
290             {
291                 d2 = dot(temp1, v3);
292 
293                 if (d2 >= 0.0f)
294                 {
295                     // Outside d1 & inside d3 ==> eliminate v2
296                     v2 = v4;
297                     v21 = v41;
298                     v22 = v42;
299                 }
300                 else
301                 {
302                     // Outside d1 & outside d3 ==> eliminate v1
303                     v1 = v4;
304                     v11 = v41;
305                     v12 = v42;
306                 }
307             }
308         }
309     }
310 
311     // Should never get here
312     assert(0);
313 }