www.pudn.com > imageFeatures_Ver3_0_source.zip > IFFeature.cpp
// IFFeature.cpp: ImageFeatures class CFeature
// After we have done the convolution and steered the filters,
// the routines in this module process the 1D results to
// find the location, type, strength and width of the feature.
/*
** 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 "IFFeature.h"
#include "quad_dis.h"
#include "profile.h"
/*------------------------------------------------------------------------*/
CFeature::CFeature(
int x_rf, int y_rf, // sampling location
int diam, // diameter of small filters
float angle, // angle of maximum response (radians)
// the following four correlations are steered to angle
float corrEven, // correlation with small even filter
float corrOdd, // correlation with small odd filter
float corrBigEven, // correlation with big even filter
float corrBigOdd, // correlation with big odd filter
float base_strength, // relative strength of small filters
float big_base_strength, // relative strength of big filters
float big_diam, // diameter of big filters
float center_lobe, // width of big even filter
float width_tuning) // average width of big and small
{
// copy information found by CLocation
m_x_rf = x_rf;
m_y_rf = y_rf;
m_diam_rf = diam;
/* adjust the center position for an even diameter. */
if ((m_diam_rf & 1) == 0)
{
m_x_rf -= (float) 0.5;
m_y_rf -= (float) 0.5;
}
m_main.m_x = m_x_rf;
m_main.m_y = m_y_rf;
m_angle = angle;
m_main.m_corrEven = corrEven;
m_main.m_corrOdd = corrOdd;
m_main.m_corrBigEven = corrBigEven;
m_main.m_corrBigOdd = corrBigOdd;
m_main.m_degrees = m_angle * (float) 180.0 / PI_1;
while (m_main.m_degrees < 0)
m_main.m_degrees += (float) 360;
while (m_main.m_degrees >= (float) 360)
m_main.m_degrees -= (float) 360;
m_base_strength = base_strength;
m_big_base_strength = big_base_strength;
m_scale = (float) ((float)(diam + 1) / (2.0 * EFF_LIM));
m_big_scale = (float) ((float)(big_diam + 1) /
(2.0 * EFF_LIM));
m_center_lobe = center_lobe;
m_width_tuning = width_tuning;
m_lateral = 0;
for (int i = 0; i < NEIGHBORS; i++)
m_pNbr[i] = 0;
}
//----------------------------------------------------------------------------------------
// copy constructor
CFeature::CFeature(const CFeature &right)
{
Copy( right );
}
/*------------------------------------------------------------------------*/
const CFeature &CFeature::operator=( const CFeature &right )
{
if (&right != this)
{
Copy( right );
}
return *this;
}
//----------------------------------------------------------------------------------------
// copy method
void CFeature::Copy(const CFeature &right)
{
m_x_rf = right.m_x_rf;
m_y_rf = right.m_y_rf;
m_diam_rf = right.m_diam_rf;
m_angle = right.m_angle;
m_length = right.m_length;
m_base_strength = right.m_base_strength;
m_big_base_strength = right.m_big_base_strength;
m_scale = right.m_scale;
m_big_scale = right.m_big_scale;
m_center_lobe = right.m_center_lobe;
m_width_tuning = right.m_width_tuning;
m_main = right.m_main;
m_90 = right.m_90;
// m_45 = right.m_45;
// m_135 = right.m_135;
m_corner = right.m_corner;
m_lateral = right.m_lateral;
for (int i = 0; i < NEIGHBORS; i++)
m_pNbr[i] = right.m_pNbr[i];
}
//----------------------------------------------------------------------------------------
// copy constructor
Data_1D::Data_1D(const Data_1D &right)
{
Copy( right );
}
/*------------------------------------------------------------------------*/
const Data_1D &Data_1D::operator=( const Data_1D &right )
{
if (&right != this)
{
Copy( right );
}
return *this;
}
//----------------------------------------------------------------------------------------
// copy method
void Data_1D::Copy(const Data_1D &right)
{
m_degrees = right.m_degrees;
m_corrEven = right.m_corrEven;
m_corrOdd = right.m_corrOdd;
m_corrBigEven = right.m_corrBigEven;
m_corrBigOdd = right.m_corrBigOdd;
m_strength = right.m_strength;
m_x = right.m_x;
m_y = right.m_y;
m_pos = right.m_pos;
m_type = right.m_type;
m_width = right.m_width;
m_cos_th = right.m_cos_th;
m_sin_th = right.m_sin_th;
}
//----------------------------------------------------------------------------------------
// copy constructor
Data_2D::Data_2D(const Data_2D &right)
{
Copy( right );
}
/*------------------------------------------------------------------------*/
const Data_2D &Data_2D::operator=( const Data_2D &right )
{
if (&right != this)
{
Copy( right );
}
return *this;
}
//----------------------------------------------------------------------------------------
// copy method
void Data_2D::Copy(const Data_2D &right)
{
m_x = right.m_x;
m_y = right.m_y;
m_type = right.m_type;
}
/*------------------------------------------------------------------------*/
// This allows us to steer to another angle and save the corelations
void CFeature::AddSecondary
( int i,
float angle, // angle of secondary response (radians)
float corrEven, float corrOdd,
float corrBigEven, float corrBigOdd)
{
Data_1D* pData;
switch (i)
{
default:
case 0:
pData = &m_main;
break;
case 1:
// pData = &m_45;
break;
case 2:
pData = &m_90;
break;
case 3:
// pData = &m_135;
break;
}
pData->m_degrees = angle * (float) 180.0 / PI_1;
while (pData->m_degrees < 0)
pData->m_degrees += (float) 360;
while (pData->m_degrees >= (float) 360- EPSILON)
pData->m_degrees -= (float) 360;
if (pData->m_degrees < 0)
pData->m_degrees = 0;
pData->m_corrEven = corrEven;
pData->m_corrOdd = corrOdd;
pData->m_corrBigEven = corrBigEven;
pData->m_corrBigOdd = corrBigOdd;
}
/*------------------------------------------------------------------------*/
void CFeature::multi_find( int i, bool edgesOnly)
{
switch (i)
{
default:
case 0:
find_pos( &m_main, edgesOnly );
m_corner.m_type = m_main.m_type;
m_corner.m_x = m_main.m_x;
m_corner.m_y = m_main.m_y;
m_angle = m_main.m_degrees * PI_1 / 180;
break;
case 1:
// find_pos( &m_45, edgesOnly );
break;
case 2:
find_pos( &m_90, edgesOnly );
break;
case 3:
// find_pos( &m_135, edgesOnly );
break;
}
}
/*------------------------------------------------------------------------*/
void CFeature::find_pos( Data_1D* pData, bool edgesOnly )
{
PROFILE("find_pos");
/* find_pos: given a convolution of an image with even G2 and odd H2
filters at two different scales, find the type of edge, the edge position
edge strength and bar width.
The position is in pixels from the center of the receptive field.
The results are placed in a structure. */
// float m_corrEven, /* convolution with G2 at small scale */
// float m_corrOdd, /* convolution with H2 at small scale */
// float m_corrBigEven, /* convolution with G2 at twice scale */
// float m_corrBigOdd, /* convolution with H2 at twice scale */
// float m_angle,
// float m_base_strength, /* strength is relative to this. */
// float m_big_base_strength, /* or else relative to this. */
// float scale,
// float big_scale,
// float center_lobe,
// float width_tuning, /* use m_corrEven for bar widths less than this */
float magn_small, magn_big;
double phase_big, offset, phase_small;
float small_white_pos, small_black_pos, small_step_pos;
float big_white_pos, big_black_pos, big_step_pos;
float white_diff, black_diff, step_diff;
float pos, strength, big_strength, bar_width;
float ratio, skinnyFactor;
/* code was originally written for DARK_BAR = low, high, low, etc.
However, white = 0 and black = 255.
Thus the names are reversed. */
magn_small = (float) sqrt((double)(pData->m_corrEven * pData->m_corrEven +
pData->m_corrOdd * pData->m_corrOdd));
magn_big = (float) sqrt((double)(pData->m_corrBigEven * pData->m_corrBigEven +
pData->m_corrBigOdd * pData->m_corrBigOdd));
if (magn_small < EPSILON || magn_big < EPSILON)
{
pData->m_type = eNO_FEATURE;
return;
}
phase_big = atan2 ((double)pData->m_corrBigOdd, (double)pData->m_corrBigEven);
phase_small = atan2 ((double)pData->m_corrOdd, (double)pData->m_corrEven);
/* From 2/3/00 to 5/27/03, used
an alternative method based on assuming
-PI/2 < white bar phase < PI/2
| black bar phase | > PI/2
Systematic tests of the two methods have not been done,
but it seems that this assumption is not always valid,
especially for narrow bars.
*/
/* 3-way position test */
small_white_pos = m_scale * find_bar (phase_small);
offset = phase_small < 0.0 ? -PI_1 : PI_1;
small_black_pos = m_scale * find_bar (phase_small - offset);
big_white_pos = m_big_scale * find_bar (phase_big);
offset = phase_big < (float)0 ? -PI_1 : PI_1;
big_black_pos = m_big_scale * find_bar (phase_big - offset);
white_diff = small_white_pos - big_white_pos;
white_diff = FL_ABS(white_diff);
black_diff = small_black_pos - big_black_pos;
black_diff = FL_ABS(black_diff);
small_step_pos = m_scale * find_step (phase_small);
big_step_pos = m_big_scale * find_step (phase_big);
step_diff = small_step_pos - big_step_pos;
step_diff = FL_ABS(step_diff);
/* 3-way test */
if (white_diff < black_diff)
{
if (white_diff < step_diff)
pData->m_type = eWHITE_LINE;
else
pData->m_type = eEDGE;
}
else
{
if (black_diff < step_diff)
pData->m_type = eBLACK_LINE;
else
pData->m_type = eEDGE;
}
/* A phase < 0.43 radians can only be produced by a white bar.
Phases this small for light bar or edge occur when magnitude is close to zero. */
if (fabs(phase_small) < 0.43 || fabs(phase_big) < 0.43)
pData->m_type = eWHITE_LINE;
/* A phase > 2.7 radians can only be produced by a black bar.
Phases this small for white bar or edge occur when magnitude is close to zero. */
if (fabs(phase_small) > 2.7 || fabs(phase_big) > 2.7)
pData->m_type = eBLACK_LINE;
if (edgesOnly)
{
pData->m_type = eEDGE;
if (fabs(phase_small) < 0.43 || fabs(phase_big) < 0.43 ||
fabs(phase_small) > 2.7 || fabs(phase_big) > 2.7)
{
pData->m_type = eNO_FEATURE;
return;
}
}
// loop back a second time if initial classification is wrong
for (bool done = false; !done; )
{
switch (pData->m_type)
{
case eEDGE:
strength =
find_strength (magn_big, (float) phase_big, pData->m_type)
/ m_big_base_strength;
bar_width = (float) 0.0;
if (fabs(phase_small) < PI * 0.25 ||
fabs(phase_small) > PI * 0.75)
{ // use big filter at edge of receptive field
pos = big_step_pos;
if (phase_big < 0.0 && pData == &m_main)
{ // A light to dark step.
// If the angle had been 90, change to 270.
// An angle of 180 changes to 0.
pData->m_degrees += 180;
while (pData->m_degrees >= (float) 360)
pData->m_degrees -= (float) 360;
// At this point m_angle (in radians) no longer
// matches m_degrees!
// Of course it isn't supposed to match if we are not working with m_main
}
}
else
{ // use small filter in center of receptive field
pos = small_step_pos;
if (phase_small < 0.0)
{ // A light to dark step.
// If the angle had been 90, change to 270.
// An angle of 180 changes to 0.
pData->m_degrees += 180;
while (pData->m_degrees >= (float) 360)
pData->m_degrees -= (float) 360;
// At this point m_angle (in radians) no longer
// matches m_degrees!
}
}
break;
case eBLACK_LINE:
ratio = magn_small / magn_big;
bar_width = find_width (ratio, m_center_lobe, small_black_pos, &skinnyFactor);
strength =
find_strength (magn_small, (float) phase_small, pData->m_type)
/ (m_base_strength * skinnyFactor);
big_strength =
find_strength (magn_big, (float) phase_big, pData->m_type)
/ (m_big_base_strength * skinnyFactor);
if (big_strength < EPSILON)
{
pData->m_type = eNO_FEATURE;
return;
}
// ratio = strength / big_strength;
// bar_width = find_width (ratio, m_center_lobe);
if (fabs(phase_small) > PI * 0.25 || bar_width > m_width_tuning)
{ // at the edge of the receptive field, use big filter
pos = big_black_pos;
strength = big_strength;
}
else
pos = small_black_pos;
/* If the bar is not centered, find the width from the edge position. */
/* When the phase is smaller than PI - 0.55, the edge position is usable
and the bar is not centered. */
if (fabs(phase_small) < 2.59)
bar_width = 2 * fabs(pos - small_step_pos);
break;
case eWHITE_LINE:
ratio = magn_small / magn_big;
bar_width = find_width (ratio, m_center_lobe, small_white_pos, &skinnyFactor);
strength =
find_strength (magn_small, (float) phase_small, pData->m_type)
/ (m_base_strength * skinnyFactor);
big_strength =
find_strength (magn_big, (float) phase_big, pData->m_type)
/ (m_big_base_strength * skinnyFactor);
if (big_strength < EPSILON)
{
pData->m_type = eNO_FEATURE;
return;
}
// ratio = strength / big_strength;
// bar_width = find_width (ratio, m_center_lobe);
if (fabs(phase_small) < PI * 0.75 || bar_width > m_width_tuning)
{ // at the edge of the receptive field, use big filter
pos = big_white_pos;
strength = big_strength;
}
else
pos = small_white_pos;
/* If the bar is not centered, find the width from the edge position. */
/* When the phase is bigger than 0.55, the edge position is usable
and the bar is not centered. */
if (fabs(phase_small) > 0.55)
bar_width = 2 * fabs(pos - small_step_pos);
break;
case eNO_FEATURE:
default:
break;
} // end of switch on m_type
done = true;
if ((pData->m_type == eWHITE_LINE || pData->m_type == eBLACK_LINE) &&
(pos - bar_width * 0.5 < -(m_diam_rf + 1) * 0.5 ||
pos + bar_width * 0.5 > (m_diam_rf + 1) * 0.5 ))
{ // bar with one side not in small filter
done = false;
pData->m_type = eEDGE;
}
} // while !done
pData->m_pos = pos;
if (FL_ABS(pos) > m_scale * EFF_LIM)
{ /* feature is out of range of smaller filter */
pData->m_type = eNO_FEATURE;
}
else
{
/* The graphics origin is at the upper left.
X axis points right and Y axis points down.
0 degrees points down and a positive angle is counter-clockwise.
The edge location is at (pos, 0) and is rotated ccw..
*/
if (m_x_rf == 61.5 && m_y_rf == 94.5)
{
pData->m_width = bar_width; // debug breakpoint
}
pData->m_cos_th = (float) cos(pData->m_degrees * PI / 180.0);
pData->m_sin_th = (float) sin(pData->m_degrees * PI / 180.0);
if (pData->m_degrees < 180)
{ // for vertical edge, dark side is to left, light to right
pData->m_x = m_x_rf + pos * pData->m_cos_th;
pData->m_y = m_y_rf + pos * pData->m_sin_th;
}
else // phase shift by 180 degrees
{ // for vertical edge, light side is to left, dark to right
pData->m_x = m_x_rf - pos * pData->m_cos_th;
pData->m_y = m_y_rf - pos * pData->m_sin_th;
}
pData->m_width = bar_width;
}
pData->m_strength = strength;
}
/*------------------------------------------------------------------------*/
/* find_strength: given the type of edge, its phase and magnitude,
compute what the magnitude would have been if the edge had been in
the center of the receptive field. */
float CFeature::find_strength
( float magn,
float phase,
enum feature_type type
)
{
/* step edge roll-off profile */
#define ROLL_OFF_ENTRIES 32
static double phases [ROLL_OFF_ENTRIES] = { 0.0,
0.036591, 0.072608, 0.109822, 0.148189, 0.187701, 0.228361, 0.270176,
0.313153, 0.357297, 0.402608, 0.449082, 0.496713, 0.545485, 0.59538,
0.646372, 0.698431, 0.751518, 0.805593, 0.860606, 0.916505, 0.97323,
1.03072, 1.088908, 1.147723, 1.207094, 1.266944, 1.327196, 1.387771,
1.448591, 1.509573, 1.570797};
static double mags [ROLL_OFF_ENTRIES] = { 0.0,
0.028628, 0.042934, 0.061335, 0.084061, 0.111176, 0.142581, 0.178033,
0.21716, 0.259494, 0.30449, 0.351561, 0.400094, 0.449474, 0.499103,
0.548411, 0.596867, 0.643985, 0.689329, 0.732513, 0.773201, 0.811105,
0.845982, 0.877626, 0.905871, 0.930581, 0.951647, 0.968984, 0.982528,
0.992231, 0.998061, 1.0};
/*Data for G1-H1 bar response; Position 0.46 to 0
Phase */
#define BAR_ENTRIES 48
static double bar_phase [BAR_ENTRIES] = {3.1416,
2.709503974, 2.697405066, 2.684777187, 2.671573382, 2.657745581, 2.643243604,
2.628014402, 2.61200147, 2.595144345, 2.577378163, 2.558633276, 2.538834912,
2.517902907, 2.495751511, 2.472289327, 2.447419396, 2.421039521, 2.393042872,
2.363318825, 2.331752406, 2.298216642, 2.262552309, 2.224538371, 2.183863909,
2.140110826, 2.092750398, 2.041151844, 1.984600117, 1.92232223, 1.853524781,
1.777448495, 1.693446865, 1.601093519, 1.500314211, 1.391522873, 1.275720046,
1.154498225, 1.029910499, 0.904205544, 0.779494823, 0.657454115, 0.539142619,
0.424964191, 0.314739594, 0.207835336, 0.103302471, 0.0};
static double bar_mags [BAR_ENTRIES] = {0.0,
0.00173125, 0.004018206, 0.00815523, 0.014898029, 0.025015096, 0.03921154,
0.058062909, 0.081965835, 0.111108148, 0.145457911, 0.184768988, 0.228599778,
0.276341541, 0.327252943, 0.380497838, 0.435183859, 0.490399867, 0.545250812,
0.598888709, 0.650536447, 0.69949525, 0.745127825, 0.786826269, 0.823988638,
0.856024887, 0.882397714, 0.90268979, 0.916682086, 0.924427583, 0.926307491,
0.923060311, 0.915776282, 0.905850871, 0.89489219, 0.884580925, 0.876490181,
0.871887226, 0.871555052, 0.875678763, 0.883830297, 0.895056179, 0.908041494,
0.921306347, 0.933394247, 0.943027408, 0.949220495, 0.951355027};
/*Data for G1-H1 edge response; position -0.37 to 0
Phase */
#define EDGE_ENTRIES 39
static double edge_phase [EDGE_ENTRIES] = {0.0,
0.410874286, 0.418921308, 0.42932989, 0.440681409, 0.452713616, 0.465384309,
0.47870702, 0.492714419, 0.507448578, 0.522958036, 0.539296881, 0.556524524,
0.574705726, 0.593910734, 0.614215438, 0.635701501, 0.658456414, 0.682573444,
0.708151388, 0.73529409, 0.764109624, 0.794709041, 0.827204564, 0.861707105,
0.898322957, 0.937149556, 0.978270205, 1.021747747, 1.067617264, 1.115878073,
1.166485447, 1.219342788, 1.27429516, 1.331125276, 1.389553042, 1.44923954,
1.509795892, 1.570796788};
static double edge_mags [EDGE_ENTRIES] = {0.0,
0.000106765, 0.000401429, 0.001177729, 0.002873849, 0.006085781, 0.011525575,
0.019955297, 0.032111702, 0.048635061, 0.070011394, 0.096532849, 0.128277275,
0.165105526, 0.206673605, 0.252456201, 0.30177812, 0.353850448, 0.407808742,
0.462751049, 0.517774065, 0.572006179, 0.624636484, 0.674939173, 0.722292909,
0.766194948, 0.806269892, 0.842273028, 0.874088269, 0.90172079, 0.925284532,
0.944984881, 0.961097056, 0.973940998, 0.983853892, 0.991161859, 0.996152632,
0.999051303, 1.0};
float roll_off;
phase = FL_ABS(phase);
while (phase > PI_1)
{
if (phase < (float) (PI_1 + 0.01))
{
phase = PI_1; /* avoid round-off errors */
break;
}
else
phase -= PI_1;
}
if (ODD_ANGLES == 4)
{
if ((type == eWHITE_LINE && phase <= (float) (0.5 * PI_1)) ||
(type == eBLACK_LINE && phase >= (float) (0.5 * PI_1)))
return (magn);
/* use the step edge roll-off profile */
if (phase > (float)(0.5 * PI_1))
phase = PI_1 - phase;
roll_off = (float)look_up( (double) phase, ROLL_OFF_ENTRIES, phases, mags);
}
else /* 2 odd orientations */
{
switch (type) {
case eBLACK_LINE:
phase = PI_1 - phase; /* Then handle like white bar */
case eWHITE_LINE:
roll_off = (float)look_up( (double) phase, BAR_ENTRIES, bar_phase, bar_mags);
roll_off /= bar_mags [BAR_ENTRIES-1];
break;
default:
case eEDGE:
if (phase > PI_HALF)
phase = PI - phase;
roll_off = (float)look_up( (double) phase, EDGE_ENTRIES, edge_phase, edge_mags);
break;
}
}
if (roll_off > MIN_ROLL_OFF)
magn /= roll_off;
else
magn = (float) 0;
return (magn);
}
/*------------------------------------------------------------------------*/
/* find_width: given the ratio of filters and the position in the receptive
field, find the bar width.
The return value is the bar width as a fraction of the center lobe
width of the larger even filter.
It also returns *faint, which is the factor by which the strength of
skinny bars appears reduced.
*/
float CFeature::find_width (float ratio, float center_lobe, float pos, float *faint)
{
#define WIDTH_DATA 6
#define WIDTH2 13
#define ECCENTRICITIES 6
/* Bar width relative to width of central lobe of larger filter */
static double width [WIDTH_DATA] =
{ 0.02, 0.06, 0.14, 0.3, 0.46, 0.82 };
static double width2 [WIDTH2] =
{0.01, 0.02, 0.02, 0.03, 0.04, 0.06, 0.07, 0.14, 0.15, 0.23, 0.3, 0.41, 0.46};
static double eccentric [ECCENTRICITIES] =
{ 0.0, 0.1, 0.2, 0.25, 0.3, 0.4};
static double ratios [ECCENTRICITIES] [WIDTH_DATA] =
{
/* Ratio of smaller filter to larger filter magnitude at 0.0 eccentricity */
{2.716531, 2.705151, 2.648604, 2.406888, 2.017341, 0.8790814},
/* Ratio of smaller filter to larger filter magnitude at 0.1 eccentricity */
{2.190965, 2.185212, 2.156598, 2.033694, 1.833451, 1.2196309},
/* Ratio of smaller filter to larger filter magnitude at 0.2 eccentricity */
{2.33697, 2.325586, 2.269978, 2.05062, 1.760186, 1.2660889},
/* Ratio of smaller filter to larger filter magnitude at 0.25 eccentricity */
{1.563844, 1.559338, 1.537893, 1.461619, 1.373351, 1.1737574},
/* Ratio of smaller filter to larger filter magnitude at 0.3 eccentricity */
{0.538405, 0.546242, 0.583861, 0.720047, 0.868897, 0.8688967},
/* Ratio of smaller filter to larger filter magnitude at 0.4 eccentricity */
{1.04E-05, 6.13E-05, 0.001029, 0.022527, 0.109461, 0.1094611}
};
/* Reduction in strength */
static double reduction [WIDTH2] =
{0.029849, 0.059675, 0.059675, 0.089447, 0.119164, 0.178188, 0.207528,
0.406118, 0.433697, 0.636638, 0.783941, 0.958885, 1.0};
float eccentricity = FL_ABS(pos) * 2 / m_diam_rf * EFF_LIM;
for (int ec = 1; ec < ECCENTRICITIES && eccentricity >= eccentric[ec]; ec++)
;
float lo_width = (float) look_up( (double) ratio, WIDTH_DATA, ratios[ec-1], width);
float hi_width = (float) look_up( (double) ratio, WIDTH_DATA, ratios[ec], width);
float wide = lo_width + (hi_width - lo_width) *
(eccentricity - eccentric[ec-1]) / (eccentric[ec] - eccentric[ec-1]);
*faint = (float) look_up( (double) wide, WIDTH2, width2, reduction);
wide *= center_lobe;
return (wide);
}
/*------------------------------------------------------------------------*/
/* find_width: given the ratio of filters, find the bar width.
The return value is the bar width as a fraction of the center lobe
width of the larger even filter. [Obsolete]
*/
float CFeature::find_width (float ratio, float center_lobe)
{
/* bar width profile */
#define WIDTH_ENTRIES 22
static double ratios [WIDTH_ENTRIES] = { 2.0,
1.991631, 1.958355, 1.892638, 1.796176, 1.671565, 1.522355, 1.353106,
1.169414, 0.97789, 0.78604, 0.601984, 0.433956, 0.289524, 0.174502,
0.091701, 0.039829, 0.013206, 0.003083, 0.000703, 0.000486, 0.0};
static double widths [WIDTH_ENTRIES] = { 0.01,
0.040058, 0.120173, 0.200288, 0.280404, 0.360519, 0.440635, 0.52075,
0.600865, 0.680981, 0.761096, 0.841211, 0.921327, 1.001442, 1.081557,
1.161673, 1.241788, 1.321904, 1.402019, 1.482134, 1.56225, 1.64};
float width;
width = (float)look_up( (double) ratio, WIDTH_ENTRIES, ratios, widths);
width *= center_lobe;
return (width);
}
/*------------------------------------------------------------------------*/
/* find_bar: given the phase, find the bar position.
The return value is the bar position.
*/
float CFeature::find_bar (double phase)
{
#define BAR_ENTRIES4 87
#define BAR_ENTRIES2 113
/* Bar phase: 4 odd filters -pi to pi matches -0.42 to 0.42 */
static double bar_phase4 [BAR_ENTRIES4] = {
-3.1416, -3.12663, -3.09061, -3.05356, -3.01546, -2.97628, -2.93599,
-2.89457, -2.85202, -2.80831, -2.76345, -2.71744, -2.67027, -2.62196,
-2.57252, -2.52197, -2.47034, -2.41763, -2.3638, -2.30869, -2.25204,
-2.19347, -2.13251, -2.06863, -2.00134, -1.93014, -1.85464, -1.77451,
-1.68954, -1.59962, -1.5048, -1.40522, -1.30117, -1.19306, -1.08138,
-0.96669, -0.84957, -0.73057, -0.61021, -0.48893, -0.36706, -0.24485,
-0.12246, 0.0, 0.122461, 0.244849, 0.367061, 0.488932, 0.610215,
0.730571, 0.849567, 0.966692, 1.081383, 1.193061, 1.301171, 1.405216,
1.504796, 1.599623, 1.689538, 1.774511, 1.854641, 1.930142, 2.001336,
2.068632, 2.132506, 2.193472, 2.252044, 2.308691, 2.363798, 2.417633,
2.470341, 2.52197, 2.572516, 2.621956, 2.670268, 2.717437, 2.763454,
2.808314, 2.852017, 2.894571, 2.935986, 2.976276, 3.015462, 3.053564,
3.090611, 3.126632, 3.1416};
/* Phases for 2 odd filter, 3 even bar response for small filter from -0.43 to 0.43 */
static double bar_phase2 [BAR_ENTRIES2] = { -3.1416,
-2.801149082, -2.792524659, -2.783437598, -2.773858169, -2.763864316, -2.753613695,
-2.743116001, -2.732306681, -2.721121223, -2.709503974, -2.697405066, -2.684777187,
-2.671573382, -2.657745581, -2.643243604, -2.628014402, -2.61200147, -2.595144345, -2.577378163,
-2.558633276, -2.538834912, -2.517902907, -2.495751511, -2.472289327, -2.447419396, -2.421039521,
-2.393042872, -2.363318825, -2.331752406, -2.298216642, -2.262552309, -2.224538371, -2.183863909,
-2.140110826, -2.092750398, -2.041151844, -1.984600117, -1.92232223, -1.853524781, -1.777448495,
-1.693446865, -1.601093519, -1.500314211, -1.391522873, -1.275720046, -1.154498225, -1.029910499,
-0.904205544, -0.779494823, -0.657454115, -0.539142619, -0.424964191, -0.314739594, -0.207835336,
-0.103302471, 0.0, 0.103302471, 0.207835336, 0.314739594, 0.424964191, 0.539142619,
0.657454115, 0.779494823, 0.904205544, 1.029910499, 1.154498225, 1.275720046, 1.391522873,
1.500314211, 1.601093519, 1.693446865, 1.777448495, 1.853524781, 1.92232223, 1.984600117,
2.041151844, 2.092750398, 2.140110826, 2.183863909, 2.224538371, 2.262552309, 2.298216642,
2.331752406, 2.363318825, 2.393042872, 2.421039521, 2.447419396, 2.472289327, 2.495751511,
2.517902907, 2.538834912, 2.558633276, 2.577378163, 2.595144345, 2.61200147, 2.628014402,
2.643243604, 2.657745581, 2.671573382,
2.684777187, 2.697405066, 2.709503974, 2.721121223, 2.732306681, 2.743116001,
2.753613695, 2.763864316, 2.773858169, 2.783437598, 2.792524659, 2.801149082,
3.1416};
static double bar_pos [BAR_ENTRIES2] = {
-1.0, -0.55, -0.54, -0.53, -0.52, -0.51, -0.5,
-0.49, -0.48, -0.47, -0.46, -0.45, -0.44,
/* 4 orientations starts here, index [13] */
-0.43, -0.42, -0.41, -0.4, -0.39, -0.38, -0.37,
-0.36, -0.35, -0.34, -0.33, -0.32, -0.31, -0.3,
-0.29, -0.28, -0.27, -0.26, -0.25, -0.24, -0.23,
-0.22, -0.21, -0.2, -0.19, -0.18, -0.17, -0.16,
-0.15, -0.14, -0.13, -0.12, -0.11, -0.1, -0.09,
-0.08, -0.07, -0.06, -0.05, -0.04, -0.03, -0.02,
-0.01, 0.0, 0.01, 0.02, 0.03, 0.04, 0.05,
0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12,
0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19,
0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26,
0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33,
0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4,
0.41, 0.42, 0.43,
0.44, 0.45, 0.46, 0.47, 0.48, 0.49,
0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 1.0};
float pos;
if (ODD_ANGLES == 4)
pos = (float) look_up (phase, BAR_ENTRIES4, bar_phase4, &bar_pos[13]);
else
pos = (float) look_up (phase, BAR_ENTRIES2, bar_phase2, bar_pos);
return (pos);
}
/*------------------------------------------------------------------------*/
/* find_step: given the phase, find the step position.
The return value is the edge position.
*/
float CFeature::find_step (double phase)
{
#define STEP_ENTRIES4 65
#define STEP_ENTRIES2 80
static double step_pos [STEP_ENTRIES2] = {
-1.0, -0.37, -0.36, -0.35, -0.34, -0.33,
/* for 4 orientations, start here at index [6] */
-0.32, -0.31, -0.3,
-0.29, -0.28, -0.27, -0.26, -0.25, -0.24, -0.23,
-0.22, -0.21, -0.2, -0.19, -0.18, -0.17, -0.16,
-0.15, -0.14, -0.13, -0.12, -0.11, -0.1, -0.09,
-0.08, -0.07, -0.06, -0.05, -0.04, -0.03, -0.02,
-0.01, 0.0, 0.01, 0.02, 0.03, 0.04, 0.05,
0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12,
0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19,
0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26,
0.27, 0.28, 0.29, 0.3, 0.31, 0.32,
/* end for 4 orientations */
0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 1.0};
/* Step phase: 0 to pi matches -0.31 to 0.31 */
static double step_phase4 [STEP_ENTRIES4] = {
0.0, 0.001827, 0.036898,
0.073045, 0.110287, 0.148642, 0.188126, 0.228755, 0.270539, 0.313488,
0.357606, 0.402895, 0.449352, 0.496966, 0.545724, 0.595607, 0.646589,
0.698639, 0.75172, 0.805788, 0.860796, 0.91669, 0.973412, 1.030898,
1.089083, 1.147896, 1.207264, 1.267112, 1.327362, 1.387936, 1.448754,
1.509735, 1.570796327, 1.631859, 1.692839, 1.753657, 1.814231, 1.874482,
1.934329, 1.993698, 2.052511, 2.110695, 2.168182, 2.224904, 2.280797,
2.335805, 2.389874, 2.442954, 2.495004, 2.545986, 2.595869, 2.644628,
2.692242, 2.738698, 2.783988, 2.828106, 2.871055, 2.91284, 2.953468,
2.992953, 3.031308, 3.06855, 3.104696, 3.139766, 3.1416};
static double step_phase2 [STEP_ENTRIES2] = { 0.0,
0.410874286, 0.418921308, 0.42932989, 0.440681409, 0.452713616, 0.465384309, 0.47870702,
0.492714419, 0.507448578, 0.522958036, 0.539296881, 0.556524524, 0.574705726, 0.593910734,
0.614215438, 0.635701501, 0.658456414, 0.682573444, 0.708151388, 0.73529409, 0.764109624,
0.794709041, 0.827204564, 0.861707105, 0.898322957, 0.937149556, 0.978270205, 1.021747747,
1.067617264, 1.115878073, 1.166485447, 1.219342788, 1.27429516, 1.331125276, 1.389553042,
1.44923954, 1.509795892, 1.570796788, 1.631797683, 1.692354033, 1.752040527, 1.810468289,
1.8672984, 1.922250767, 1.975108104, 2.025715477, 2.073976286, 2.119845808, 2.163323357,
2.20444402, 2.243270639, 2.279886519, 2.314389096, 2.346884667, 2.377484146, 2.40629976,
2.433442563, 2.459020637, 2.483137835, 2.505892969, 2.527379322, 2.547684419, 2.566889967,
2.585071931, 2.602300683, 2.618641203, 2.634153298, 2.648891831, 2.662906939, 2.676244282,
2.688945371, 2.701048208, 2.712588774, 2.723604716, 2.734144316, 2.744286265, 2.75417223,
2.764009427, 3.1416};
float pos;
double abs_phase;
abs_phase = fabs(phase);
if (ODD_ANGLES == 4)
pos = (float) look_up (abs_phase, STEP_ENTRIES4, step_phase4, &step_pos[6]);
else
pos = (float) look_up (abs_phase, STEP_ENTRIES2, step_phase2, step_pos);
if (phase < 0.0)
pos = -pos;
return (pos);
}
/*------------------------------------------------------------------------*/
/* Look up a value in a table and do linear interpolation.
Assume that values in x_table are monotonic. */
double CFeature::look_up( double x, int points, double *x_table, double *y_table)
{
int last, i, i_low, i_high;
double denom, y;
last = points - 1;
if (x_table[0] < x_table[last])
{
if (x <= x_table[0])
return (y_table[0]);
if (x >= x_table[last])
return (y_table[last]);
i_low = 0;
i_high = last;
while (i_high - i_low > 1)
{
if ((denom = x_table[i_high] - x_table[i_low]) == 0.0)
return (y_table [i_low]);
i = (int) (i_low + (i_high - i_low) * (x - x_table[i_low]) / denom);
if (i <= i_low)
i = i_low + 1;
if (i >= i_high)
i = i_high - 1;
if (x_table[i] < x)
i_low = i;
else
i_high = i;
}
}
else
{
if (x >= x_table[0])
return (y_table[0]);
if (x <= x_table[last])
return (y_table[last]);
i_low = last;
i_high = 0;
while (i_low - i_high > 1)
{
if ((denom = x_table[i_high] - x_table[i_low]) == 0.0)
return (y_table [i_low]);
i = (int) (i_high + (i_low - i_high) * (x_table[i_high] - x) /
denom);
if (i >= i_low)
i = i_low - 1;
if (i <= i_high)
i = i_high + 1;
if (x_table[i] < x)
i_low = i;
else
i_high = i;
}
}
y = y_table [i_low] + (y_table [i_high] - y_table [i_low]) *
(x - x_table[i_low]) / (x_table[i_high] - x_table[i_low]);
return (y);
}
/*------------------------------------------------------------------------*/
void CFeature::FindSecondary()
// assumes all steering has been done
{
}
/*------------------------------------------------------------------------*/
void CFeature::FindResidual()
{
}
/*------------------------------------------------------------------------*/
// read a feature from a file
/*
void CFeature::read(ifstream& inFile)
{
inFile >> m_x_rf >> m_y_rf >> m_diam_rf;
inFile >> m_strength >> m_angle >> m_x >> m_y;
int featureType;
inFile >> featureType;
m_type = (feature_type) featureType;
if (m_type != eEDGE)
{
inFile >> m_width;
}
if (m_type == eBLOB || m_type == eDARK_BLOB ||
m_type == eLIGHT_BLOB)
{
inFile >> m_length;
}
}
*/
/*------------------------------------------------------------------------*/
// Return true if the line thru (px + m_pos + delta, py) rotated m_degrees
// intersects the circle centered at (px,py) in one or two points.
bool Data_1D::DoesIntersect
( float px, float py,
float radius, // [in] of circle
float delta // [in] displacement
)
{
float qx, qy, pos;
/* This routine is used with m_answer, which defines an absolute position
for (m_x, m_y) but has m_pos == 0, since the center can be anywhere.
Thus we need to compute pos relative to the center (px, py).
*/
if (m_degrees < 180)
{
if (FL_ABS(m_cos_th) > FL_ABS(m_sin_th))
pos = (m_x - px) / m_cos_th;
else
pos = (m_y - py) / m_sin_th;
qx = (pos + delta) * m_cos_th;
qy = (pos + delta) * m_sin_th;
}
else
{
if (FL_ABS(m_cos_th) > FL_ABS(m_sin_th))
pos = -(m_x - px) / m_cos_th;
else
pos = -(m_y - py) / m_sin_th;
qx = -(pos + delta) * m_cos_th;
qy = -(pos + delta) * m_sin_th;
}
/* Q = (qx, qy) is a point on the line
The direction of the line is V = (-sin_th, cos_th)
The equation for the line is S(t) = Q + t V
The equation for the circle is (X-P).(X-P) = radius*radius
where P = (px, py)
The intersection of the line and circle is found by setting X = S(t)
and solving for t with the quadratic formula.
Two solutions exist if
(V . (Q-P))^2 + radius^2 > (Q-P) . (Q-P)
(Q-P) is the distance from the center of the circle to a point on the line
Its projection onto the line and the pythagorean theorem give a
geometric interpretation.
*/
float q_p = qx * qx + qy * qy;
float projection = -m_sin_th * qx + m_cos_th * qy;
return (projection * projection + radius * radius >= q_p);
}
/*------------------------------------------------------------------------*/
// Return true if the line thru (px + m_pos + delta, py) rotated m_degrees
// intersects the circle in one or two points.
bool Data_1D::FindIntersect
( float px, float py,
float radius, // [in] of circle
float delta, // [in] displacement
float *ax, float *ay, // [out] 1st intersection point
float *bx, float *by // [out] 2nd intersection point
)
{
float qx, qy, pos;
if (m_degrees < 180)
{
if (FL_ABS(m_cos_th) > FL_ABS(m_sin_th))
pos = (m_x - px) / m_cos_th;
else
pos = (m_y - py) / m_sin_th;
qx = (pos + delta) * m_cos_th;
qy = (pos + delta) * m_sin_th;
}
else
{
if (FL_ABS(m_cos_th) > FL_ABS(m_sin_th))
pos = -(m_x - px) / m_cos_th;
else
pos = -(m_y - py) / m_sin_th;
qx = -(pos + delta) * m_cos_th;
qy = -(pos + delta) * m_sin_th;
}
/* Q = (qx, qy) is a point on the line
The direction of the line is V = (-sin_th, cos_th)
The equation for the line is S(t) = Q + t V
The equation for the circle is (X-P).(X-P) = radius*radius
where P = (px, py)
The intersection of the line and circle is found by setting X = S(t)
and solving for t with the quadratic formula.
Two solutions exist if
radical = (V . (Q-P))^2 + radius^2 - (Q-P) . (Q-P) > 0
Solution is at t = -(V . (Q-P)) +/- sqrt(radical)
(Q-P) is the distance from the center of the circle to a point on the line
Its projection onto the line and the pythagorean theorem give a
geometric interpretation.
*/
float q_p = qx * qx + qy * qy;
float projection = -m_sin_th * qx + m_cos_th * qy;
float radical = (projection * projection + radius * radius - q_p);
if (radical < 0)
return (false);
float t1 = sqrt(radical);
float t2 = -t1 - projection;
t1 -= projection;
*ax = px + qx - t1 * m_sin_th;
*bx = px + qx - t2 * m_sin_th;
*ay = py + qy + t1 * m_cos_th;
*by = py + qy + t2 * m_cos_th;
return (true);
}
/*------------------------------------------------------------------------*/
void Data_1D::write
( ofstream& outFile,
bool names, // If names is true, give a text name,
// otherwise write numbers for feature_type;
int verbosity // 0 for main only, 1 for brief, 2 for longer ...
)
{
if (m_type == eNO_FEATURE || m_type == eNOT_AVAILABE)
outFile << '\t' << '\t' << '\t' << '\t';
else
{
outFile << m_degrees << '\t';
// 6/11/03 Strength output changed to pixels.
// Previously, it had been a percent.
outFile << (MAX_PIXEL-MIN_PIXEL)/2 * m_strength << '\t' << m_x << '\t' << m_y << '\t';
}
if (names)
{
outFile << Name(m_type) << '\t';
}
else
{
outFile << m_type << '\t';
}
outFile << m_pos << '\t';
if (m_type != eEDGE && m_type != eNO_FEATURE && m_type != eNOT_AVAILABE)
{
outFile << m_width;
}
outFile << '\t';
if (verbosity > 1)
{
if (m_type == eNO_FEATURE || m_type == eNOT_AVAILABE)
outFile << '\t' << '\t' << '\t' << '\t';
else
{
outFile << m_corrOdd << '\t' << m_corrEven << '\t' << m_corrBigOdd <<
'\t' << m_corrBigEven << '\t';
}
}
}
/*------------------------------------------------------------------------*/
void Data_2D::write
( ofstream& outFile,
bool names, // If names is true, give a text name,
// otherwise write numbers for feature_type;
int verbosity // 0 for main only, 1 for brief, 2 for longer ...
)
{
if (m_type == eNO_FEATURE || m_type == eNOT_AVAILABE)
outFile << '\t' << '\t';
else
{
outFile << m_x << '\t' << m_y << '\t';
}
if (names)
{
outFile << Name(m_type) << '\t';
}
else
{
outFile << m_type << '\t';
}
}
/*------------------------------------------------------------------------*/
// return a string describing the feature type.
const char* Data_1D::Name( feature_type featureType)
{
switch (featureType)
{
case eNO_FEATURE:
default:
return ("None");
case eEDGE:
return ("Edge");
case eWHITE_LINE:
return ("LiteBar");
case eBLACK_LINE:
return ("DarkBar");
}
}
/*------------------------------------------------------------------------*/
// return a string describing the feature type.
const char* Data_2D::Name( feature_type featureType)
{
switch (featureType)
{
case eNO_FEATURE:
default:
return ("None");
case eEDGE:
return ("Edge");
case eWHITE_LINE:
return ("LiteBar");
case eBLACK_LINE:
return ("DarkBar");
case eCORNER:
return ("Corner");
case eDARK_CORNER:
return ("Dark Corner");
case eLIGHT_CORNER:
return ("Light Corner");
case eBLOB:
return ("Blob");
case eDARK_BLOB:
return ("Dark Blob");
case eLIGHT_BLOB:
return ("Light Blob");
}
}