BendingEnergyDescriptor.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;


public class BendingEnergyDescriptor implements RegionDescriptor {

        public double measure(short[] pixels, int width, int height) {

                // Gaussian 1d-kernels for sigma = 3.0:
                // standard gaussian
                // final double[] d0xGaussian = {0.0807 , 0.1065 , 0.1258 , 0.1330,
                // 0.1258, 0.1065, 0.0807};
                // derivated 1 time
                final double[] d1xGaussian = { 0.0269, 0.0237, 0.0140, 0, -0.0140,
                                -0.0237, -0.0269 };
                // derivated 2 times
                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);
                if (points == null)
                        return -1; // return impossible value as error indication

                int vectorLength = points.size() + 2 * radius;
                double x[] = new double[vectorLength];
                double y[] = new double[vectorLength];

                Point p = (Point) points.get(0);
                double xMin = p.x, yMin = p.y;
                for (int k = 0; k < points.size(); k++) {
                        p = (Point) points.get(k);
                        x[k + radius] = p.x;
                        y[k + radius] = p.y;

                        if (p.x < xMin)
                                xMin = p.x;
                        if (p.y < yMin)
                                yMin = p.y;

                }

                // 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 bendingEnergySum = 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, sign doesn't matter
                        if (Math.abs(dividend) < 0.001) {
                                dividend = 0.001;
                        }
                        bendingEnergySum += (factor / dividend) * (factor / dividend);

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

        public String getName() {
                return "Bending energy";
        }

        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 String getAbout() {
                // TODO Auto-generated method stub
                return null;
        }

        public boolean does3D() {
                // TODO Auto-generated method stub
                return false;
        }

        public boolean isNumeric() {
                return true;
        }

}