Add source of Orientation_6 for FPG

2025-04-24 09:19:00 +01:00 · 2025-04-24 09:19:00 +01:00 · 3d0a5e54a2
parent 90a376de17
commit 3d0a5e54a2
1 changed files with 211 additions and 0 deletions
--- a/NewKernel_d/include/CGAL/NewKernel_d/orientationC6.h
+++ b/NewKernel_d/include/CGAL/NewKernel_d/orientationC6.h
@ -0,0 +1,211 @@
+// Copyright (c) 2025  GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Andreas Fabri
+
+typedef double RT;
+
+int
+orientation(
+ RT a01,RT a02,RT a03,RT a04,RT a05,RT a06,
+ RT a11,RT a12,RT a13,RT a14,RT a15,RT a16,
+ RT a21,RT a22,RT a23,RT a24,RT a25,RT a26,
+ RT a31,RT a32,RT a33,RT a34,RT a35,RT a36,
+ RT a41,RT a42,RT a43,RT a44,RT a45,RT a46,
+ RT a51,RT a52,RT a53,RT a54,RT a55,RT a56,
+ RT a61,RT a62,RT a63,RT a64,RT a65,RT a66)
+{
+// First compute the det2x2
+  const RT m01 = a11 - a01;
+  const RT m02 = a21 - a01;
+  const RT m03 = a31 - a01;
+  const RT m04 = a41 - a01;
+  const RT m05 = a51 - a01;
+  const RT m06 = a61 - a01;
+
+  const RT m12 = a21 - a11;
+  const RT m13 = a31 - a11;
+  const RT m14 = a41 - a11;
+  const RT m15 = a51 - a11;
+  const RT m16 = a61 - a11;
+
+  const RT m23 = a31 - a21;
+  const RT m24 = a41 - a21;
+  const RT m25 = a51 - a21;
+  const RT m26 = a61 - a21;
+
+  const RT m34 = a41 - a31;
+  const RT m35 = a51 - a31;
+  const RT m36 = a61 - a31;
+
+  const RT m45 = a51 - a41;
+  const RT m46 = a61 - a41;
+
+  const RT m56 = a61 - a51;
+
+// Now compute the minors of rank 3
+  const RT m012 = m01*a22 - m02*a12 + m12*a02;
+  const RT m013 = m01*a32 - m03*a12 + m13*a02;
+  const RT m014 = m01*a42 - m04*a12 + m14*a02;
+  const RT m015 = m01*a52 - m05*a12 + m15*a02;
+  const RT m016 = m01*a62 - m06*a12 + m16*a02;
+
+  const RT m023 = m02*a32 - m03*a22 + m23*a02;
+  const RT m024 = m02*a42 - m04*a22 + m24*a02;
+  const RT m025 = m02*a52 - m05*a22 + m25*a02;
+  const RT m026 = m02*a62 - m06*a22 + m26*a02;
+  const RT m034 = m03*a42 - m04*a32 + m34*a02;
+  const RT m035 = m03*a52 - m05*a32 + m35*a02;
+  const RT m036 = m03*a62 - m06*a32 + m36*a02;
+
+  const RT m045 = m04*a52 - m05*a42 + m45*a02;
+  const RT m046 = m04*a62 - m06*a42 + m46*a02;
+
+  const RT m056 = m05*a62 - m06*a52 + m56*a02;
+
+
+  const RT m123 = m12*a32 - m13*a22 + m23*a12;
+  const RT m124 = m12*a42 - m14*a22 + m24*a12;
+  const RT m125 = m12*a52 - m15*a22 + m25*a12;
+  const RT m126 = m12*a62 - m16*a22 + m26*a12;
+
+  const RT m134 = m13*a42 - m14*a32 + m34*a12;
+  const RT m135 = m13*a52 - m15*a32 + m35*a12;
+  const RT m136 = m13*a62 - m16*a32 + m36*a12;
+
+  const RT m145 = m14*a52 - m15*a42 + m45*a12;
+  const RT m146 = m14*a62 - m16*a42 + m46*a12;
+
+  const RT m156 = m15*a62 - m16*a52 + m56*a12;
+
+  const RT m234 = m23*a42 - m24*a32 + m34*a22;
+  const RT m235 = m23*a52 - m25*a32 + m35*a22;
+  const RT m236 = m23*a62 - m26*a32 + m36*a22;
+
+  const RT m245 = m24*a52 - m25*a42 + m45*a22;
+  const RT m246 = m24*a62 - m26*a42 + m46*a22;
+
+  const RT m256 = m25*a62 - m26*a52 + m56*a22;
+
+  const RT m345 = m34*a52 - m35*a42 + m45*a32;
+  const RT m346 = m34*a62 - m36*a42 + m46*a32;
+
+  const RT m356 = m35*a62 - m36*a52 + m56*a32;
+
+  const RT m456 = m45*a62 - m46*a52 + m56*a42;
+
+// Now compute the minors of rank 4
+  const RT m0123 = m012*a33 - m013*a23 + m023*a13 - m123*a03;
+
+  const RT m0124 = m012*a43 - m014*a23 + m024*a13 - m124*a03;
+  const RT m0125 = m012*a53 - m015*a23 + m025*a13 - m125*a03;
+  const RT m0126 = m012*a63 - m016*a23 + m026*a13 - m126*a03;
+
+  const RT m0134 = m013*a43 - m014*a33 + m034*a13 - m134*a03;
+  const RT m0135 = m013*a53 - m015*a33 + m035*a13 - m135*a03;
+  const RT m0136 = m013*a63 - m016*a33 + m036*a13 - m136*a03;
+
+  const RT m0145 = m014*a53 - m015*a43 + m045*a13 - m145*a03;
+  const RT m0146 = m014*a63 - m016*a43 + m046*a13 - m146*a03;
+
+  const RT m0156 = m015*a63 - m016*a53 + m056*a13 - m156*a03;
+
+  const RT m0234 = m023*a43 - m024*a33 + m034*a23 - m234*a03;
+  const RT m0235 = m023*a53 - m025*a33 + m035*a23 - m235*a03;
+  const RT m0236 = m023*a63 - m026*a33 + m036*a23 - m236*a03;
+
+  const RT m0245 = m024*a53 - m025*a43 + m045*a23 - m245*a03;
+  const RT m0246 = m024*a63 - m026*a43 + m046*a23 - m246*a03;
+
+  const RT m0256 = m025*a63 - m026*a53 + m056*a23 - m256*a03;
+
+  const RT m0345 = m034*a53 - m035*a43 + m045*a33 - m345*a03;
+  const RT m0346 = m034*a63 - m036*a43 + m046*a33 - m346*a03;
+
+  const RT m0356 = m035*a63 - m036*a53 + m056*a33 - m356*a03;
+
+  const RT m0456 = m045*a63 - m046*a53 + m056*a43 - m456*a03;
+
+  const RT m1234 = m123*a43 - m124*a33 + m134*a23 - m234*a13;
+  const RT m1235 = m123*a53 - m125*a33 + m135*a23 - m235*a13;
+  const RT m1236 = m123*a63 - m126*a33 + m136*a23 - m236*a13;
+
+  const RT m1245 = m124*a53 - m125*a43 + m145*a23 - m245*a13;
+  const RT m1246 = m124*a63 - m126*a43 + m146*a23 - m246*a13;
+
+  const RT m1256 = m125*a63 - m126*a53 + m156*a23 - m256*a13;
+
+  const RT m1345 = m134*a53 - m135*a43 + m145*a33 - m345*a13;
+  const RT m1346 = m134*a63 - m136*a43 + m146*a33 - m346*a13;
+
+  const RT m1356 = m135*a63 - m136*a53 + m156*a33 - m356*a13;
+  const RT m1456 = m145*a63 - m146*a53 + m156*a43 - m456*a13;
+
+  const RT m2345 = m234*a53 - m235*a43 + m245*a33 - m345*a23;
+  const RT m2346 = m234*a63 - m236*a43 + m246*a33 - m346*a23;
+
+  const RT m2356 = m235*a63 - m236*a53 + m256*a33 - m356*a23;
+  const RT m2456 = m245*a63 - m246*a53 + m256*a43 - m456*a23;
+
+  const RT m3456 = m345*a63 - m346*a53 + m356*a43 - m456*a33;
+
+
+  // Now compute the minors of rank 5
+  const RT m01234 = m0123*a44 - m0124*a34 + m0134*a24 - m0234*a14 + m1234*a04;
+
+  const RT m01235 = m0123*a54 - m0125*a34 + m0135*a24 - m0235*a14 + m1235*a04;
+
+  const RT m01236 = m0123*a64 - m0126*a34 + m0136*a24 - m0236*a14 + m1236*a04;
+
+  const RT m01245 = m0124*a54 - m0125*a44 + m0145*a24 - m0245*a14 + m1245*a04;
+  const RT m01246 = m0124*a64 - m0126*a44 + m0146*a24 - m0246*a14 + m1246*a04;
+
+  const RT m01256 = m0125*a64 - m0126*a54 + m0156*a24 - m0256*a14 + m1256*a04;
+
+  const RT m01345 = m0134*a54 - m0135*a44 + m0145*a34 - m0345*a14 + m1345*a04;
+  const RT m01346 = m0134*a64 - m0136*a44 + m0146*a34 - m0346*a14 + m1346*a04;
+
+  const RT m01356 = m0135*a64 - m0136*a54 + m0156*a34 - m0356*a14 + m1356*a04;
+  const RT m01456 = m0145*a64 - m0146*a54 + m0156*a44 - m0456*a14 + m1456*a04;
+
+  const RT m02345 = m0234*a54 - m0235*a44 + m0245*a34 - m0345*a24 + m2345*a04;
+  const RT m02346 = m0234*a64 - m0236*a44 + m0246*a34 - m0346*a24 + m2346*a04;
+
+  const RT m02356 = m0235*a64 - m0236*a54 + m0256*a34 - m0356*a24 + m2356*a04;
+  const RT m02456 = m0245*a64 - m0246*a54 + m0256*a44 - m0456*a24 + m2456*a04;
+  const RT m03456 = m0345*a64 - m0346*a54 + m0356*a44 - m0456*a34 + m3456*a04;
+
+  const RT m12345 = m1234*a54 - m1235*a44 + m1245*a34 - m1345*a24 + m2345*a14;
+  const RT m12346 = m1234*a64 - m1236*a44 + m1246*a34 - m1346*a24 + m2346*a14;
+
+
+  const RT m12356 = m1235*a64 - m1236*a54 + m1256*a34 - m1356*a24 + m2356*a14;
+  const RT m12456 = m1245*a64 - m1246*a54 + m1256*a44 - m1456*a24 + m2456*a14;
+  const RT m13456 = m1345*a64 - m1346*a54 + m1356*a44 - m1456*a34 + m3456*a14;
+
+  const RT m23456 = m2345*a64 - m2346*a54 + m2356*a44 - m2456*a34 + m3456*a24;
+
+// Now compute the minors of rank 6
+  const RT m012345 = m01234*a55 - m01235*a45 + m01245*a35 - m01345*a25 + m02345*a15 - m12345*a05;
+  const RT m012346 = m01234*a65 - m01236*a45 + m01246*a35 - m01346*a25 + m02346*a15 - m12346*a05;
+  const RT m012356 = m01235*a65 - m01236*a55 + m01256*a35 - m01356*a25 + m02356*a15 - m12356*a05;
+  const RT m012456 = m01245*a65 - m01246*a55 + m01256*a45 - m01456*a25 + m02456*a15 - m12456*a05;
+  const RT m013456 = m01345*a65 - m01346*a55 + m01356*a45 - m01456*a35 + m03456*a15 - m13456*a05;
+  const RT m023456 = m02345*a65 - m02346*a55 + m02356*a45 - m02456*a35 + m03456*a25 - m23456*a05;
+  const RT m123456 = m12345*a65 - m12346*a55 + m12356*a45 - m12456*a35 + m13456*a25 - m23456*a15;
+
+
+  // Now compute the minors of rank 7
+  const RT m0123456 = m012345 * a66 - m012346 * a56 + m012356 * a46 -  m012456 * a36 + m013456 * a26 - m023456 * a16 + m123456 * a06;
+
+  if(m0123456 > 0) return 1;
+  if(m0123456 < 0) return -1;
+  return 0;
+}