iffeature.cpp

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

    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");
    }
}