diff --git a/Geometrica/Algebra/Matrix4.cs b/Geometrica/Algebra/Matrix4.cs
index baefa5c..3f89b80 100644
--- a/Geometrica/Algebra/Matrix4.cs
+++ b/Geometrica/Algebra/Matrix4.cs
@@ -1,4 +1,5 @@
-namespace Geometrica.Algebra
+
+namespace Geometrica.Algebra
 {
     public class Matrix4
     {
@@ -23,8 +24,8 @@ public Matrix4()
         }
 
         public Matrix4(
-            double m00, double m01, double m02, double m03, 
-            double m10, double m11, double m12, double m13, 
+            double m00, double m01, double m02, double m03,
+            double m10, double m11, double m12, double m13,
             double m20, double m21, double m22, double m23,
             double m30, double m31, double m32, double m33)
         {
@@ -66,14 +67,32 @@ public override string ToString()
 
         public double Determinant()
         {
-            return _matrix[0, 0] * _matrix[1, 1] * _matrix[2, 2] * _matrix[3, 3] +
-                   _matrix[1, 0] * _matrix[2, 1] * _matrix[3, 2] * _matrix[0, 3] +
-                   _matrix[2, 0] * _matrix[3, 1] * _matrix[0, 2] * _matrix[1, 3] +
-                   _matrix[3, 0] * _matrix[0, 1] * _matrix[1, 2] * _matrix[2, 3] -
-                   _matrix[0, 3] * _matrix[1, 2] * _matrix[2, 1] * _matrix[3, 0] -
-                   _matrix[0, 2] * _matrix[1, 1] * _matrix[2, 0] * _matrix[3, 3] -
-                   _matrix[0, 1] * _matrix[1, 0] * _matrix[2, 3] * _matrix[3, 2] -
-                   _matrix[0, 0] * _matrix[1, 3] * _matrix[2, 2] * _matrix[3, 1];
+            return _matrix[0, 0] * SubMatrix3(0, 0).Determinant() - _matrix[1, 0] * SubMatrix3(1, 0).Determinant() + _matrix[2, 0] * SubMatrix3(2, 0).Determinant() - _matrix[3, 0] * SubMatrix3(3, 0).Determinant();
+        }
+
+        public Matrix3 SubMatrix3(int r, int c)
+        {
+            var m3 = new Matrix3();
+            int x = 0, y = 0;
+
+            for (int i = 0; i < 4; i++)
+            {
+                y = 0;
+                if (i == r)
+                    continue;
+
+                for (int j = 0; j < 4; j++)
+                {
+                    if (j == c)
+                        continue;
+
+                    m3[x, y] = _matrix[i, j];
+                    y++;
+                }
+                x++;
+            }
+
+            return m3;
         }
     }
 }
diff --git a/Geometrica/DataStructures/DelaunayTriangulation.cs b/Geometrica/DataStructures/DelaunayTriangulation.cs
index 95ade73..974aaac 100644
--- a/Geometrica/DataStructures/DelaunayTriangulation.cs
+++ b/Geometrica/DataStructures/DelaunayTriangulation.cs
@@ -1,6 +1,7 @@
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices.ComTypes;
 using System.Xml.Linq;
+using Geometrica.Algebra;
 using Geometrica.Primitives;
 
 namespace Geometrica.DataStructures;
@@ -308,6 +309,17 @@ public static int GetNearestVertexIndexY(Triangle triangle, Point2 p)
 
         return distY2 < distY3 ? 1 : 2;
     }
+
+    public static bool InCircle(Point2 p1, Point2 p2, Point2 p3, Point2 p4)
+    {
+        var m = new Matrix4(
+            p1.X, p1.Y, p1.X * p1.X + p1.Y * p1.Y, 1,
+            p2.X, p2.Y, p2.X * p2.X + p2.Y * p2.Y, 1,
+            p3.X, p3.Y, p3.X * p3.X + p3.Y * p3.Y, 1,
+            p4.X, p4.Y, p4.X * p4.X + p4.Y * p4.Y, 1);
+
+        return m.Determinant() < -double.Epsilon;
+    }
 }
 
 public static class TriangleExtensions
diff --git a/GeometricaTests/DelaunayTriangulationTests.cs b/GeometricaTests/DelaunayTriangulationTests.cs
index d48f863..f3751cd 100644
--- a/GeometricaTests/DelaunayTriangulationTests.cs
+++ b/GeometricaTests/DelaunayTriangulationTests.cs
@@ -389,6 +389,28 @@ public void LegalizeTriangle_SwapsEdge_ForNonDelaunayPair()
             Assert.AreEqual("", actual);
         }
 
+        [TestMethod]
+        public void IncircleTest_Delaunay_ReturnsTrue()
+        {
+            var p1 = new Point2(0, 0);
+            var p2 = new Point2(1, 0);
+            var p3 = new Point2(0, 1);
+            var p4 = new Point2(1.1, 1.1);
+
+            Assert.IsTrue(DelaunayTriangulation.InCircle(p1, p2, p3, p4));
+        }
+
+        [TestMethod]
+        public void IncircleTest_NonDelaunay_ReturnsFalse()
+        {
+            var p1 = new Point2(0, 0);
+            var p2 = new Point2(1, 0);
+            var p3 = new Point2(0, 1);
+            var p4 = new Point2(0.9, 0.9);
+
+            Assert.IsFalse(DelaunayTriangulation.InCircle(p1, p2, p3, p4));
+        }
+
         private static Triangle[] PrepareRegularGrid(int resolutionX, int resolutionY)
         {
             var pts = new Point2[resolutionX, resolutionY];
diff --git a/GeometricaTests/Matrix4Tests.cs b/GeometricaTests/Matrix4Tests.cs
index 0e145ce..3e720e4 100644
--- a/GeometricaTests/Matrix4Tests.cs
+++ b/GeometricaTests/Matrix4Tests.cs
@@ -5,6 +5,35 @@ namespace GeometricaTests
     [TestClass]
     public class Matrix4Tests
     {
+        [TestMethod]
+        public void Submatrix3_RemovesGiven_ColAndRow()
+        {
+            var m = new Matrix4(
+                1, 2, 3, 4,
+                5, 6, 7, 8,
+                9, 10, 11, 12,
+                13, 14, 15, 16);
+
+            var m3 = m.SubMatrix3(0, 0);
+            var s = "[6;7;8][10;11;12][14;15;16]";
+
+            Assert.AreEqual(s, m3.ToString());
+        }
+
+        [TestMethod]
+        public void Matrix4_Returns_Negative022()
+        {
+            var m = new Matrix4(
+                0, 0, 0, 1,
+                1, 0, 1, 1,
+                0, 1, 1, 1,
+                1.1, 1.1, 2.42, 1);
+
+            var det = m.Determinant();
+
+            Assert.AreEqual(-0.22, det, 1e-3);
+        }
+
         [TestMethod]
         public void Matrix4_Members_CanBeEditedWithIndexers()
         {