7

I have a two dimensional euclidean space. Three points are given.

For example (p2 is the middle point):

Point2D p1 = new Point2D.Double(177, 289);
Point2D p2 = new Point2D.Double(178, 290);
Point2D p3 = new Point2D.Double(178, 291);

Now i want to calculate the curvature for these three points.

double curvature = calculateCurvature(p1, p2, p3);

How to do this? Ist there a existing method (no java external libraries)?

Spenhouet
  • 6,556
  • 12
  • 51
  • 76

5 Answers5

13

For the Menger Curvature, the formula is right there in the Wikipedia article :

curvature = 4*triangleArea/(sideLength1*sideLength2*sideLength3)

Which code did you try exactly?

It shouldn't be too hard to calculate those 4 values given your 3 points.

Here are some helpful methods : 

/**
 * Returns twice the signed area of the triangle a-b-c.
 * @param a first point
 * @param b second point
 * @param c third point
 * @return twice the signed area of the triangle a-b-c
 */
public static double area2(Point2D a, Point2D b, Point2D c) {
    return (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x);
}

/**
 * Returns the Euclidean distance between this point and that point.
 * @param that the other point
 * @return the Euclidean distance between this point and that point
 */
public double distanceTo(Point2D that) {
    double dx = this.x - that.x;
    double dy = this.y - that.y;
    return Math.sqrt(dx*dx + dy*dy);
}

There's not much more to do. Warning : area2 returns a signed double, depending on the orientation of your points (clockwise or anticlockwise).

Eric Duminil
  • 52,989
  • 9
  • 71
  • 124
  • I don't know what Point2D that is, but the Java Point2D doesn't have these methods: https://docs.oracle.com/javase/7/docs/api/java/awt/geom/Point2D.html – Spenhouet Dec 14 '16 at 14:41
  • So define them, and replace `b.x` with `b.getX()` ;) – Eric Duminil Dec 14 '16 at 14:44
  • I'm about to do so. Thank you :) – Spenhouet Dec 14 '16 at 14:44
  • 1
    I'd like to note also that Heron's formula is a numerically stable way to get the (unsigned) area of the triangle. In fact, it naturally gives you 4*area, and its inputs are the lengths of ths sides -- exactly what you already need for the demoninator of the menger curvature. – Paul Du Bois Jan 09 '19 at 19:37
5

As already pointed out by Eric Duminil in his answer, the computation is

curvature = 4*triangleArea/(sideLength0*sideLength1*sideLength2)

I wasted some time with creating this interactive example that contains a computeCurvature method that does the whole computation at once:

Curvature

import java.awt.Color;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.event.MouseEvent;
import java.awt.event.MouseListener;
import java.awt.event.MouseMotionListener;
import java.awt.geom.Ellipse2D;
import java.awt.geom.Point2D;
import java.util.ArrayList;
import java.util.List;

import javax.swing.JFrame;
import javax.swing.JPanel;
import javax.swing.SwingUtilities;

public class CurvatureFromThreePoints
{
    public static void main(String[] args)
    {
        SwingUtilities.invokeLater(new Runnable()
        {
            @Override
            public void run()
            {
                createAndShowGUI();
            }
        });
    }
    
    private static void createAndShowGUI()
    {
        JFrame f = new JFrame();
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.getContentPane().add(new CurvatureFromThreePointsPanel());
        f.setSize(800,800);
        f.setLocationRelativeTo(null);
        f.setVisible(true);
    }

}

