iffilter.cpp

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

// IFFilter.cpp: ImageFeatures class CFilter

/*
** 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 "IFFilter.h"
#include "quad_dis.h"
#include "profile.h"

CFilter::CFilter( int scale, int spacing)  // constructor
{
    // Use positive diameters for even filters,
    // negative diameters for odd.
    // Scale = 0 for the default constructor, but this
    // is temporary and not a usable value.
    if (scale > 0)
    {
        m_isEven = true;
        m_diam = scale;
        m_orientations = 3;    // Number of even filters
        m_coef[0] =  (float) -502.48;      /* x*x term */
        m_coef[1] =  (float) 0.0;          /* x   term */
        m_coef[2] =  (float) 7.8287;       /* constant */
    }
    else
    {
        m_isEven = false;
        m_diam = -scale;        // diameter of receptive field (pixels)
        m_orientations = ODD_ANGLES;    // Number of odd filters
        if (ODD_ANGLES == 4)
        {
            m_coef[0] =   (float) -925.81;      /* x*x*x term */
            m_coef[1] =   (float) 0.0;          /* x*x term */
            m_coef[2] =   (float) 97.7913;      /* x   term */
            m_coef[3] =   (float) 0.0;          /* constant */
        }
        else
        {
            m_coef[0] =   (float) 72.0232;      /* x   term */
            m_coef[1] =   (float) 0.0;          /* constant */
        }
    }
    // The overlap terms are not filled in until needed
    m_center_lobe = 0;  // indicates no overlap terms.
	m_sampleSpacing = (spacing>0 && 4*spacing >= m_diam)? spacing:1;  // 1 for no subsampling
//  ASSERT( m_diam > 1);
    if (m_diam > 1)
    {
		int diam = m_diam / m_sampleSpacing;
        kern = new float[m_orientations * diam * diam];
        setKernel();
        fl_correlate (*this); // self overlaps
        if (m_isEven)
        {   // compute the width of the even kernel.
			 
            int mid = (diam & 1)?
				diam * diam / 2:  // odd numbered diameter
				diam * (diam + 1) / 2 - 1;  // even numbered diameter
            float center_value = kern[mid];
            for (int i = mid; center_value * kern[i] > (float) 0; i++)
                ;
			int right = i;  // the index where kernel sign no longer matches center
            /* kern[] changes signs between (right-1) and right */
            /* Now interpolate the distance to the sign change. */
            float delta = kern[right-1] / (kern[right-1] - kern[right]);
            for (i = mid; center_value * kern[i] > (float) 0; i--)
                ;
			/* Add in the interpolated distance at the other end */
			delta += kern[i+1] / (kern[i+1] - kern[i]);
            m_center_lobe = (delta + (right - i - 2)) * m_sampleSpacing;
        }
    }
}
/*------------------------------------------------------------------------*/
CFilter::~CFilter( )  // destructor
{
    if (m_diam > 1)
    {
        delete[] kern;
    }
}
//----------------------------------------------------------------------------------------
// copy constructor 
CFilter::CFilter(const CFilter &right)
{
    if (right.m_diam > 1)
    {
        kern = new float[right.m_orientations * 
            (right.m_diam / right.m_sampleSpacing) * 
			(right.m_diam / right.m_sampleSpacing)];
    }
    Copy( right );
}
/*------------------------------------------------------------------------*/
const CFilter &CFilter::operator=( const CFilter &right )
{
    if (&right != this)
    {    
        if (m_diam > 1)
        {
            delete[] kern;
        }

        if (right.m_diam > 1)
        {
            kern = new float[right.m_orientations * 
                    (right.m_diam / right.m_sampleSpacing) * 
					(right.m_diam / right.m_sampleSpacing)];
        }
        Copy( right );
    }
    return *this;
}
//----------------------------------------------------------------------------------------
//  copy method
void CFilter::Copy(const CFilter &right)
{
    m_diam = right.m_diam;
    m_orientations = right.m_orientations;
    m_isEven = right.m_isEven;
    m_center_lobe = right.m_center_lobe;
	m_sampleSpacing = right.m_sampleSpacing;
	m_maxResponse = right.m_maxResponse;
    int i;
    for (i = 0; i < m_orientations; i++)
    {
        m_coef[i] = right.m_coef[i];
        m_self_ovlap[i] = right.m_self_ovlap[i];
        m_cross_ovlap[i] = right.m_cross_ovlap[i];   
    }
    int size = m_orientations * (m_diam/m_sampleSpacing) * (m_diam/m_sampleSpacing);
    if (m_diam > 1)
    {
        for (i = 0; i < size; i++)
        {
            kern[i] = right.kern[i];
        }
    }
}
/*----------------------------------------------------------------------*/
// construct the kernels
void CFilter::setKernel()
{
	PROFILE("setKernel");

    /* kernel is a square of "m_diam" by "m_diam" */
    int i, j, k, orn;
    float poly;
    float u [MAX_COEFS], v [MAX_COEFS], x, y, z;
    float gaus, gaus1, gaus2;
    float th [MAX_COEFS];
    float incr_x, incr_y, x_min, y_min;

    for (orn = 0; orn < m_orientations; orn++)
    {
        th [orn] = (float) (orn * PI_1 / m_orientations);
        u[orn] = (float) cos( (double) th[orn]);
        v[orn] = (float) sin( (double) th[orn]);
    }
    incr_x = (float) (2.0 * EFF_LIM * m_sampleSpacing / (float)( m_diam + 1));
    /* It is assumed that the filter values are essentially zero at
       +/-EFF_LIM; thus values are not computed there. */
    x_min = -EFF_LIM + incr_x;
    incr_y = (float) (2.0 * EFF_LIM * m_sampleSpacing / (float)(m_diam + 1));
//  y_max = EFF_LIM - incr_y;
	y_min = -EFF_LIM + incr_y;
    int index = 0;

// new version 5/19/03 with x in inner loop
    for (orn = 0; orn < m_orientations; orn++)
    {
        y = y_min;  
        for (i = 0; i < m_diam; i += m_sampleSpacing)
        {
            gaus1 = bump(y);
            x = x_min;
            for (j = 0; j < m_diam; j += m_sampleSpacing)
            {
                z = (float) (u[orn] * x + v[orn] * y);
                poly = (float) 0;
                for (k = 0; k < m_orientations - 1; k++)
                    poly = z * (poly + m_coef[k]);
                poly += m_coef[k];

                gaus2 = bump(x);
                gaus = gaus1 * gaus2;
                kern[index++] = (float) gaus * poly;
                x += incr_x;
            }
            y += incr_y;
        }
    }
}
/*------------------------------------------------------------------------*/
/* This is used for correlating a kernel with another kernel. */
/* Given a pair of quadrature phase kernels at different orientations,
   compute how each phase and orientation responds to an image identical
   to a kernel.
   Results are placed either in m_self_ovlap or m_cross_ovlap
 */
