www.pudn.com > imageFeatures_Ver3_0_source.zip > IFFilter.cpp
// 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);
}