v2.cpp

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

// V2.cpp: Intermediate level methods for class CFeature

// The previous processing performed in class CFeature (file IFFeature.cpp)
// roughly corresponds to primary visual cortex, or area V1.
// The routines in this file can be thought of as doing the next higher level processing.
// No attempt is made to emulate actual biological mechanisms.


/*
** Copyright (C) 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"

const char *featureName[eFEATURE_TYPE_SIZE] =
{	"none", "NONE",	"Edge",
	"BAR", "bar",	  // Uppercase for dark; lower case for light
	"Corner",	"CORNER",	"corner",
	"Blob",	"BLOB",	"blob",	"END",	"end"
};

/*------------------------------------------------------------------------*/
// Fill in CFeature::m_pNbr with pointers to adjacent features
bool CFeature::IsNeighbor(CFeature *other)
{
	int i, far_index;
	float dist_sqr, farthest, dSqr;
	float xyd1 = m_x_rf * m_x_rf
			   + m_y_rf * m_y_rf
			   - m_diam_rf * m_diam_rf / 4;
	float twiceX1 = m_x_rf * 2;
	float twiceY1 = m_y_rf * 2;

	// Neighbor if the distance between centers is less than the sum of radii.
	float radius = other->m_diam_rf / 2;
	if (xyd1 - other->m_x_rf * twiceX1 + other->m_x_rf * other->m_x_rf
		     - other->m_y_rf * twiceY1 + other->m_y_rf * other->m_y_rf
		< (m_diam_rf + radius) * radius )
	{  // these are neighbors
		if (m_pNbr[NEIGHBORS-1] == 0)
		{	// room for more 
			for (i = 0; i < NEIGHBORS; i++)
			{
				if (m_pNbr[i] == 0)
				{
					m_pNbr[i] = other;
					break;
				}
			}
		}
		else
		{  // replace the farthest neighbor
			dist_sqr = (m_x_rf - other->m_x_rf) * (m_x_rf - other->m_x_rf) 
			         + (m_y_rf - other->m_y_rf) * (m_y_rf - other->m_y_rf);
			farthest = 0;
			for (i = 0; i < NEIGHBORS; i++)
			{
				if (m_pNbr[i] == 0)
				{
					m_pNbr[i] = other;
					far_index = -1;
					break;
				}
				dSqr = (m_x_rf - m_pNbr[i]->m_x_rf) * (m_x_rf - m_pNbr[i]->m_x_rf) 
			         + (m_y_rf - m_pNbr[i]->m_y_rf) * (m_y_rf - m_pNbr[i]->m_y_rf);
				if (dSqr == 0)
				{   // do not count a field with the same center as a neighbor
					m_pNbr[i] = other;
					far_index = -1;
					break;
				}
				if (dSqr > farthest)
				{
					farthest = dSqr;
					far_index = i;
				}
			}
			if (far_index >= 0 && dist_sqr < farthest)
				m_pNbr[far_index] = other;
		}

		if (other->m_pNbr[NEIGHBORS-1] == 0)
		{	// room for more 
			for (i = 0; i < NEIGHBORS; i++)
			{
				if (other->m_pNbr[i] == 0)
				{
					other->m_pNbr[i] = this;
					break;
				}
			}
		}
		else
		{  // replace the farthest neighbor
			dist_sqr = (m_x_rf - other->m_x_rf) * (m_x_rf - other->m_x_rf) 
			         + (m_y_rf - other->m_y_rf) * (m_y_rf - other->m_y_rf);
			farthest = 0;
			for (i = 0; i < NEIGHBORS; i++)
			{
				if (other->m_pNbr[i] == 0)
				{
					other->m_pNbr[i] = this;
					far_index = -1;
					break;
				}
				dSqr = (m_x_rf - other->m_pNbr[i]->m_x_rf) * (m_x_rf - other->m_pNbr[i]->m_x_rf) 
			         + (m_y_rf - other->m_pNbr[i]->m_y_rf) * (m_y_rf - other->m_pNbr[i]->m_y_rf);
				if (dSqr == 0)
				{   // do not count a field with the same center as a neighbor
					other->m_pNbr[i] = this;
					far_index = -1;
					break;
				}
				if (dSqr > farthest)
				{
					farthest = dSqr;
					far_index = i;
				}
			}
			if (far_index >= 0 && dist_sqr < farthest)
				other->m_pNbr[far_index] = this;
		}

		return true;
	}
	return false;
}
/*------------------------------------------------------------------------*/
// Lateral antagonism: increase or decrease the strength of neighbors.
void CFeature::LateralAntagonism()
{
	/* 
	If an edge passes though a neighbor receptive field in the same direction
	 that neighbor has its strength increased.
	Else, if the normal to this edge passes though a neighbor, 
	 that neighbor has its strength decreased.
	The amount of increase or decrease is proportional to the strength of this cell.

	*/

	PROFILE("LateralAntagonism");
	if (m_corner.m_type < eEDGE || m_corner.m_type >= eBLOB)
		return;

    float vx = m_main.m_cos_th;
    float vy = m_main.m_sin_th;
	if (m_angle > PI_1)
	{
		vx = -vx;
		vy = -vy;
	}
	float half_width = m_main.m_width / 2;
    float qy  = m_y_rf - (m_main.m_pos + half_width) * vx;
    float qx  = m_x_rf + (m_main.m_pos + half_width) * vy;
    float qy2 = m_y_rf - (m_main.m_pos - half_width) * vx;
    float qx2 = m_x_rf + (m_main.m_pos - half_width) * vy;
	int i;

	if (m_corner.m_type == eCORNER)
	{  // Use neighbors to decide if light or dark corner
		float light = 0;
		float dark = 0;
		for (i = 0; i < NEIGHBORS; i++)
		{
			if (m_pNbr[i] == 0)
				continue;
			float wx = m_pNbr[i]->m_corner.m_x - m_corner.m_x;
			float wy = m_pNbr[i]->m_corner.m_y - m_corner.m_y;
			float wlength = (float) sqrt(wx*wx + wy*wy);
			wx /= wlength;
			wy /= wlength;
			// dot product of unit vector in corner direction with
			// unit vector from corner to neighboring feature
			float cos_phi = vx * wx + vy * wy;
			float normal = -vy * wx + vx * wy;
			float diff = m_pNbr[i]->m_main.m_degrees - m_main.m_degrees;
			while (diff >= 360)
				diff -= 360;
			while (diff < 0)
				diff += 360;
			if ((cos_phi >= ROOT_2 / 2) ||
				(normal >=  ROOT_2 / 2))
			{  // zone for a dark corner
				if (135 < diff && diff < 315)
					dark -= m_pNbr[i]->m_main.m_strength;
				else
					dark += m_pNbr[i]->m_main.m_strength;
			}
			else
			{  // zone for a light corner
				if (135 < diff && diff < 315)
					light += m_pNbr[i]->m_main.m_strength;
				else
					light -= m_pNbr[i]->m_main.m_strength;
			}
		}  // loop on neighbors
		if (light > dark)
		{
			m_corner.m_type = eLIGHT_CORNER;
		}
		else if (dark > light)
		{
			m_corner.m_type = eDARK_CORNER;
		}
	}
/*
	// new method needs to be debugged
	for (i = 0; i < NEIGHBORS; i++)
	{
		if (m_pNbr[i] == 0)
			continue;
		float wx = m_pNbr[i]->m_corner.m_x - m_corner.m_x;
		float wy = m_pNbr[i]->m_corner.m_y - m_corner.m_y;
		float wlength = (float) sqrt(wx*wx + wy*wy);
		wx /= wlength;
		wy /= wlength;
		// dot product of unit vector in corner direction with
		// unit vector from corner to neighboring feature
		float cos_phi = vx * wx + vy * wy;
		float normal = -vy * wx + vx * wy;
		if (cos_phi >= ROOT_2 / 2)
		{	// feature point is in the same direction
			if (m_corner.m_type != eLIGHT_CORNER)
				m_pNbr[i]->m_lateral += m_main.m_strength;
		}
		else if (cos_phi <= -ROOT_2 / 2)
		{	// neighbor is in opposite direction
			if (m_corner.m_type != eDARK_CORNER)
				m_pNbr[i]->m_lateral += m_main.m_strength;
		}
		else if (normal >= ROOT_2 / 2)
		{	// neighbor is in the normal direction
			if (m_corner.m_type == eDARK_CORNER ||
				m_corner.m_type == eCORNER)
				m_pNbr[i]->m_lateral += m_main.m_strength;
			else if (m_corner.m_type == eWHITE_LINE ||
				m_corner.m_type == eBLACK_LINE ||
				m_corner.m_type == eEDGE)
				m_pNbr[i]->m_lateral -= m_main.m_strength;
		}
		else if (normal <= -ROOT_2 / 2)
		{	// neighbor in the negative normal direction
			if (m_corner.m_type == eLIGHT_CORNER ||
				m_corner.m_type == eCORNER)
				m_pNbr[i]->m_lateral += m_main.m_strength;
			else if (m_corner.m_type == eWHITE_LINE ||
				m_corner.m_type == eBLACK_LINE ||
				m_corner.m_type == eEDGE)
				m_pNbr[i]->m_lateral -= m_main.m_strength;
		}
	}
*/
/*
	If any part of this edge passes though a neighbor receptive field,
	 that neighbor has its strength increased.
	Else, if the normal to this edge passes though a neighbor, 
	 that neighbor has its strength decreased.
	The amount of increase or decrease is proportional to the strength of this cell.

    Consider a line passing thru the point (qx, qy) in the direction of the 
	unit vector (vx, vy) and a circle of radius r centered at (x, y).
	They intersect in two points if
	r >  vx * (qy - y) - vy * (qx - x).
	The normal line thru (qx, qy) intersects the circle if
	r > -vy * (qy - y) - vx * (qx - x).
*/
	// old method
	for (i = 0; i < NEIGHBORS; i++)
	{
		if (m_pNbr[i] == 0)
			continue;
		if (m_pNbr[i]->m_diam_rf /2 > -vy * (qy  - m_pNbr[i]->m_y_rf)
									 + vx * (qx  - m_pNbr[i]->m_x_rf)  ||
			m_pNbr[i]->m_diam_rf /2 > -vy * (qy2 - m_pNbr[i]->m_y_rf)
									 + vy * (qx2 - m_pNbr[i]->m_x_rf) )
			m_pNbr[i]->m_lateral += m_main.m_strength;
		else if (m_pNbr[i]->m_diam_rf / 2 > vx * (qy - m_pNbr[i]->m_y_rf)
										  - vy * (qx - m_pNbr[i]->m_x_rf))
			m_pNbr[i]->m_lateral -= m_main.m_strength;
	}
}
/*------------------------------------------------------------------------*/
/* Find Interest Points.
   Edges are 1 dimensional features.  
     They have nothing in the direction perpendicular to the main orientation.
   Interest points are 2 dimensional.
      They have a significant component perpendicular to the main orientation.
	  These points are most significant for object recognition.
	  Interest points can be subclassified into:
		Corners
		X crossings
		T junctions
		End-stopped bars
		Blobs
    Slightly curved lines can morph into rounded corners as the curvature increases.
	Thus it is hard to make a sharp distinction between Straight edges and interest points.
	Interest points correspond to vertices and are the terminators for edges.
	In computer graphics, a vertex must lie on three edges.
*/
void CFeature::corners()
{
	if (m_90.m_strength < CORNER_THRESH * m_main.m_strength)
		return;  // An edge

	switch(m_main.m_type) 
	{
	case eEDGE:
		switch(m_90.m_type) 
		{
		case eEDGE:
			// Need to find out what this really look like.  Maybe a double corner?
			m_corner.m_type = eCORNER;
			break;
		case eWHITE_LINE:
			m_corner.m_type = eLIGHT_CORNER;
			break;
		case eBLACK_LINE:
			m_corner.m_type = eDARK_CORNER;
			break;
		default:
			return;
		}
		break;

	case eWHITE_LINE:
		switch(m_90.m_type) 
		{
		case eEDGE:
			m_corner.m_type = eSTOPPED_BAR_LIGHT;
			break;
		case eWHITE_LINE:
			m_corner.m_type = eLIGHT_BLOB;  // or an X crossing
			break;
		case eBLACK_LINE:
			// Need to find out what this really look like.  Maybe a double corner?
			m_corner.m_type = eCORNER;
			break;
		default:
			return;
		}
		break;

	case eBLACK_LINE:
		switch(m_90.m_type) 
		{
		case eEDGE:
			m_corner.m_type = eSTOPPED_BAR_DARK;
			break;
		case eWHITE_LINE:
			// Need to find out what this really look like.  Maybe a double corner?
			m_corner.m_type = eCORNER;
			break;
		case eBLACK_LINE:
			m_corner.m_type = eDARK_BLOB;  // or an X crossing
			break;
		default:
			return;
		}
		break;

	default:
		return;
	}

	// If we haven't returned yet, it is some kind of interest point.
	// Set its location to the intersection of the two edges.
	float x1 = m_main.m_x - m_x_rf;
	float y1 = m_main.m_y - m_y_rf;
	float x2 = m_90.m_x - m_x_rf;
	float y2 = m_90.m_y - m_y_rf;
	/* denom can be zero if the lines from the origin to (x1,y1)
	 and to (y2, y2) are parallel; but we know that these lines are perpendicular. 
	 It could also be zero if (x1, y1) or (x2, y2) is at the origin.
	*/
	float denom = x1 * y2 - x2 * y1;
	if (FL_ABS(denom) > EPSILON)
	{
		m_corner.m_x = (y2 * m_main.m_pos * m_main.m_pos -
						y1 * m_90.m_pos * m_90.m_pos) / denom + m_x_rf;