Find largest convex area within a given shape

MX455687a.java

/*
 * Proof of concept for
 * http://stackoverflow.com/a/17967201/1468366
 * http://math.stackexchange.com/a/456233/35416
 *
 * Note: This code was written ONLY to get this single example image.
 * It makes a number of implicit assumptions, and therefore won't be
 * fit for general applications without some modifications.
 * Furthermore, speedy coding was considered more important than
 * performance at run time.
 *
 * 2013-07-31 Martin von Gagern
 */

import java.awt.*;
import java.awt.geom.*;
import javax.swing.*;
import java.util.*;
import java.util.List;

class MX455687a extends JPanel {

    Path2D.Double p, fp, convexParts, concaveParts;
    List<Point2D.Double> pts;
    List<Integer> inflectionPoints;

    MX455687a() {
        setPreferredSize(new Dimension(600, 600));
        initPath();
        flatten();
        findInflections();
        composeParts();
        opt();
    }

    void initPath() {
        double x, y;
        p = new Path2D.Double();

        // These coordinates come from Inkscape.

        p.moveTo(x = 237.68304, y = 996.48825-452.36218);
        p.curveTo(x + 6.51879, y - 1.8326,
                  x + 21.80476, y - 7.8128,
                  x = x + 33.96883, y = y - 13.2891);
        p.curveTo(x + 35.39892, y - 15.9369,
                  x + 57.77364, y - 20.8704,
                  x = x + 95.533, y = y - 21.0644);
        p.curveTo(x + 17.94058, y - 0.092,
                  x + 46.49802, y + 1.5691,
                  x = x + 65.58485, y = y + 3.8151);
        p.curveTo(x + 42.7117, y + 5.0265,
                  x + 67.20049, y + 3.9623,
                  x = x + 83.80193, y = y - 3.6408);
        p.curveTo(x + 16.30881, y - 7.4692,
                  x + 26.86538, y - 17.33109,
                  x = x + 33.85875, y = y - 31.63035);
        p.curveTo(x + 8.9218, y - 18.24244,
                  x + 9.37729, y - 36.65979,
                  x = x + 1.83938, y = y - 74.37299);
        p.curveTo(x - 4.18125, y - 20.91934,
                  x - 6.16273, y - 37.79353,
                  x = x - 6.28951, y = y - 53.56106);
        p.curveTo(x - 0.16743, y - 20.82869,
                  x + 0.39748, y - 24.26443,
                  x = x + 5.93532, y = y - 36.09556);
        p.curveTo(x + 3.66125, y - 7.82194,
                  x + 11.82786, y - 18.98987,
                  x = x + 20.32519, y = y - 27.79488);
        p.curveTo(x + 18.95025, y - 19.63646,
                  x + 23.84825, y - 30.72434,
                  x = x + 23.58954, y = y - 53.40104);

        p.curveTo(595.26225, 635.6563-452.36218,
                  553.20055, 577.42-452.36218,
                  x = 499.75728, y = 552.43576-452.36218);

        p.curveTo(x - 16.60094, y - 7.7608,
                  x - 23.50652, y - 8.25737,
                  x = x - 83.26282, y = y - 5.98713);
        p.curveTo(x - 28.66827, y + 1.08916,
                  x - 44.1929, y - 3.69973,
                  x = x - 77.79958, y = y - 23.99885);
        p.curveTo(x - 32.27854, y - 19.49689,
                  x - 40.26394, y - 22.01741,
                  x = x - 53.60573, y = y - 16.9202);
        p.curveTo(x - 19.38047, y + 7.40425,
                  x - 23.91467, y + 17.04982,
                  x = x - 34.82266, y = y + 74.07759);
        p.curveTo(x - 6.2155, y + 32.49513,
                  x - 13.92691, y + 50.61982,
                  x = x - 28.24379, y = y + 66.38345);
        p.curveTo(x - 19.25485, y + 21.2006,
                  x - 38.95949, y + 30.19284,
                  x = x - 110.92011, y = y + 50.61855);
        p.curveTo(x - 27.286991, y + 7.74529,
                  x - 56.347896, y + 16.65699,
                  x = x - 64.579815, y = y + 19.80378);

        p.curveTo(16.76242, 727.7893-452.36218,
                  1.8467711, 745.44327-452.36218,
                  4.4516545, 766.20783-452.36218);
        p.curveTo(7.6378213, 791.60621-452.36218,
                  16.690734, 807.26349-452.36218,
                  74.13503, 886.72752-452.36218);
        p.curveTo(150.96578, 993.00915-452.36218,
                  180.22006, 1012.6424-452.36218,
                  237.68304, 996.48825-452.36218);

        p.closePath();
    }

    void flatten() {
        fp = new Path2D.Double();
        fp.append(p.getPathIterator(null, 0.1), false);
        pts = new ArrayList<Point2D.Double>();
        double[] coords = new double[6];
        for (PathIterator pi = fp.getPathIterator(null);
             !pi.isDone(); pi.next()) {
            if (pi.currentSegment(coords) != PathIterator.SEG_CLOSE)
                pts.add(new Point2D.Double(coords[0], coords[1]));
        }
    }

    void findInflections() {
        inflectionPoints = new ArrayList<Integer>();
        Point2D.Double a = pts.get(pts.size() - 2);
        Point2D.Double b = pts.get(0), c;
        double oldDet = -1, det;
        // negative det means convex due to the orientation of the path
        // the initial vertex has negative det due to the starting point
        for (int i = 1; i != pts.size(); ++i) {
            c = pts.get(i);
            det = a.x*b.y + b.x*c.y + c.x*a.y - a.y*b.x - b.y*c.x - c.y*a.x;
            if (det * oldDet < 0) {
                System.out.println("At " + i + "/" + pts.size() +
                                   " from " + oldDet + " to " + det);
                inflectionPoints.add(det < 0 ? i - 1 : i);
            }
            a = b;
            b = c;
            oldDet = det;
        }
    }

