Calculate curvature for 3 Points (x,y)

Question

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)?

• Curvature: https://en.wikipedia.org/wiki/Curvature
• Menger Curvature: https://en.wikipedia.org/wiki/Menger_curvature

Egidius

please add more information about what you have tried or what is not working in your code – Actualize

Answer 1

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).

Phylogeny

I don't know what Point2D that is, but the Java Point2D doesn't have these methods: docs.oracle.com/javase/7/docs/api/java/awt/geom/Point2D.html – Egidius

So define them, and replace b.x with b.getX() ;) – Phylogeny

I'm about to do so. Thank you :) – Egidius

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. – Appease

I like the approach involving Heron's formula. a bit of code: s = (a + b + c) / 2; area = sqrt(s * (s - a) * (s - b) * (s - c)); return 4 * area / (a * b * c); (where a,b,c are the side lengths) – Allerie

Answer 2

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:

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
    }


}

Rozamond

Impressive work. @Spen : This man should get the accepted answer, he deserves it! – Phylogeny

Minor nitpick : sideLength0 shouldn't appear three times, shoud it? – Phylogeny

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! – Phylogeny

@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. – Rozamond

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. – Moe

@JoshTauberer I'll verify that ASAP and update if needed (sounds plausible for now, but will check...) – Rozamond

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); – Omniumgatherum

@NjålArneGjermundshaug and JoshTauberer: Fixed the area computation – Rozamond

Answer 3

From the wiki you referenced, the curvature is defined as

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!

Actualize

Answer 4

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;
}

Vermont

Answer 5

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

Columbary

Answer 6

Paul Du Bois mentioned Heron's formula for the area. I find this formula to be elegant, so here is a bit of Java that uses it:

import java.awt.geom.Point2D;

public static double menger_curvature(Point2D p1, Point2D p2, Point2D p3)
{
    // side lengths
    double a = p1.distance(p2);
    double b = p2.distance(p3);
    double c = p3.distance(p1);

    // semi-perimeter
    double s = (a + b + c) / 2;

    // Heron's formula for the area of a triangle
    double area = Math.sqrt(s * (s - a) * (s - b) * (s - c));

    double abc = a * b * c; // becomes 0 if points coincide

    return (abc == 0) ? 0 : 4 * area / abc;
}

Allerie