com.stromberglabs.jopensurf.Surf.java Source code

Java tutorial

Introduction

Here is the source code for com.stromberglabs.jopensurf.Surf.java

Source

/*
 * 
 * See also "LICENCE" file in this package.
 * 
This work was derived from Chris Evan's opensurf project and re-licensed as the
3 clause BSD license with permission of the original author. Thank you Chris! 
    
Copyright (c) 2010, Andrew Stromberg
All rights reserved.
    
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
 * Neither Andrew Stromberg nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
    
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL Andrew Stromberg BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
package com.stromberglabs.jopensurf;

import com.heeere.androjsurf.IntensityProvider;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;
import java.util.Map;
import java.util.HashMap;
//import static org.apache.commons.math.util.FastMath.*;
import static java.lang.Math.*;

/**
 * A class to calculate the upright or free oriented interest points of an image lazily (will not calculate until you ask for them)
 * @author astromberg
 *
 */
public class Surf implements Serializable {

    private static final long serialVersionUID = 1L;
    private static final int HESSIAN_OCTAVES = 5;
    //private static final int HESSIAN_INTERVALS = 4;
    private static final int HESSIAN_INIT_SAMPLE = 2;
    private static final float HESSIAN_THRESHOLD = 0.0004F;
    private static final float HESSIAN_BALANCE_VALUE = 0.81F;
    private FastHessian mHessian;
    private List<SURFInterestPoint> mFreeOrientedPoints;
    private List<SURFInterestPoint> mUprightPoints;
    private List<SURFInterestPoint> mDescriptorFreeInterestPoints;
    //octaves = steps
    private int mNumOctaves = HESSIAN_OCTAVES;
    //private int mNumIntervals = HESSIAN_INTERVALS;
    private float mThreshold = HESSIAN_THRESHOLD;
    private float mBalanceValue = HESSIAN_BALANCE_VALUE;
    private IntegralImage mIntegralImage;

    public Surf(IntensityProvider image) {
        this(image, HESSIAN_BALANCE_VALUE, HESSIAN_THRESHOLD, HESSIAN_OCTAVES);
    }

    public Surf(IntensityProvider image, float balanceValue, float threshold, int octaves) {
        mNumOctaves = octaves;
        mBalanceValue = balanceValue;
        mThreshold = threshold;

        //Calculate the integral image
        mIntegralImage = new IntegralImage(image);

        //Calculate the fast hessian
        mHessian = new FastHessian(mIntegralImage, mNumOctaves, HESSIAN_INIT_SAMPLE, mThreshold, mBalanceValue);

        //Calculate the descriptor and orientation free interest points
        mDescriptorFreeInterestPoints = mHessian.getIPoints();
    }

    public List<SURFInterestPoint> getUprightInterestPoints() {
        return getPoints(true);
    }

    public List<SURFInterestPoint> getFreeOrientedInterestPoints() {
        return getPoints(false);
    }

    private List<SURFInterestPoint> getPoints(boolean upright) {
        List<SURFInterestPoint> points = upright ? mUprightPoints : mFreeOrientedPoints;
        if (points == null) {
            points = getDescriptorFreeInterestPoints();
            // cache for next time through
            if (upright) {
                mUprightPoints = points;
            } else {
                mFreeOrientedPoints = points;
            }
            for (SURFInterestPoint point : points) {
                getOrientation(point);
                getMDescriptor(point, upright);
            }
        }
        return points;
    }

    private List<SURFInterestPoint> getDescriptorFreeInterestPoints() {
        List<SURFInterestPoint> points = new ArrayList<SURFInterestPoint>(mDescriptorFreeInterestPoints.size());
        for (SURFInterestPoint point : mDescriptorFreeInterestPoints) {
            try {
                points.add((SURFInterestPoint) point.clone());
            } catch (CloneNotSupportedException e) {
                e.printStackTrace();
            }
        }
        return points;
    }

    //TODO: Why the multiply by 1?
    private float haarX(int row, int column, int s) {
        return ImageTransformUtils.BoxIntegral(mIntegralImage, row - s / 2, column, s, s / 2)
                - 1 * ImageTransformUtils.BoxIntegral(mIntegralImage, row - s / 2, column - s / 2, s, s / 2);
    }