class CurvatureFromThreePointsPanel extends JPanel 
    implements MouseListener, MouseMotionListener
{
    private final List<Point2D> pointList;
    private Point2D draggedPoint;
    
    public CurvatureFromThreePointsPanel()
    {
        this.pointList = new ArrayList<Point2D>();
        
        pointList.add(new Point2D.Double(132,532));
        pointList.add(new Point2D.Double(275,258));
        pointList.add(new Point2D.Double(395,267));

        addMouseListener(this);
        addMouseMotionListener(this);
    }
    
    private static double computeCurvature(Point2D p0, Point2D p1, Point2D p2)
    {
        double dx1 = p1.getX() - p0.getX();
        double dy1 = p1.getY() - p0.getY();
        double dx2 = p2.getX() - p0.getX();
        double dy2 = p2.getY() - p0.getY();
        double area = 0.5 * (dx1 * dy2 - dy1 * dx2;
        double len0 = p0.distance(p1);
        double len1 = p1.distance(p2);
        double len2 = p2.distance(p0);
        return 4 * area / (len0 * len1 * len2);
    }
    
    // Adapted from https://stackoverflow.com/a/4103418
    private static Point2D computeCircleCenter(
        Point2D p0, Point2D p1, Point2D p2)
    {
        double x0 = p0.getX();
        double y0 = p0.getY();
        double x1 = p1.getX();
        double y1 = p1.getY();
        double x2 = p2.getX();
        double y2 = p2.getY();
        double offset = x1 * x1 + y1 * y1;
        double bc = (x0 * x0 + y0 * y0 - offset) / 2.0;
        double cd = (offset - x2 * x2 - y2 * y2) / 2.0;
        double det = (x0 - x1) * (y1 - y2) - (x1 - x2) * (y0 - y1);
        double invDet = 1 / det;
        double cx = (bc * (y1 - y2) - cd * (y0 - y1)) * invDet;
        double cy = (cd * (x0 - x1) - bc * (x1 - x2)) * invDet;
        return new Point2D.Double(cx, cy);
    }
    
    @Override
    protected void paintComponent(Graphics gr)
    {
        super.paintComponent(gr);
        Graphics2D g = (Graphics2D)gr;
        
        g.setColor(Color.RED);
        for (Point2D p : pointList)
        {
            double r = 5;
            g.draw(new Ellipse2D.Double(p.getX()-r, p.getY()-r, r+r, r+r));
        }
        
        g.setColor(Color.BLACK);
        //g.draw(Paths.fromPoints(spline.getInterpolatedPoints(), false));
        
        Point2D p0 = pointList.get(0);
        Point2D p1 = pointList.get(1);
        Point2D p2 = pointList.get(2);
        double curvature = computeCurvature(p0, p1, p2);
        g.drawString("Curvature: "+curvature, 10,  20);
        
        Point2D center = computeCircleCenter(p0, p1, p2);
        double radius = center.distance(p0);
        g.draw(new Ellipse2D.Double(
            center.getX() - radius, center.getY() - radius,
            radius + radius, radius + radius));
    }
    
    @Override
    public void mouseDragged(MouseEvent e)
    {
        if (draggedPoint != null)
        {
            draggedPoint.setLocation(e.getX(), e.getY());
            repaint();
            
            System.out.println("Points: ");
            for (Point2D p : pointList)
            {
                System.out.println("    "+p);
            }
        }
    }


    @Override
    public void mousePressed(MouseEvent e)
    {
        final double thresholdSquared = 10 * 10;
        Point2D p = e.getPoint();
        Point2D closestPoint = null;
        double minDistanceSquared = Double.MAX_VALUE;
        for (Point2D point : pointList)
        {
            double dd = point.distanceSq(p);
            if (dd < thresholdSquared && dd < minDistanceSquared)
            {
                minDistanceSquared = dd;
                closestPoint = point;
            }
        }
        draggedPoint = closestPoint;
    }

    @Override
    public void mouseReleased(MouseEvent e)
    {
        draggedPoint = null;
    }

    @Override
    public void mouseMoved(MouseEvent e)
    {
        // Nothing to do here
    }


    @Override
    public void mouseClicked(MouseEvent e)
    {
        // Nothing to do here
    }

    @Override
    public void mouseEntered(MouseEvent e)
    {
        // Nothing to do here
    }


    @Override
    public void mouseExited(MouseEvent e)
    {
        // Nothing to do here
    }


}
Marco13
  • 53,703
  • 9
  • 80
  • 159
  • Impressive work. @Spen : This man should get the accepted answer, he deserves it! – Eric Duminil Dec 14 '16 at 15:42
  • Minor nitpick : `sideLength0` shouldn't appear three times, shoud it? – Eric Duminil Dec 14 '16 at 15:44
  • Finally, since curvature depends on scale, it could be a good idea to show the unit length on the screen. But again : impressive and fun work! – Eric Duminil Dec 14 '16 at 15:45
  • @EricDuminil Most of my answer is just glitter around what you already wrote, and your answer is in this sense "more to the point" (except for the actual `computeCurvature` function which one would still have to assemble from your snippets - you may add this to your answer if you like, because most people will first look at the accepted answer). Regarding the typos: I'll fix this now. – Marco13 Dec 14 '16 at 15:46
  • computeCurvature appears to be returning twice the curvature since the `area` formula is computing 2x the triangle area. Looks like the `4` should be changed to a `2` in the return statement. – Josh Tauberer Jul 15 '21 at 16:16
  • @JoshTauberer I'll verify that ASAP and update if needed (sounds plausible for now, but will check...) – Marco13 Jul 18 '21 at 00:52
  • `double area = dx1 * dy2 - dy1 * dx2;` returns TWICE the signed area for a triangle. So you need to update this code to `area = 0.5 * (dx1 * dy2 - dy1 * dx2);` – Njål Arne Gjermundshaug Nov 11 '21 at 09:26
  • @NjålArneGjermundshaug and JoshTauberer: Fixed the area computation – Marco13 Nov 12 '21 at 13:38
3

From the wiki you referenced, the curvature is defined as Curvature

where A is the area enclosed by the triangle formed by the three points, x, y and z (p1, p2, p3 in your case) and |x-y| is the distance between points x and y.

Translate the formula to code and you're done!

Hennio
  • 433
  • 2
  • 18
1

C/C++

// https://www.mathopenref.com/coordtrianglearea.html
float getAreaOfTriangle(Point2f A, Point2f B, Point2f C)
{
    return fabs(
            (A.x * (B.y - C.y) + B.x * (C.y - A.y) + C.x * (A.y - B.y)) / 2);
}

float getDistFromPtToPt(Point2f pt1, Point2f pt2)
{
    return sqrt((pt2.x - pt1.x) * (pt2.x - pt1.x) +
                (pt2.y - pt1.y) * (pt2.y - pt1.y));
}


// https://en.wikipedia.org/wiki/Menger_curvature
float
getCurvatureUsingTriangle(Point2f pt1, Point2f pt2, Point2f pt3, bool bDebug)
{
    float fAreaOfTriangle = getAreaOfTriangle(pt1, pt2, pt3);
    float fDist12 = getDistFromPtToPt(pt1, pt2);
    float fDist23 = getDistFromPtToPt(pt2, pt3);
    float fDist13 = getDistFromPtToPt(pt1, pt3);
    float fKappa = 4 * fAreaOfTriangle / (fDist12 * fDist23 * fDist13);
    return fKappa;
}
marikhu
  • 51
  • 5
1

If someone wants this in python:

def getCurvature(x1, y1, x2, y2, x3, y3):
    point1 = (x1, y1)
    point2 = (x2, y2)
    point3 = (x3, y3)

    # Calculating length of all three sides
    len_side_1 = round( math.dist(point1, point2), 2)
    len_side_2 = round( math.dist(point2, point3), 2)
    len_side_3 = round( math.dist(point1, point3), 2)

    # sp is semi-perimeter
    sp = (len_side_1 + len_side_2 + len_side_3) / 2

    # Calculating area using Herons formula
    area = math.sqrt(sp * (sp - len_side_1) * (sp - len_side_2) * (sp - len_side_3))

    # Calculating curvature using Menger curvature formula
    curvature = (4 * area) / (len_side_1 * len_side_2 * len_side_3)

    return curvature
Mohit Kumar
  • 11
  • 2