AvgCurvature.java

Go to the documentation of this file.

package theba.descriptors;

import java.awt.Point;
import java.util.ArrayList;

import theba.core.ImageFunctions;
import theba.core.RegionDescriptor;
import theba.core.RegionMask;
import theba.core.math.NonUniformInterpolation;


public class AvgCurvature implements RegionDescriptor {

        public Object measure(RegionMask vmask) {
                short[] mask = new short[vmask.getWidth() * vmask.getHeight()];
                for (int x = 0; x < vmask.getWidth(); x++)
                        for (int y = 0; y < vmask.getHeight(); y++)
                                if (vmask.isSet(x, y, 0))
                                        mask[x + y * vmask.getWidth()] = 1;
                return measure(mask, vmask.getWidth(), vmask.getHeight());
        }

        public double measure(short[] pixels, int width, int height) {
                        
                final double[] d1xGaussian = { 0.0269, 0.0237, 0.0140, 0, -0.0140,
                                -0.0237, -0.0269 };
                final double[] d2xGaussian = { 0, -0.0066, -0.0124, -0.0148, -0.0124,
                                -0.0066, 0 };

                final int radius = 3;

                ArrayList points = ImageFunctions.borderTrace(pixels, width, height);
                int[] chaincode = ImageFunctions.getChaincode(points);
                double[] s = ImageFunctions.getArcLenghts(chaincode);
                int vectorLength = points.size() + 2 * radius;

                double x[] = new double[points.size()];
                double y[] = new double[points.size()];
                Point p = null;
                for (int i = 0; i < points.size(); i++) {
                        p = (Point) points.get(i);
                        x[i] = p.x;
                        y[i] = p.y;
                }

                NonUniformInterpolation nx = new NonUniformInterpolation(s, x);
                NonUniformInterpolation ny = new NonUniformInterpolation(s, y);
                double step = s[s.length - 1] / s.length; // average distance between
                                                                                                        // points

                double xi[] = new double[vectorLength];
                double yi[] = new double[vectorLength];
                double xMin = Double.MAX_VALUE;
                double yMin = Double.MAX_VALUE;
                double xMax = Double.MIN_VALUE;
                double yMax = Double.MIN_VALUE;
                double arcLength = 0;
                for (int i = 0; i < s.length; i++) {
                        xi[i + radius] = nx.interpolate(arcLength);
                        yi[i + radius] = ny.interpolate(arcLength);

                        if (xi[i + radius] < xMin)
                                xMin = xi[i + radius];
                        if (yi[i + radius] < yMin)
                                yMin = yi[i + radius];

                        if (xi[i + radius] > xMax)
                                xMax = xi[i + radius];
                        if (yi[i + radius] > yMax)
                                yMax = yi[i + radius];

                        arcLength += step;
                }

                // for simplicity:
                x = xi;
                y = yi;

                // pad with real values to avoid hit by zeros
                // get real values by wrap-around
                for (int k = 0; k < radius; k++) {
                        x[k] = x[x.length - 2 * radius + k];
                        y[k] = y[x.length - 2 * radius + k];
                }
                for (int k = x.length - radius; k < x.length; k++) {
                        x[k] = x[k - x.length + 2 * radius];
                        y[k] = y[k - x.length + 2 * radius];
                }

                // Make translation invariant
                for (int k = 0; k < x.length; k++) {
                        x[k] -= xMin;
                        y[k] -= yMin;
                }

                double[] dx = new double[vectorLength];
                double[] dy = new double[vectorLength];
                double[] dxx = new double[vectorLength];
                double[] dyy = new double[vectorLength];
                // Calculate df/dx, df/dy, df/dxx, df/dyy by convolution
                for (int k = radius; k < vectorLength - radius; k++) {

                        for (int m = -radius; m <= radius; m++) {
                                dx[k] += x[k + m] * d1xGaussian[m + radius];
                                dxx[k] += x[k + m] * d2xGaussian[m + radius];

                                dy[k] += y[k + m] * d1xGaussian[m + radius];
                                dyy[k] += y[k + m] * d2xGaussian[m + radius];
                        }
                }

                double curvatureSum = 0;
                for (int k = radius; k < vectorLength - radius; k++) {
                        double factor = dx[k] * dyy[k] - dy[k] * dxx[k];
                        double dividend = Math.pow(dx[k] * dx[k] + dy[k] * dy[k], 1.5);
                        // guard against small dividend, its always positive
                        if (dividend < 0.001) {
                                dividend = 0.001;
                        }
                        curvatureSum += Math.abs(factor / dividend);

                }
                // Normalize results by dividing by number of points of border and
                // return
                return curvatureSum / (points.size() * points.size());

        }

        public String getName() {
                return "Average, scale-independent curvature";
        }

        public String getAbout() {
                return "Returns the average, weighted curvature of a 4-connected 2D region";
        }

        public boolean does3D() {
                return false;
        }

        public boolean isNumeric() {
                return true;
        }

}