1 /*
2 Copyright (c) 2014-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.jss;
30 
31 import dlib.math.vector;
32 
33 struct JohnsonSimplexSolver
34 {
35     void initialize()
36     {
37         foreach(ref v; supportPointsOnA)
38             v = Vector3f(0, 0, 0);
39 
40         foreach(ref v; supportPointsOnB)
41             v = Vector3f(0, 0, 0);
42 
43         foreach(ref e; edges)
44         foreach(ref v; e)
45             v = Vector3f(0, 0, 0);
46 
47         foreach(ref row; determinants)
48         foreach(ref v; row)
49             v = 0.0f;
50     }
51 
52     void calcClosestPoint(out Vector3f v)
53     {
54         float maxVertexSqrd = 0;
55         float deltaX = 0;
56 
57         Vector3f closestPoint = Vector3f(0, 0, 0);
58 
59         for (byte i = 0; i < 4; ++i)
60         {
61             if (isReducedSimplexPoint(i))
62             {
63                 float determinant = determinants[XBits][i];
64                 Vector3f point = edges[i][i];
65 
66                 assert(determinant > 0, "Negative or zero determinant found");
67 
68                 closestPoint += point * determinant;
69                 deltaX += determinant;
70 
71                 if (determinants[0][i] > maxVertexSqrd)
72                 {
73                     maxVertexSqrd = determinants[0][i];
74                     maxVertexIndex = i;
75                 }
76             }
77         }
78     
79         v = closestPoint / deltaX;
80     }
81 
82     void calcClosestPoint(byte set, out Vector3f v)
83     {
84         float deltaX = 0;
85 
86         Vector3f closestPoint = Vector3f(0, 0, 0);
87 
88         for (byte i = 0; i < 4; ++i)
89         {
90             if (set & (1 << i))
91             {
92                 float determinant = determinants[set][i];
93                 Vector3f point = edges[i][i];
94 
95                 assert(determinant > 0, "Negative or zero determinant found");
96 
97                 closestPoint += point * determinant;
98 
99                 deltaX += determinant;
100             }
101         }
102 
103         v = closestPoint / deltaX;
104     }
105 
106     void backupCalcClosestPoint(out Vector3f v)
107     {
108         // We don't need to update maxVertexIndex because this method is called when the algorithm terminates
109         float closestPointSqrd = float.max;
110 
111         for (byte subset = YBits; subset > 0; --subset)
112         {
113             if (isSubsetOrEqualTo(YBits, subset) && isProperSet(subset))
114             {
115                 Vector3f point;
116 
117                 calcClosestPoint(subset, point);
118 
119                 float pointSqrd = dot(point, point);
120                 if (pointSqrd < closestPointSqrd)
121                 {
122                     XBits = subset;
123 
124                     closestPointSqrd = pointSqrd;
125                     v = point;
126                 }
127             }
128         }
129     }
130 
131     void calcClosestPoints(out Vector3f pClosestPointOnA, out Vector3f pClosestPointOnB)
132     {
133         float deltaX = 0;
134 
135         Vector3f p = Vector3f(0,0,0);
136         Vector3f q = Vector3f(0,0,0);
137 
138         for (byte i = 0; i < 4; ++i)
139         {
140             if (isReducedSimplexPoint(i))
141             {
142                 float determinant = determinants[XBits][i];
143 
144                 assert(determinant > 0, "Negative or zero determinant found");
145 
146                 p += supportPointsOnA[i] * determinant;
147                 q += supportPointsOnB[i] * determinant;
148 
149                 deltaX += determinant;
150             }
151         }
152 
153         pClosestPointOnA = p / deltaX;
154         pClosestPointOnB = q / deltaX;
155     }
156 
157     bool reduceSimplex()
158     {
159         for (byte subset = YBits; subset > 0; --subset)
160         {
161             if ((subset & (1 << lastVertexIndex)) && isSubsetOrEqualTo(YBits, subset) && isValidSet(subset))
162             {
163                 XBits = subset;
164                 return true;
165             }
166         }
167 
168         // failed to reduce simplex
169         return false;
170     }
171 
172 
173     void addPoint(Vector3f point)
174     {
175         // Find a free slot for new point
176         byte freeSlot = findFreeSlot();
177         if (freeSlot == 4)
178         {
179             // Simplex is full
180             return;
181         }
182 
183         lastVertexIndex = freeSlot;
184         YBits = cast(byte)(XBits | (1 << lastVertexIndex));
185         edges[lastVertexIndex][lastVertexIndex] = point;
186 
187         // Update edges and determinants for new point
188         updateEdges(lastVertexIndex);
189         updateDeterminants(lastVertexIndex);
190 
191         // Update max vertex
192         float pointSqrd = dot(point, point);
193         if (pointSqrd > getMaxVertexSqrd())
194             maxVertexIndex = lastVertexIndex;
195 
196         determinants[0][lastVertexIndex] = pointSqrd;
197     }
198 
199     void addPoint(Vector3f point, Vector3f supportOnA, Vector3f supportOnB)
200     {
201         // Find a free slot for new point
202         byte freeSlot = findFreeSlot();
203         if (freeSlot == 4)
204         {
205             // Simplex is full
206             return;
207         }
208 
209         lastVertexIndex = freeSlot;
210         YBits = cast(byte)(XBits | (1 << lastVertexIndex));
211 
212         edges[lastVertexIndex][lastVertexIndex] = point;
213         supportPointsOnA[lastVertexIndex] = supportOnA;
214         supportPointsOnB[lastVertexIndex] = supportOnB;
215 
216         // Update edges and determinants for new point
217         updateEdges(lastVertexIndex);
218         updateDeterminants(lastVertexIndex);
219 
220         // Update max vertex
221         float pointSqrd = dot(point, point);
222         if (pointSqrd > getMaxVertexSqrd())
223             maxVertexIndex = lastVertexIndex;
224 
225         determinants[0][lastVertexIndex] = pointSqrd;
226     }
227 
228     void removeAllPoints()
229     {
230         XBits = 0;
231         YBits = 0;
232 
233         maxVertexIndex = 3;
234         lastVertexIndex = 3;
235     }
236 
237     void loadSimplex(byte simplexBits, Vector3f[] simplexPoints)
238     {
239         // Reset simplex
240         removeAllPoints();
241 
242         float maxPointSqrd = -float.max;
243 
244         // Add new vertices
245         for (byte i = 0; i < 4; ++i)
246         {
247             if (simplexBits & (1 << i))
248             {
249                 edges[i][i] = simplexPoints[i];
250                 determinants[0][i] = dot(edges[i][i], edges[i][i]);
251 
252                 if (determinants[0][i] > maxPointSqrd)
253                 {
254                     maxPointSqrd = determinants[0][i];
255                     maxVertexIndex = i;
256                 }
257 
258                 lastVertexIndex = i;
259             }
260         }
261   
262         assert(maxVertexIndex <= 3, "Invalid max vertex index");
263 
264         XBits = simplexBits;
265         YBits = simplexBits;
266 
267         // Calculate new edges based on new vertices
268         for (byte i = 0; i < 4; ++i)
269         {
270             if (simplexBits & (1 << i))
271                 updateEdges(i);
272         }
273 
274         // Calculate new determinants based on new vertices
275         for (byte i = 0; i < 4; ++i)
276         {
277             if (simplexBits & (1 << i))
278                 updateDeterminants(i);
279         }
280     }
281 
282     void loadSimplex(byte simplexBits, Vector3f[] simplexPoints, Vector3f[] supportOnA, Vector3f[] supportOnB)
283     {
284         loadSimplex(simplexBits, simplexPoints);
285 
286         for (byte i = 0; i < 4; ++i)
287         {
288             if (simplexBits & (1 << i))
289             {
290                 supportPointsOnA[i] = supportOnA[i];
291                 supportPointsOnB[i] = supportOnB[i];
292             }
293         }
294     }
295 
296     bool isReducedSimplexPoint(byte index)
297     {
298         assert(index < 4, "Invalid index");
299         return (XBits & (1 << index)) > 0;
300     }
301 
302     bool isReducedSimplexPoint(Vector3f point)
303     {
304         for (byte i = 0; i < 4; ++i)
305         {
306             if (isReducedSimplexPoint(i) && point == getPoint(i))
307                 return true;
308         }
309 
310         return false;
311     }
312 
313     bool isSimplexPoint(byte index)
314     {
315         assert(index < 4, "Invalid index");
316         return (YBits & (1 << index)) > 0;
317     }
318 
319     bool isSimplexPoint(Vector3f point)
320     {
321         for (byte i = 0; i < 4; ++i)
322         {
323             if (isSimplexPoint(i) && point == getPoint(i))
324                 return true;
325         }
326 
327         return false;
328     }
329 
330     Vector3f getPoint(byte index)
331     {
332         assert(index < 4, "Invalid index");
333         return edges[index][index];
334     }
335 
336     void setPoint(byte index, Vector3f point)
337     {
338         assert(index < 4, "Invalid index");
339         edges[index][index] = point;
340     }
341 
342     void setPoint(byte index, Vector3f point, Vector3f supportPointOnA, Vector3f supportPointOnB)
343     {
344         assert(index < 4, "Invalid index!");
345         edges[index][index] = point;
346         supportPointsOnA[index] = supportPointOnA;
347         supportPointsOnB[index] = supportPointOnB;
348     }
349 
350     byte getSimplex(Vector3f[] simplexPoints)
351     {
352         assert(simplexPoints.length >= 4, "Input array length must be greater or equal to 4");
353 
354         byte nbPoints = 0;
355 
356         for (byte i = 0; i < 4; ++i)
357         {
358             if (isReducedSimplexPoint(i))
359             {
360                 simplexPoints[nbPoints] = getPoint(i);
361                 ++nbPoints;
362             }
363         }
364 
365         return nbPoints;
366     }
367 
368     byte getSimplex(float[] simplexPoints)
369     {
370         assert(simplexPoints.length >= 12, "Input array length must be greater or equal to 12");
371 
372         byte nbPoints = 0;
373 
374         for (byte i = 0; i < 4; ++i)
375         {
376             if (isReducedSimplexPoint(i))
377             {
378                 Vector3f point = getPoint(i);
379                 simplexPoints[nbPoints * 3]     = point.x;
380                 simplexPoints[nbPoints * 3 + 1] = point.y;
381                 simplexPoints[nbPoints * 3 + 2] = point.z;
382 
383                 ++nbPoints;
384             }
385         }
386 
387        return nbPoints;
388     }
389 
390     byte getSupportPointsOnA(Vector3f[] worldSupportPointsOnA)
391     {
392         assert(supportPointsOnA.length >= 4, "Input array length must be greater or equal to 4");
393 
394         byte nbPoints = 0;
395 
396         for (byte i = 0; i < 4; ++i)
397         {
398             if (isReducedSimplexPoint(i))
399             {
400                 supportPointsOnA[nbPoints] = worldSupportPointsOnA[i];
401                 ++nbPoints;
402             }
403         }
404 
405         return nbPoints;
406     }
407 
408     byte getSupportPointsOnB(Vector3f[] worldSupportPointsOnB)
409     {
410         assert(supportPointsOnB.length >= 4, "Input array length must be greater or equal to 4");
411 
412         byte nbPoints = 0;
413 
414         for (byte i = 0; i < 4; ++i)
415         {
416             if (isReducedSimplexPoint(i))
417             {
418                 supportPointsOnB[nbPoints] = worldSupportPointsOnB[i];
419                 ++nbPoints;
420             }
421         }
422 
423         return nbPoints;
424     }
425 
426     Vector3f getSupportPointOnA(byte index)
427     {
428         assert(index < 4, "Invalid index");
429         return supportPointsOnA[index];
430     }
431 
432     Vector3f getSupportPointOnB(byte index)
433     {
434         assert(index < 4, "Invalid index");
435         return supportPointsOnB[index];
436     }
437 
438     bool isFullSimplex()
439     {
440         return (XBits == 0xf);
441     }
442 
443     bool isEmptySimplex()
444     {
445         return (XBits == 0);
446     }
447 
448     void setMaxVertexSqrd(float maxVertexSqrd)
449     {
450         determinants[0][maxVertexIndex] = maxVertexSqrd;
451     }
452 
453     float getMaxVertexSqrd()
454     {
455         return determinants[0][maxVertexIndex];
456     }
457 
458     void updateMaxVertex()
459     {
460         float maxVertexSqrd = -float.max;
461 
462         for (byte i = 0; i < 4; ++i)
463         {
464             if (isSimplexPoint(i))
465             {
466                 float pointSqrd = dot(edges[i][i], edges[i][i]);
467 
468                 if (pointSqrd > maxVertexSqrd)
469                 {
470                     maxVertexSqrd = pointSqrd;
471                     maxVertexIndex = i;
472                 }
473             }
474         }
475 
476         determinants[0][maxVertexIndex] = maxVertexSqrd;
477     }
478 
479     bool isSubsetOrEqualTo(byte set, byte subset)
480     {
481         return ((set & subset) == subset);
482     }
483 
484     bool isValidSet(byte set)
485     {
486         for (byte i = 0; i < 4; ++i)
487         {
488             if (isSimplexPoint(i))
489             {
490                 if (set & (1 << i))
491                 {
492                     // i-th point does belong to set
493                     if (determinants[set][i] <= 0)
494                         return false;
495                 }
496                 else
497                 {
498                     // i-th point does not belong to set
499                     if (determinants[set | (1 << i)][i] > 0)
500                         return false;
501                 }
502             }
503         }
504 
505         return true;
506     }
507 
508     bool isProperSet(byte set)
509     {
510         for (byte i = 0; i < 4; ++i)
511         {
512             if ((set & (1 << i)) && determinants[set][i] <= 0)
513                 return false;
514         }
515 
516         return true;
517     }
518 
519     byte findFreeSlot()
520     {
521         for (byte i = 0; i < 4; ++i)
522         {
523             if (!isReducedSimplexPoint(i))
524                 return i;
525         }
526 
527         // Invalid slot
528         return 4;
529     }
530 
531     void updateEdges(byte index)
532     {
533         assert(index < 4, "Invalid index");
534 
535         for (byte i = 0; i < 4; ++i)
536         {
537             if ((i != index) && isSimplexPoint(i))
538             {
539                 edges[index][i] =  edges[index][index] - edges[i][i];
540                 edges[i][index] = -edges[index][i];
541             }
542         }
543     }
544 
545     void updateDeterminants(byte index)
546     {
547         assert(index < 4, "Invalid index!");
548 
549         byte indexBit = cast(byte)(1 << index);
550         determinants[indexBit][index] = 1;
551 
552         // Update determinants for all subsets that contain a "valid" point at index
553 
554         for (byte i = 0; i < 4; ++i)
555         {
556             if ((i != index) && isReducedSimplexPoint(i))
557             {
558                 // calculate all determinants for subsets of combinations of 2 points including point at index
559                 byte subset2 = cast(byte)((1 << i) | indexBit);
560                 determinants[subset2][i] = dot(edges[index][i], getPoint(index));
561                 determinants[subset2][index] = dot(edges[i][index], getPoint(i));
562 
563                 for (byte j = 0; j < i; ++j)
564                 {
565                     if ((j != index) && isReducedSimplexPoint(j))
566                     {
567                         // calculate all determinants for subsets of combinations of 3 points including point at index
568                         byte subset3 = cast(byte)(subset2 | (1 << j));
569 
570                         determinants[subset3][j] = 
571                             dot(edges[i][j], getPoint(i)) * determinants[subset2][i] +
572                             dot(edges[i][j], getPoint(index)) * determinants[subset2][index];
573 
574                         determinants[subset3][i] = 
575                             dot(edges[j][i], getPoint(j)) * determinants[(1 << j) | indexBit][j] +
576                             dot(edges[j][i], getPoint(index)) * determinants[(1 << j) | indexBit][index];
577 
578                         determinants[subset3][index] = 
579                             dot(edges[j][index], getPoint(i)) * determinants[(1 << i) | (1 << j)][i] +
580                             dot(edges[j][index], getPoint(j)) * determinants[(1 << i) | (1 << j)][j];
581                     }
582                 }
583             }
584         }
585     
586         if (YBits == 0xf)
587         {
588             // compute determinants for full simplex
589 
590             determinants[YBits][0] = 
591                 dot(edges[1][0], getPoint(1)) * determinants[0xe][1] +
592                 dot(edges[1][0], getPoint(2)) * determinants[0xe][2] +
593                 dot(edges[1][0], getPoint(3)) * determinants[0xe][3];
594 
595             determinants[YBits][1] =
596                 dot(edges[0][1], getPoint(0)) * determinants[0xd][0] + 
597                 dot(edges[0][1], getPoint(2)) * determinants[0xd][2] +
598                 dot(edges[0][1], getPoint(3)) * determinants[0xd][3];
599 
600             determinants[YBits][2] = 
601                 dot(edges[0][2], getPoint(0)) * determinants[0xb][0] +
602                 dot(edges[0][2], getPoint(1)) * determinants[0xb][1] +
603                 dot(edges[0][2], getPoint(3)) * determinants[0xb][3];
604 
605             determinants[YBits][3] = 
606                 dot(edges[0][3], getPoint(0)) * determinants[0x7][0] +
607                 dot(edges[0][3], getPoint(1)) * determinants[0x7][1] +
608                 dot(edges[0][3], getPoint(2)) * determinants[0x7][2];
609         }
610     }
611 
612     // Reduced simplex bits.
613     byte XBits = 0;
614 
615     // Y = Wk + wk
616     byte YBits = 0;
617 
618     // Max vertex index
619     byte maxVertexIndex = 3;
620 
621     // Last vertex index
622     byte lastVertexIndex = 3;
623 
624     // Support points on shape A
625     Vector3f[4] supportPointsOnA;
626 
627     // Support points on shape B
628     Vector3f[4] supportPointsOnB;
629 
630     // Cached edges. The diagonal does not hold a subtraction, it holds all four simplex points instead
631     Vector3f[4][4] edges;
632 
633     // Cached determinants. This cache supports any combination of simplex's 4 points.
634     // The first column stores the squared length for each point
635     float[4][16] determinants;
636 }
637