iffeature.cpp

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

// 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 =