iffeature.cpp

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

                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]);