iflocation.cpp

A tutorial and open source code for finding edges and corners based on the filters used in primary v

// IFLocation.cpp: ImageFeatures class CLocation

/*
** Copyright (C) 1994, 2003 Tyler C. Folsom
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/
#include "stdafx.h"
#include "IFLocation.h"
#include "quad_dis.h"

#include "profile.h"

#define FL4TH  ((float) 0.25)
#define FLHALF ((float) 0.5)
#define FL0    ((float) 0)
#define FL1    ((float) 1)
#define FL2    ((float) 2)
#define FL3    ((float) 3)
#define FL4    ((float) 4)
#define FL5    ((float) 5)
#define FL6    ((float) 6)

/* Solve for angles within 0.1 degree */
#define ANGLE_CONVEGENCE ((float) 0.00174)

#define FACTOR ((float)1.6)

#define EVEN_ORDER 2
#define ODD_ORDER  3
#define G2_ORN (EVEN_ORDER+1)
#define H2_ORN (4)
#define G1_ORN 2
#define H1_ORN 3
#define WAY_SMALL  ((float) -1.0e20)


static float deriv_g2h2( float *gl, float theta);
static float g1h1( float *gl, float theta, float *d1, float *d2);

/*------------------------------------------------------------------------*/
// constructor specifying (x,y) sampling location and
// smaller filter diameter.  Numbers are in pixels.
CLocation::CLocation( int x, int y, int diam )
{
    Initialize( x, y, diam );
}
/*------------------------------------------------------------------------*/
// Initialization specifying (x,y) sampling location and
// smaller filter diameter.  Numbers are in pixels.
void CLocation::Initialize( int x, int y, int diam )
{
    m_x_rf = x;
    m_y_rf = y;
    if (diam > 0)
    {
        m_diam = diam;
        m_usingOdd = false; 
    }
    else
    {
        m_diam = -diam;
        m_usingOdd = true; 
    }
    if ((m_diam & 1) == 0)
        m_big_diam = 2 * m_diam;
    else  /* keep odd filters centered at same point */
        m_big_diam = 2 * m_diam + 1;
    m_usingSmall = true;
}
//----------------------------------------------------------------------------------------
// copy constructor 
CLocation::CLocation(const CLocation &right)
{
    Copy( right );
}
//----------------------------------------------------------------------------------------
//  copy method
void CLocation::Copy(const CLocation &right)
{
    m_x_rf = right.m_x_rf;
    m_y_rf = right.m_y_rf;
    m_usingSmall = right.m_usingSmall;
    m_usingOdd = right.m_usingOdd;
    m_diam = right.m_diam;
    m_big_diam = right.m_big_diam;

    int size = right.corrEven.GetSize();
    corrEven.SetSize(size);
    for (int i = 0; i < size; i++)
    {
        corrEven.SetAt( i, right.corrEven[i]);
    }

    size = right.corrOdd.GetSize();
    corrOdd.SetSize(size);
    for (i = 0; i < size; i++)
    {
        corrOdd.SetAt( i, right.corrOdd[i]);
    }

    size = right.corrBigEven.GetSize();
    corrBigEven.SetSize(size);
    for (i = 0; i < size; i++)
    {
        corrBigEven.SetAt( i, right.corrBigEven[i]);
    }

    size = right.corrBigOdd.GetSize();
    corrBigOdd.SetSize(size);
    for (i = 0; i < size; i++)
    {
        corrBigOdd.SetAt( i, right.corrBigOdd[i]);
    }

}
/*------------------------------------------------------------------------*/
const CLocation &CLocation::operator=( const CLocation &right )
{
    if (&right != this)
    {
        Copy( right );
    }
    return *this;
}
/*------------------------------------------------------------------------*/
// store correlation "result" in slot n of even/odd big/small filter
void CLocation::SetResult( int n, float result)
{
    if (m_usingSmall)
    {
        if (m_usingOdd)
        {
            corrOdd.SetAtGrow( n, result);
        }
        else // small even
        {
            corrEven.SetAtGrow( n, result);
        }
    }
    else  // big filters
    {
        if (m_usingOdd)
        {
            corrBigOdd.SetAtGrow( n, result);
        }
        else // big even
        {
            corrBigEven.SetAtGrow( n, result);
        }
    }
}
/*------------------------------------------------------------------------*/
// return correlation from slot n of even/odd big/small filter
float CLocation::GetResult( int n)
{
    if (m_usingSmall)
    {
        if (m_usingOdd)
        {
            return (corrOdd[n]);
        }
        else // small even
        {
            return (corrEven[n]);
        }
    }
    else  // big filters
    {
        if (m_usingOdd)
        {
            return (corrBigOdd[n]);
        }
        else // big even
        {
            return (corrBigEven[n]);
        }
    }

}
/*------------------------------------------------------------------------*/
// returns a rectangle defining the subimage of the location
RECT CLocation::GetRect()
{
    RECT bounds;
    int low_half;
    int up_half;
    if (m_usingSmall)
    {
        low_half = m_diam / 2;
        up_half = (m_diam+1) / 2;
    }
    else  // big filters
    {
        low_half = m_big_diam / 2;
        up_half = (m_big_diam+1) / 2;
    }
    bounds.left   = m_x_rf - low_half;
    bounds.top    = m_y_rf - low_half;
    bounds.right  = m_x_rf + up_half;
    bounds.bottom = m_y_rf + up_half;
    return bounds;
}