void CFilter::fl_correlate (CFilter& other)
{
	PROFILE("fl_correlate");


    int size = (m_diam / m_sampleSpacing) * (m_diam / m_sampleSpacing);
    float *pLast = kern + size;

	float positive = 0;
	float negative = 0;
    for (float *pKern = kern; pKern < pLast; pKern++)
    {
        if(*pKern > 0)
			positive += *pKern;
		else
			negative += *pKern;
    }
    m_maxResponse = MAX_PIXEL * positive - MIN_PIXEL * negative;

	if (m_diam != other.m_diam)
    return;  // mismatched filter sizes
    
	float *pResults;
    if (m_isEven == other.m_isEven)
    {
        pResults = m_self_ovlap;
    }
    else
    {
        pResults = m_cross_ovlap;
    }
    float *pOther = other.kern;
	int maxOrn = min( m_orientations, other.m_orientations);
    for (int orn = 0; orn < maxOrn; orn ++)
    {
        *pResults = (float) 0.0;
        for (float *pKern = kern; pKern < pLast; pKern++)
        {
                *pResults += *pKern * *pOther++;
        }
        pResults++;
    }
}

/*----------------------------------------------------------------------*/
float CFilter::bump( float x)
/*
  This function generates a bump function mapping the interval
  [-sqrt(LIM_SQ),sqrt(LIM_SQ)] onto [0,1].
  If called with a line such as "y=bump(x)", where x
  has been specified, it will automatically return the corresponding value.
  bump(x) is supported on the interval specified by LIM_SQ, clamped to 0
  outside it, and normalized to peak at y(x) = 1.0 when x = 0.0

*/
{
    float xx;

    xx = x * x;
    if (xx >= (float) (LIM_SQ - EPSILON))
        return ((float) 0.0);
    xx = (float) exp( (double) (4.0 * xx / (xx - LIM_SQ)));
    return (xx);
}