    void composeParts() {
        convexParts = new Path2D.Double();
        concaveParts = new Path2D.Double();
        Point2D.Double p;
        int begin, end;
        for (int i = 0; i < inflectionPoints.size(); i += 2) {
            begin = inflectionPoints.get(i);
            end = inflectionPoints.get(i + 1);
            p = pts.get(begin);
            concaveParts.moveTo(p.x, p.y);
            for (int j = begin + 1; j <= end; ++j) {
                p = pts.get(j);
                concaveParts.lineTo(p.x, p.y);
            }
        }
        for (int i = 1; i < inflectionPoints.size() - 1; i += 2) {
            begin = inflectionPoints.get(i);
            end = inflectionPoints.get(i + 1);
            p = pts.get(begin);
            convexParts.moveTo(p.x, p.y);
            for (int j = begin + 1; j <= end; ++j) {
                p = pts.get(j);
                convexParts.lineTo(p.x, p.y);
            }
        }
        begin = inflectionPoints.get(inflectionPoints.size() - 1);
        end = pts.size();
        p = pts.get(begin);
        convexParts.moveTo(p.x, p.y);
        for (int j = begin + 1; j < end; ++j) {
            p = pts.get(j);
            convexParts.lineTo(p.x, p.y);
        }
        end = inflectionPoints.get(0);
        for (int j = 0; j <= end; ++j) {
            p = pts.get(j);
            convexParts.lineTo(p.x, p.y);
        }
    }

    double bestValue;

    Shape bestShape;

    List<Line2D> lines, bestLines;

    void opt() {
        bestValue = 0;
        lines = new ArrayList<Line2D>();
        opt(0, new Area(fp));
    }

    void opt(int i, Area a1) {
        if (i == inflectionPoints.size()) {
            double a = area(a1);
            if (a > bestValue) {
                bestValue = a;
                bestShape = a1;
                bestLines = new ArrayList<Line2D>(lines);
            }
            return;
        }
        int begin = inflectionPoints.get(i), end = inflectionPoints.get(i + 1);
        Point2D.Double a = pts.get(begin - 1), b;
        for (int j = begin; j <= end + 1; ++j) {
            b = pts.get(j);
            double f = 1200/a.distance(b);
            double dx = f*(b.x - a.x), dy = f*(b.y - a.y);
            Path2D.Double p = new Path2D.Double();
            p.moveTo(a.x - dx     , a.y - dy     );
            p.lineTo(a.x - dx + dy, a.y - dy - dx);
            p.lineTo(a.x + dx + dy, a.y + dy - dx);
            p.lineTo(a.x + dx     , a.y + dy     );
            p.closePath();
            Area a2 = new Area(a1);
            a2.intersect(new Area(p));
            lines.add(new Line2D.Double(a.x - dx, a.y - dy,
                                        a.x + dx, a.y + dy));
            opt(i + 2, a2);
            lines.remove(lines.size() - 1);
            a = b;
        }
    }

    double area(Shape s) {
        double total = 0, part = 0;
        double[] coords = new double[6];
        Point2D.Double initial = null;
        double ax = 0, ay = 0, bx, by;
        for (PathIterator pi = s.getPathIterator(null);
             !pi.isDone(); pi.next()) {
            switch (pi.currentSegment(coords)) {
            case PathIterator.SEG_MOVETO:
                ax = coords[0];
                ay = coords[1];
                initial = new Point2D.Double(ax, ay);
                part = 0;
                break;
            case PathIterator.SEG_CLOSE:
                bx = initial.x;
                by = initial.y;
                part += ax*by - bx*ay;
                total += part;
                part = 0;
                ax = bx;
                ay = by;
                break;
            case PathIterator.SEG_LINETO:
                bx = coords[0];
                by = coords[1];
                part += ax*by - bx*ay;
                ax = bx;
                ay = by;
                break;
            default:
                throw new IllegalArgumentException("Non-flat shape");
            }
        }
        if (total < 0)
            total = -total;
        return total;
    }

    @Override public void paintComponent(Graphics g) {
        Graphics2D g2d = (Graphics2D)g;
        g2d.setRenderingHint(RenderingHints.KEY_ANTIALIASING,
                             RenderingHints.VALUE_ANTIALIAS_ON);
        g2d.setRenderingHint(RenderingHints.KEY_STROKE_CONTROL,
                             RenderingHints.VALUE_STROKE_PURE);
        g2d.setColor(Color.WHITE);
        g2d.fillRect(0, 0, getWidth(), getHeight());
        g2d.setColor(new Color(128, 255, 128));
        g2d.fill(bestShape);
        g2d.setColor(Color.BLUE);
        g2d.setStroke(new BasicStroke(1.25f));
        for (Line2D l: bestLines)
            g2d.draw(l);
        g2d.setStroke(new BasicStroke(2.f));
        g2d.setColor(Color.BLACK);
        g2d.draw(convexParts);
        g2d.setColor(Color.RED);
        g2d.draw(concaveParts);
    }

    public static void main(String[] args) {
        JFrame f = new JFrame();
        f.getContentPane().add(new MX455687a());
        f.pack();
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setVisible(true);
    }

}