boofcv.alg.shapes.ellipse.SnapToEllipseEdge Maven / Gradle / Ivy

/*
 * Copyright (c) 2011-2016, Peter Abeles. All Rights Reserved.
 *
 * This file is part of BoofCV (http://boofcv.org).
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package boofcv.alg.shapes.ellipse;

import boofcv.alg.shapes.edge.BaseIntegralEdge;
import boofcv.struct.image.ImageGray;
import georegression.fitting.ellipse.FitEllipseWeightedAlgebraic;
import georegression.geometry.UtilEllipse_F64;
import georegression.metric.UtilAngle;
import georegression.struct.point.Point2D_F64;
import georegression.struct.shapes.EllipseRotated_F64;
import org.ddogleg.struct.FastQueue;
import org.ddogleg.struct.GrowQueue_F64;

/**
 * Refines an initial estimate of an elipse using a subpixel contour technique.  A local line integral around each
 * point is used to determine how important each point is.  The contour being
 *
 * @author Peter Abeles
 */
public class SnapToEllipseEdge extends BaseIntegralEdge {

	// maximum number of iterations it will performance
	protected int maxIterations = 10;
	// when the difference between two ellipses is less than this amount stop iterating
	protected double convergenceTol = 1e-6;

	// how many points along the contour it will sample
	protected int numSampleContour;

	// Determines the number of points sampled radially outwards from the line
	// Total intensity values sampled at each point along the line is radius*2+2,
	// and points added to line fitting is radius*2+1.
	protected int radialSamples;

	protected GrowQueue_F64 weights = new GrowQueue_F64();// storage for weights in line fitting
	// storage for where the points that are sampled along the line
	protected FastQueue samplePts = new FastQueue<>(Point2D_F64.class, true);

	protected FitEllipseWeightedAlgebraic fitter = new FitEllipseWeightedAlgebraic();

	protected EllipseRotated_F64 previous = new EllipseRotated_F64();

	/**
	 * Constructor with configuration
	 *
	 * @param numSampleContour Maximum number of iterations it will performance
	 * @param radialSamples When the difference between two ellipses is less than this amount stop iterating
	 * @param imageType Type of gray-scale input image
	 */
	public SnapToEllipseEdge( int numSampleContour, int radialSamples , Class imageType) {
		super(imageType);

		this.numSampleContour = numSampleContour;
		this.radialSamples = radialSamples;
	}

	/**
	 * Refines provided list by snapping it to edges found in the image
	 *
	 * @param input (Output) Close approximation of the ellipse in the image
	 * @param refined (Output) Storage for refined estimate.  Can be same instance as input
	 * @return True if a refined estimate could be found, false if it failed
	 */
	public boolean process(EllipseRotated_F64 input, EllipseRotated_F64 refined) {

		refined.set(input);
		previous.set(input);

		for (int iteration = 0; iteration < maxIterations; iteration++) {
			refined.set(previous);
			computePointsAndWeights(refined);

			if( fitter.process(samplePts.toList(),weights.data) ) {
				// Get the results in local coordinates
				UtilEllipse_F64.convert(fitter.getEllipse(),refined);
				// convert back into image coordinates
				double scale = previous.a;
				refined.center.x = refined.center.x*scale + previous.center.x;
				refined.center.y = refined.center.y*scale + previous.center.y;
				refined.a *= scale;
				refined.b *= scale;
			} else {
				return false;
			}

			// stop once the change between two iterations is insignificant
			if( change(previous,refined) <= convergenceTol) {
				return true;
			} else {
				previous.set(refined);
			}
		}
		return true;
	}

	/**
	 * Computes a numerical value for the difference in parameters between the two ellipses
	 */
	protected static double change( EllipseRotated_F64 a , EllipseRotated_F64 b ) {
		double total = 0;

		total += Math.abs(a.center.x - b.center.x);
		total += Math.abs(a.center.y - b.center.y);
		total += Math.abs(a.a - b.a);
		total += Math.abs(a.b - b.b);

		// only care about the change of angle when it is not a circle
		double weight = Math.min(4,2.0*(a.a/a.b-1.0));
		total += weight*UtilAngle.distHalf(a.phi , b.phi);

		return total;
	}

	/**
	 * Computes the location of points along the line and their weights
	 */
	void computePointsAndWeights(EllipseRotated_F64 ellipse) {

		// use the semi-major axis to scale the input points for numerical stability
		double localScale = ellipse.a;

		samplePts.reset();
		weights.reset();
		int numSamples = radialSamples * 2 + 2;
		int numPts = numSamples - 1;

		Point2D_F64 sample = new Point2D_F64();
		for (int i = 0; i < numSampleContour; i++) {

			// find a point along the ellipse at evenly spaced angles
			double theta = 2.0 * Math.PI * i / numSampleContour;

			UtilEllipse_F64.computePoint(theta, ellipse, sample);

			// compute the unit tangent along the ellipse at this point
			double tanX = sample.x - ellipse.center.x;
			double tanY = sample.y - ellipse.center.y;
			double r = Math.sqrt(tanX * tanX + tanY * tanY);
			tanX /= r;
			tanY /= r;

			// define the line it will sample along
			double x = sample.x - numSamples * tanX / 2.0;
			double y = sample.y - numSamples * tanY / 2.0;

			double lengthX = numSamples * tanX;
			double lengthY = numSamples * tanY;

			// Unless all the sample points are inside the image, ignore this point
			if (!integral.isInside(x, y) || !integral.isInside(x + lengthX, y + lengthY))
				continue;

			double sample0 = integral.compute(x, y, x + tanX, y + tanY);
			x += tanX;
			y += tanY;
			for (int j = 0; j < numPts; j++) {
				double sample1 = integral.compute(x, y, x + tanX, y + tanY);

				double w = sample0 - sample1;
				if (w < 0) w = -w;

				if (w > 0) {
					// convert into a local coordinate so make the linear fitting more numerically stable and
					// independent on position in the image
					samplePts.grow().set((x - ellipse.center.x) / localScale, (y - ellipse.center.y) / localScale);
					weights.add(w);
				}

				x += tanX;
				y += tanY;
				sample0 = sample1;
			}

		}
	}

	public int getMaxIterations() {
		return maxIterations;
	}

	public void setMaxIterations(int maxIterations) {
		this.maxIterations = maxIterations;
	}

	public double getConvergenceTol() {
		return convergenceTol;
	}

	public void setConvergenceTol(double convergenceTol) {
		this.convergenceTol = convergenceTol;
	}
}