/*------------------------------------------------------------------------*/
// returns an upper bound on the strength
float CLocation::GetStrengthBound()  
{
    float bound = 0;
    for (int i = 0; i < corrOdd.GetSize(); i++)
        bound += FL_ABS(corrOdd[i]);
    for (i = 0; i < corrEven.GetSize(); i++)
        bound += FL_ABS(corrEven[i]);
    return( bound );
}
/*------------------------------------------------------------------------*/
// Steer a filter type to the indicated angle (radians)
float CLocation::Steer( float angle )
{
	PROFILE("Steer");

    float steered = 0;
    if (m_usingSmall)
    {
        if (m_usingOdd)
        {   
            steered = (ODD_ANGLES == 4)? 
                steer_45( corrOdd, angle ): // 4 orientations
                steer_90( corrOdd, angle ); // 2 orientations
        }
        else // small even: 3 orientations
        {
            steered = steer_60( corrEven, angle );
        }
    }
    else  // big filters
    {
        if (m_usingOdd)
        {
            steered = (ODD_ANGLES == 4)? 
                steer_45( corrBigOdd, angle ): // 4 orientations
                steer_90( corrBigOdd, angle ); // 2 orientations
        }
        else // big even: 3 orientations
        {
            steered = steer_60( corrBigEven, angle );
        }
    }
    return steered;
}
/*------------------------------------------------------------------------*/
/* 
    steer_60: Given input filters sampled at 60 degrees, interpolate
    to the filter response at angle theta.

    input:  sampled - response of filters at 0, 60, 120 degrees.
            theta   - angle for desired response in radians.
    returned value: the interpolated filter response at theta.
*/
float CLocation::steer_60( 
    CArray<float, float>& sampled,
    float theta)
{
    float angle;
    float steered;
    int i_angle;

    angle = (float) (PI_1 / G2_ORN);
    steered = (float) 0.0;
    for (i_angle = 0; i_angle < G2_ORN; i_angle++)
        steered += sampled[i_angle] * ((float) 1.0 + (float) 2.0 *
         (float) cos( (double) (2.0 * theta - PI_2 * (float)i_angle / 
         (float) G2_ORN)));
    steered /= G2_ORN;
    return (steered);
}
/*------------------------------------------------------------------------*/
/* 
    steer_45: Given input filters sampled at 45 degrees, interpolate
    to the filter response at angle theta.

    input:  sampled - response of filters at 0, 45, 90, 135 degrees.
            theta   - angle for desired response in radians.
    returned value: the interpolated filter response at theta.
*/
float CLocation::steer_45( 
    CArray<float, float>& sampled,
    float theta)
{
    int angle;
    float steered;

    steered = (float) 0.0;
    for (angle = 0; angle < H2_ORN; angle++)
        steered += sampled[angle] * (float) ( 
         cos( (double) theta - PI_1 * (float)angle / (float) H2_ORN) +
         cos( theta * 3.0 - PI_3 * (float) angle / (float) H2_ORN));
    steered /= (float) 2.0;
    return (steered);
}
/*------------------------------------------------------------------------*/
/* 
    steer_90: Given input filters sampled at 90 degrees, interpolate
    to the filter response at angle theta.

    input:  sampled - response of filters at 0, 90 degrees.
            theta   - angle for desired response in radians.
    returned value: the interpolated filter response at theta.
*/
float CLocation::steer_90( 
    CArray<float, float>& sampled,
    float theta)
{
    return ( (float) (sampled[0] * cos( (double) theta)
                    + sampled[1] * sin( (double) theta)) );
}
/*------------------------------------------------------------------------*/
// Finds the angle (radians) of the feature based on filter correlations
/* 
    GetAngle: Given input filters sampled at 60 and 90 degrees, find the angle
    of the dominant response.

    input:  corrEven - response of even filters at 0, 60, 120 degrees.
            corrOdd - response of odd filters at 0 and 90 degrees.
    output: dom_resp - squared response of the filter at the dominant angle.
            *steeredEven - interpolated even response at the dominant angle.
            *steeredOdd - interpolated odd response at the dominant angle.
    returned value: angle of dominant response in radians.

    There are two methods here.  The older, obsolete, method is based on using
     corrOdd - response of odd filters at 0, 45, 90, 135 degrees.
	The solve() and deriv_g2h2() routines find the angle where the derivative is zero.

    The newer, more efficient, method uses solve_max() and g1h1() to find the 
	angle that gives maximum response.  For documentation, see the header for
	g1h1()  and  Folsom & Pinter, "Primitive Features by Steering, Quadrature
	and Scale", IEEE Trans PAMI, Nov 1998, pp. 1161-1173.

*/
float CLocation::GetAngle( 
    float *dom_resp,
    float *steeredEven,
    float *steeredOdd,
	bool  noSteering) // no steering if true