    //TODO: Why the multiply by 1?
    private float haarY(int row, int column, int s) {
        return ImageTransformUtils.BoxIntegral(mIntegralImage, row, column - s / 2, s / 2, s)
                - 1 * ImageTransformUtils.BoxIntegral(mIntegralImage, row - s / 2, column - s / 2, s / 2, s);
    }

    private void getOrientation(SURFInterestPoint input) {
        double gauss;
        float scale = input.getScale();

        int s = (int) round(scale);
        int r = (int) round(input.getY());
        int c = (int) round(input.getX());

        List<Double> xHaarResponses = new ArrayList<Double>();
        List<Double> yHaarResponses = new ArrayList<Double>();
        List<Double> angles = new ArrayList<Double>();
        //System.out.println("s = " + s + ", r = " + r + ", c = " + c + ", scale = " + scale);

        //calculate haar responses for points within radius of 6*scale
        for (int i = -6; i <= 6; ++i) {
            for (int j = -6; j <= 6; ++j) {
                if (i * i + j * j < 36) {
                    gauss = GaussianConstants.Gauss25[abs(i)][abs(j)];
                    //System.out.println("i = " + i + ", j = " + j + ", gauss = " + gauss);
                    double xHaarResponse = gauss * haarX(r + j * s, c + i * s, 4 * s);
                    double yHaarResponse = gauss * haarY(r + j * s, c + i * s, 4 * s);
                    xHaarResponses.add(xHaarResponse);
                    yHaarResponses.add(yHaarResponse);
                    angles.add(getAngle(xHaarResponse, yHaarResponse));
                    //System.out.format("gauss = %.8f, haarX = %.8f, haarY = %.8f, getAngle(x,y) = %.8f",gauss,xHaarResponse,yHaarResponse,getAngle(xHaarResponse,yHaarResponse));
                    //System.out.println();
                }
            }
        }

        // calculate the dominant direction
        float sumX = 0, sumY = 0;
        float ang1, ang2, ang;
        float max = 0;
        float orientation = 0;

        // loop slides pi/3 window around feature point
        for (ang1 = 0; ang1 < 2 * PI; ang1 += 0.15f) {
            ang2 = (float) (ang1 + PI / 3.0f > 2 * Math.PI ? ang1 - 5.0f * PI / 3.0f : ang1 + PI / 3.0f);
            sumX = sumY = 0;
            for (int k = 0; k < angles.size(); k++) {
                ang = angles.get(k).floatValue();

                if (ang1 < ang2 && ang1 < ang && ang < ang2) {
                    sumX += xHaarResponses.get(k).floatValue();
                    sumY += yHaarResponses.get(k).floatValue();
                } else if (ang2 < ang1 && ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2 * PI))) {
                    sumX += xHaarResponses.get(k).floatValue();
                    sumY += yHaarResponses.get(k).floatValue();
                }
            }
            // if the vector produced from this window is longer than all
            // previous vectors then this forms the new dominant direction
            if (sumX * sumX + sumY * sumY > max) {
                // store largest orientation
                max = sumX * sumX + sumY * sumY;
                orientation = (float) getAngle(sumX, sumY);
            }
        }
        input.setOrientation(orientation);
        //System.out.println("orientation = " + orientation);
        //return input;
    }

    private void getMDescriptor(SURFInterestPoint point, boolean upright) {
        int y, x, count = 0;
        int sample_x, sample_y;
        double scale, dx, dy, mdx, mdy, co = 1F, si = 0F;
        float desc[] = new float[64];
        double gauss_s1 = 0.0D, gauss_s2 = 0.0D, xs = 0.0D, ys = 0.0D;
        double rx = 0.0D, ry = 0.0D, rrx = 0.0D, rry = 0.0D, len = 0.0D;
        int i = 0, ix = 0, j = 0, jx = 0;

        float cx = -0.5f, cy = 0.0f; //Subregion centers for the 4x4 gaussian weighting

        scale = point.getScale();
        x = round(point.getX());
        y = round(point.getY());
        //System.out.println("x = " + point.getX() + ", y = " + point.getY());
        //System.out.println("x = " + x + ", y = " + y);
        if (!upright) {
            co = cos(point.getOrientation());
            si = sin(point.getOrientation());
        }
        //System.out.println("co = " + co + ", sin = " + si);
        i = -8;
        //Calculate descriptor for this interest point
        //Area of size 24 s x 24 s
        //***********************************************
        while (i < 12) {
            j = -8;
            i = i - 4;

            cx += 1.0F;
            cy = -0.5F;

            while (j < 12) {
                dx = dy = mdx = mdy = 0.0F;
                cy += 1.0F;

                j = j - 4;

                ix = i + 5;
                jx = j + 5;

                xs = round(x + (-jx * scale * si + ix * scale * co));
                ys = round(y + (jx * scale * co + ix * scale * si));

                for (int k = i; k < i + 9; ++k) {
                    for (int l = j; l < j + 9; ++l) {
                        //Get coords of sample point on the rotated axis
                        sample_x = (int) round(x + (-1D * l * scale * si + k * scale * co));
                        sample_y = (int) round(y + (l * scale * co + k * scale * si));

                        //Get the gaussian weighted x and y responses
                        gauss_s1 = gaussian(xs - sample_x, ys - sample_y, 2.5F * scale);

                        rx = haarX(sample_y, sample_x, (int) (2 * round(scale)));
                        ry = haarY(sample_y, sample_x, (int) (2 * round(scale)));

                        //Get the gaussian weighted x and y responses on rotated axis
                        rrx = gauss_s1 * (-rx * si + ry * co);
                        rry = gauss_s1 * (rx * co + ry * si);

                        dx += rrx;
                        dy += rry;

                        mdx += abs(rrx);
                        mdy += abs(rry);
                    }
                }

                //Add the values to the descriptor vector
                gauss_s2 = gaussian(cx - 2.0f, cy - 2.0f, 1.5f);

                //Casting from a double to a float, might be a terrible idea
                //but doubles are expensive
                desc[count++] = (float) (dx * gauss_s2);
                desc[count++] = (float) (dy * gauss_s2);

                desc[count++] = (float) (mdx * gauss_s2);
                desc[count++] = (float) (mdy * gauss_s2);

                //Accumulate length for vector normalisation
                len += (dx * dx + dy * dy + mdx * mdx + mdy * mdy) * (gauss_s2 * gauss_s2);

                j += 9;
            }
            i += 9;
        }

        len = sqrt(len);

        for (i = 0; i < 64; i++) {
            desc[i] /= len;
        }

        point.setDescriptor(desc);
        //        for ( double v : desc ){
        //           System.out.printf("%.7f",v);
        //           System.out.print(",");
        //        }
        //        System.out.println();
    }

    private double getAngle(double xHaarResponse, double yHaarResponse) {
        if (xHaarResponse >= 0 && yHaarResponse >= 0) {
            return atan(yHaarResponse / xHaarResponse);
        }

        if (xHaarResponse < 0 && yHaarResponse >= 0) {
            return PI - atan(-yHaarResponse / xHaarResponse);
        }

        if (xHaarResponse < 0 && yHaarResponse < 0) {
            return PI + atan(yHaarResponse / xHaarResponse);
        }

        if (xHaarResponse >= 0 && yHaarResponse < 0) {
            return 2 * PI - atan(-yHaarResponse / xHaarResponse);
        }

        return 0;
    }

    public static Map<SURFInterestPoint, SURFInterestPoint> getMatchingPoints(List<SURFInterestPoint> set1,
            List<SURFInterestPoint> set2, double smallestOverSecondSmallestThreshold) {
        //System.out.println("Finding matching points..");
        Map<SURFInterestPoint, SURFInterestPoint> matchingPoints = new HashMap<SURFInterestPoint, SURFInterestPoint>();
        List<SURFInterestPoint> points = set1;

        for (SURFInterestPoint a : points) {
            double smallestDistance = Float.MAX_VALUE;
            double nextSmallestDistance = Float.MAX_VALUE;
            SURFInterestPoint possibleMatch = null;

            for (SURFInterestPoint b : set2) {
                double distance = a.getDistance(b);
                //System.out.println("Distance = " + distance);
                if (distance < smallestDistance) {
                    nextSmallestDistance = smallestDistance;
                    smallestDistance = distance;
                    possibleMatch = b;
                } else if (distance < nextSmallestDistance) {
                    nextSmallestDistance = distance;
                }
            }

            // If match has a d1:d2 ratio < 0.65 ipoints are a match
            //if ( smallestDistance/nextSmallestDistance < 0.65d ){
            if (smallestDistance / nextSmallestDistance < smallestOverSecondSmallestThreshold) {
                //not storing change in position
                matchingPoints.put(a, possibleMatch);
            }
        }
        //System.out.println("Done finding matching points");
        return matchingPoints;
    }

    public Map<SURFInterestPoint, SURFInterestPoint> getMatchingPoints(Surf descriptor, boolean upright) {
        //System.out.println("Finding matching points..");
        Map<SURFInterestPoint, SURFInterestPoint> matchingPoints = new HashMap<SURFInterestPoint, SURFInterestPoint>();
        List<SURFInterestPoint> points = upright ? getUprightInterestPoints() : getFreeOrientedInterestPoints();

        for (SURFInterestPoint a : points) {
            double smallestDistance = Float.MAX_VALUE;
            double nextSmallestDistance = Float.MAX_VALUE;
            SURFInterestPoint possibleMatch = null;

            for (SURFInterestPoint b : upright ? descriptor.getUprightInterestPoints()
                    : descriptor.getFreeOrientedInterestPoints()) {
                double distance = a.getDistance(b);
                //System.out.println("Distance = " + distance);
                if (distance < smallestDistance) {
                    nextSmallestDistance = smallestDistance;
                    smallestDistance = distance;
                    possibleMatch = b;
                } else if (distance < nextSmallestDistance) {
                    nextSmallestDistance = distance;
                }
            }

            // If match has a d1:d2 ratio < 0.65 ipoints are a match
            //if ( smallestDistance/nextSmallestDistance < 0.65d ){
            if (smallestDistance / nextSmallestDistance < 0.75d) {
                //not storing change in position
                matchingPoints.put(a, possibleMatch);
            }
        }
        //System.out.println("Done finding matching points");
        return matchingPoints;
    }

    public String getStringRepresentation(boolean freeOriented) {
        StringBuffer buffer = new StringBuffer();
        for (SURFInterestPoint point : freeOriented ? getFreeOrientedInterestPoints()
                : getUprightInterestPoints()) {
            for (double val : point.getDescriptor()) {
                buffer.append(Double.doubleToLongBits(val) + ",");
            }
        }
        buffer.substring(0, buffer.length() - 1);
        return buffer.toString();
    }

    public void setStringRepresentation(String str) {
    }

    //   private double gaussian(int x, int y, double sig){
    //      return (1.0f/(2.0f*PI*sig*sig)) * exp( -(x*x+y*y)/(2.0f*sig*sig));
    //   }
    private double gaussian(double x, double y, double sig) {
        return (1.0f / (2.0f * PI * sig * sig)) * exp(-(x * x + y * y) / (2.0f * sig * sig));
    }

    //   public static void main(String[] args){
    //      try {
    //         BufferedImage image = ImageIO.read(new File("/home/astromberg/test_images/cover.jpg"));
    //         Surf board = new Surf(image);
    //         List<InterestPoint> matchingPoints = board.getMatchingPoints(board,false);
    //         System.out.println("between the first image and itself there are " + matchingPoints.size() + " of " + board.getFreeOrientedInterestPoints().size() + " possible matching points");
    //         
    //         image = ImageIO.read(new File("/home/astromberg/test_images/photo.jpg"));
    //         Surf board2 = new Surf(image);
    //         matchingPoints = board.getMatchingPoints(board2,false);
    //         System.out.println("between the first and second image there are " + matchingPoints.size() + " of " + board.getFreeOrientedInterestPoints().size() + " possible matching points");
    //         
    //         image = ImageIO.read(new File("/data/work/OpenSURF/Images/img2.jpg"));
    //         Surf board3 = new Surf(image);
    //         matchingPoints = board.getMatchingPoints(board3,false);
    //         System.out.println("between the first and third image there are " + matchingPoints.size() + " of " + board.getFreeOrientedInterestPoints().size() + " possible matching points");
    //
    //         image = ImageIO.read(new File("/home/astromberg/test_images/photo2.jpg"));
    //         Surf board4 = new Surf(image);
    //         matchingPoints = board.getMatchingPoints(board4,false);
    //         System.out.println("between the first and third image there are " + matchingPoints.size() + " of " + board.getFreeOrientedInterestPoints().size() + " possible matching points");
    //
    //         image = ImageIO.read(new File("/home/astromberg/test_images/photo3.jpg"));
    //         Surf board5 = new Surf(image);
    //         matchingPoints = board.getMatchingPoints(board5,false);
    //         System.out.println("between the first and third image there are " + matchingPoints.size() + " of " + board.getFreeOrientedInterestPoints().size() + " possible matching points");
    //         
    //         image = ImageIO.read(new File("/home/astromberg/test_images/photo4.jpg"));
    //         Surf board6 = new Surf(image);
    //         matchingPoints = board.getMatchingPoints(board6,false);
    //         System.out.println("between the first and third image there are " + matchingPoints.size() + " of " + board.getFreeOrientedInterestPoints().size() + " possible matching points");
    //         
    //         image = ImageIO.read(new File("/home/astromberg/test_images/cover2.jpg"));
    //         Surf board7 = new Surf(image);
    //         matchingPoints = board.getMatchingPoints(board7,false);
    //         System.out.println("between the first and third image there are " + matchingPoints.size() + " of " + board.getFreeOrientedInterestPoints().size() + " possible matching points");
    //         
    //         matchingPoints = board6.getMatchingPoints(board7,false);
    //         System.out.println("between the first and third image there are " + matchingPoints.size() + " of " + board6.getFreeOrientedInterestPoints().size() + " possible matching points");
    //      } catch (Exception e){
    //         e.printStackTrace();
    //      }
    //   }
    public boolean isEquivalentTo(Surf surf) {
        List<SURFInterestPoint> pointsA = surf.getFreeOrientedInterestPoints();
        List<SURFInterestPoint> pointsB = getFreeOrientedInterestPoints();
        if (pointsA.size() != pointsB.size()) {
            return false;
        }
        for (int i = 0; i < pointsA.size(); i++) {
            SURFInterestPoint pointA = pointsA.get(i);
            SURFInterestPoint pointB = pointsB.get(i);
            if (!pointA.isEquivalentTo(pointB)) {
                return false;
            }
        }
        return true;
    }

    public static void saveToFile(Surf surf, String file) {
        try {
            ObjectOutputStream stream = new ObjectOutputStream(new FileOutputStream(file));
            stream.writeObject(surf);
            stream.close();
        } catch (IOException e) {
            e.printStackTrace();
        }
    }

    public static Surf readFromFile(String location) {
        File file = new File(location);
        if (file != null && file.exists()) {
            try {
                ObjectInputStream str = new ObjectInputStream(new FileInputStream(file));
                return (Surf) str.readObject();
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
        return null;
    }

    private void writeObject(ObjectOutputStream out) throws IOException {
        if (mFreeOrientedPoints == null) {
            getFreeOrientedInterestPoints();
        }
        if (mUprightPoints == null) {
            getUprightInterestPoints();
        }
        out.defaultWriteObject();
    }

    /*
    public static void main(String args[]) {
    try {
        BufferedImage image = ImageIO.read(new File("test.png"));
        Surf board = new Surf(image);
        saveToFile(board, "C:\\surf_test.bin");
        Surf boarder = readFromFile("C:\\surf_test.bin");
        List<SURFInterestPoint> points = boarder.getFreeOrientedInterestPoints();
        System.out.println("Found " + points.size() + " interest points");
    } catch (Exception e) {
        e.printStackTrace();
    }
    }*/
}