diff --git a/src/main/scala/scalismo/mesh/TetrahedralMesh.scala b/src/main/scala/scalismo/mesh/TetrahedralMesh.scala index 8dbb8318..3b789445 100644 --- a/src/main/scala/scalismo/mesh/TetrahedralMesh.scala +++ b/src/main/scala/scalismo/mesh/TetrahedralMesh.scala @@ -16,11 +16,11 @@ package scalismo.mesh import breeze.linalg.DenseVector -import scalismo.common._ -import scalismo.geometry._ +import scalismo.common.* +import scalismo.geometry.* +import scalismo.numerics.Determinant import scalismo.transformations.Transformation import scalismo.utils.Random -import vtk.vtkTetra import scala.language.implicitConversions @@ -160,8 +160,18 @@ case class TetrahedralMesh3D(pointSet: UnstructuredPoints[_3D], tetrahedralizati val c = pointSet.point(tetrahedron.ptId3) val d = pointSet.point(tetrahedron.ptId4) - // note: replace call to vtk with own implementation - val signedVolume = new vtkTetra().ComputeVolume(a.toArray, b.toArray, c.toArray, d.toArray) + val signedVolume = Determinant.det3x3( + b.x - a.x, + c.x - a.x, + d.x - a.x, + b.y - a.y, + c.y - a.y, + d.y - a.y, + b.z - a.z, + c.z - a.z, + d.z - a.z + ) / 6 + math.abs(signedVolume) } diff --git a/src/main/scala/scalismo/mesh/boundingSpheres/BSIntersection.scala b/src/main/scala/scalismo/mesh/boundingSpheres/BSIntersection.scala index bc73495b..83cebcd3 100644 --- a/src/main/scala/scalismo/mesh/boundingSpheres/BSIntersection.scala +++ b/src/main/scala/scalismo/mesh/boundingSpheres/BSIntersection.scala @@ -17,6 +17,7 @@ package scalismo.mesh.boundingSpheres import scalismo.geometry.{_3D, EuclideanVector, Point} import scalismo.mesh.{BarycentricCoordinates, BarycentricCoordinates4} +import scalismo.numerics.Determinant /** * Helper function to calculate intersection points. @@ -29,48 +30,48 @@ private[boundingSpheres] object BSIntersection { b: EuclideanVector[_3D], c: EuclideanVector[_3D] ): (Boolean, Point[_3D]) = { - val det = Determinantes.det3x3(a.x - b.x, - a.x - c.x, + val det = Determinant.det3x3(a.x - b.x, + a.x - c.x, + direction.x, + a.y - b.y, + a.y - c.y, + direction.y, + a.z - b.z, + a.z - c.z, + direction.z + ) + + val beta = Determinant.det3x3(a.x - point.x, + a.x - c.x, + direction.x, + a.y - point.y, + a.y - c.y, + direction.y, + a.z - point.z, + a.z - c.z, + direction.z + ) / det + + val gamma = Determinant.det3x3(a.x - b.x, + a.x - point.x, direction.x, a.y - b.y, - a.y - c.y, + a.y - point.y, direction.y, a.z - b.z, - a.z - c.z, + a.z - point.z, direction.z - ) - - val beta = Determinantes.det3x3(a.x - point.x, - a.x - c.x, - direction.x, - a.y - point.y, - a.y - c.y, - direction.y, - a.z - point.z, - a.z - c.z, - direction.z ) / det - val gamma = Determinantes.det3x3(a.x - b.x, - a.x - point.x, - direction.x, - a.y - b.y, - a.y - point.y, - direction.y, - a.z - b.z, - a.z - point.z, - direction.z - ) / det - - val t = Determinantes.det3x3(a.x - b.x, - a.x - c.x, - a.x - point.x, - a.y - b.y, - a.y - c.y, - a.y - point.y, - a.z - b.z, - a.z - c.z, - a.z - point.z + val t = Determinant.det3x3(a.x - b.x, + a.x - c.x, + a.x - point.x, + a.y - b.y, + a.y - c.y, + a.y - point.y, + a.z - b.z, + a.z - c.z, + a.z - point.z ) / det if (beta >= 0.0 && gamma >= 0.0 && beta + gamma <= 1.0) { @@ -87,37 +88,37 @@ private[boundingSpheres] object BSIntersection { b: EuclideanVector[_3D], c: EuclideanVector[_3D] ): (Boolean, BarycentricCoordinates) = { - val det = Determinantes.det3x3(a.x - b.x, - a.x - c.x, + val det = Determinant.det3x3(a.x - b.x, + a.x - c.x, + direction.x, + a.y - b.y, + a.y - c.y, + direction.y, + a.z - b.z, + a.z - c.z, + direction.z + ) + + val beta = Determinant.det3x3(a.x - point.x, + a.x - c.x, + direction.x, + a.y - point.y, + a.y - c.y, + direction.y, + a.z - point.z, + a.z - c.z, + direction.z + ) / det + + val gamma = Determinant.det3x3(a.x - b.x, + a.x - point.x, direction.x, a.y - b.y, - a.y - c.y, + a.y - point.y, direction.y, a.z - b.z, - a.z - c.z, + a.z - point.z, direction.z - ) - - val beta = Determinantes.det3x3(a.x - point.x, - a.x - c.x, - direction.x, - a.y - point.y, - a.y - c.y, - direction.y, - a.z - point.z, - a.z - c.z, - direction.z - ) / det - - val gamma = Determinantes.det3x3(a.x - b.x, - a.x - point.x, - direction.x, - a.y - b.y, - a.y - point.y, - direction.y, - a.z - b.z, - a.z - point.z, - direction.z ) / det if (beta >= 0.0 && gamma >= 0.0 && beta + gamma <= 1.0) { @@ -207,57 +208,3 @@ private[boundingSpheres] object BSIntersection { } } - -private[boundingSpheres] object Determinantes { - - @inline - def det2x2(a1: Double, a2: Double, b1: Double, b2: Double): Double = { - (+a1 * b2 - - b1 * a2) - } - - @inline - def det3x3(a1: Double, a2: Double, a3: Double, b1: Double, b2: Double, b3: Double, c1: Double, c2: Double, c3: Double) - : Double = { - (+a1 * det2x2(b2, b3, c2, c3) - - b1 * det2x2(a2, a3, c2, c3) - + c1 * det2x2(a2, a3, b2, b3)) - } - - @inline - def det4x4(a1: Double, - a2: Double, - a3: Double, - a4: Double, - b1: Double, - b2: Double, - b3: Double, - b4: Double, - c1: Double, - c2: Double, - c3: Double, - c4: Double, - d1: Double, - d2: Double, - d3: Double, - d4: Double - ): Double = { - (+a1 * det3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4) - - b1 * det3x3(c2, c3, c4, d2, d3, d4, a2, a3, a4) - + c1 * det3x3(d2, d3, d4, a2, a3, a4, b2, b3, b4) - - d1 * det3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4) - - a2 * det3x3(b3, b4, b1, c3, c4, c1, d3, d4, d1) - + b2 * det3x3(c3, c4, c1, d3, d4, d1, a3, a4, a1) - - c2 * det3x3(d3, d4, d1, a3, a4, a1, b3, b4, b1) - + d2 * det3x3(a3, a4, a1, b3, b4, b1, c3, c4, c1) - + a3 * det3x3(b4, b1, b2, c4, c1, c2, d4, d1, d2) - - b3 * det3x3(c4, c1, c2, d4, d1, d2, a4, a1, a2) - + c3 * det3x3(d4, d1, d2, a4, a1, a2, b4, b1, b2) - - d3 * det3x3(a4, a1, a2, b4, b1, b2, c4, c1, c2) - - a4 * det3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3) - + b4 * det3x3(c1, c2, c3, d1, d2, d3, a1, a2, a3) - - c4 * det3x3(d1, d2, d3, a1, a2, a3, b1, b2, b3) - + d4 * det3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3)) - } - -} diff --git a/src/main/scala/scalismo/numerics/Determinant.scala b/src/main/scala/scalismo/numerics/Determinant.scala new file mode 100644 index 00000000..417bfcef --- /dev/null +++ b/src/main/scala/scalismo/numerics/Determinant.scala @@ -0,0 +1,54 @@ +package scalismo.numerics + +object Determinant { + + @inline + def det2x2(a1: Double, a2: Double, b1: Double, b2: Double): Double = { + (+a1 * b2 + - b1 * a2) + } + + @inline + def det3x3(a1: Double, a2: Double, a3: Double, b1: Double, b2: Double, b3: Double, c1: Double, c2: Double, c3: Double) + : Double = { + (+a1 * det2x2(b2, b3, c2, c3) + - b1 * det2x2(a2, a3, c2, c3) + + c1 * det2x2(a2, a3, b2, b3)) + } + + @inline + def det4x4(a1: Double, + a2: Double, + a3: Double, + a4: Double, + b1: Double, + b2: Double, + b3: Double, + b4: Double, + c1: Double, + c2: Double, + c3: Double, + c4: Double, + d1: Double, + d2: Double, + d3: Double, + d4: Double + ): Double = { + (+a1 * det3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4) + - b1 * det3x3(c2, c3, c4, d2, d3, d4, a2, a3, a4) + + c1 * det3x3(d2, d3, d4, a2, a3, a4, b2, b3, b4) + - d1 * det3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4) + - a2 * det3x3(b3, b4, b1, c3, c4, c1, d3, d4, d1) + + b2 * det3x3(c3, c4, c1, d3, d4, d1, a3, a4, a1) + - c2 * det3x3(d3, d4, d1, a3, a4, a1, b3, b4, b1) + + d2 * det3x3(a3, a4, a1, b3, b4, b1, c3, c4, c1) + + a3 * det3x3(b4, b1, b2, c4, c1, c2, d4, d1, d2) + - b3 * det3x3(c4, c1, c2, d4, d1, d2, a4, a1, a2) + + c3 * det3x3(d4, d1, d2, a4, a1, a2, b4, b1, b2) + - d3 * det3x3(a4, a1, a2, b4, b1, b2, c4, c1, c2) + - a4 * det3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3) + + b4 * det3x3(c1, c2, c3, d1, d2, d3, a1, a2, a3) + - c4 * det3x3(d1, d2, d3, a1, a2, a3, b1, b2, b3) + + d4 * det3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3)) + } +}