{
	PROFILE("CLocation::GetAngle");

    float gl[7];  /* gamma and lamda coefficients */
    float coef[8];  /* combinations of gamma and lambda for trig terms*/
    float angles[4];  /* 4 solutions */
    float dom2;
    float dom1;
    float theta1, theta;
    float theta_rad;
    float max_mag, mag_sq;

	if (m_x_rf == 77 && m_y_rf == 87)
	{
		max_mag = 10;  // break here to debug
	}
    max_mag = WAY_SMALL;
    /* Gamma1, gamma2, gamma3 */
    gl[0] = (corrEven[0] + corrEven[1] + corrEven[2]) / FL3;
    gl[1] = (FL2 * corrEven[0] - corrEven[1] - corrEven[2]) / FL3;
    gl[2] = (corrEven[1] - corrEven[2]) / ROOT_3;
    if (ODD_ANGLES == 4)
    {
        /* Lambda1, .. lambda3 */
        gl[3] = (corrOdd[0] + (corrOdd[1] - corrOdd[3]) / ROOT_2) / FL2;
        gl[4] = (corrOdd[2] + (corrOdd[1] + corrOdd[3]) / ROOT_2) / FL2;
        gl[5] = (corrOdd[0] - (corrOdd[1] - corrOdd[3]) / ROOT_2) / FL2;
        gl[6] = (-corrOdd[2] + (corrOdd[1] + corrOdd[3]) / ROOT_2) / FL2;
        coef[0] = FL2 * (FL2*gl[0]*gl[2] + gl[3]*gl[4] - gl[4]*gl[5] + gl[3]*gl[6]);
        coef[1] = -(FL4*gl[0]*gl[1] + gl[3]*gl[3] - gl[4]*gl[4] + FL2*gl[3]*gl[5]
         + FL2*gl[4]*gl[6]);
        coef[2] = FL4 * (gl[1]*gl[2] + gl[4]*gl[5] + gl[3]*gl[6]);
        coef[3] = FL2 * (-gl[1]*gl[1] + gl[2]*gl[2] - FL2*gl[3]*gl[5]
         + FL2*gl[4]*gl[6]);
        coef[4] = FL6 * gl[5]*gl[6];
        coef[5] = FL3*(-gl[5]*gl[5] + gl[6]*gl[6]);

        /* try arbitrary initial guess for theta */
        for (theta = FL0; theta < PI_1;)
        {
            theta1 = theta;
		    if (!noSteering)
		    {	// Find the strongest angle near this guess
			    solve (coef, &theta1, deriv_g2h2);
		    }
            dom1 = steer_45 (corrOdd, theta1);
            dom2 = steer_60 (corrEven, theta1);
            mag_sq = dom1 * dom1 + dom2 * dom2;
            if (mag_sq >= max_mag)
            {
                max_mag = mag_sq;
                theta_rad = theta1;
                *steeredEven = dom2;
                *steeredOdd = dom1;
            }
		    if (noSteering)
		    {	// skip if not looking for the angle.
			    break;
		    }
            /* start next guess about 1 degree beyond previous root. */
            theta = theta1 + (float) 0.02;
        }
        *dom_resp = max_mag;
        return (theta_rad);
    }
    else
    {