www.pudn.com > imageFeatures_Ver3_0_source.zip > V2.cpp


// 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; 
		m_corner.m_y = (x1 * m_90.m_pos * m_90.m_pos - 
						x2 * m_main.m_pos * m_main.m_pos) / denom + m_y_rf; 
	} 
	else if (FL_ABS(m_main.m_pos) < EPSILON) 
	{	// main feature is at origin 
		m_corner.m_x = m_90.m_x; 
		m_corner.m_y = m_90.m_y; 
	} 
	// else perpendicular feature is at origin. 
	// But we defaulted to using (m_main.m_x, m_main.m_y), which works. 
} 
 
/*------------------------------------------------------------------------*/ 
// write text feature information to a file. 
 
void CFeature::write 
(   ofstream& outFile,  
    bool names,     // If names is true, give a text name, 
                    // otherwise write numbers for feature_type; 
    int verbosity   // 0 for brief, 1 for longer ... 
)   
{ 
    if (m_corner.m_type == eNO_FEATURE || m_corner.m_type == eNOT_AVAILABE) 
    { 
        outFile << endl; 
        return; 
    } 
 
    // set to show 1 decimal place 
    int oldPrecision = outFile.precision(1); 
    outFile.setf( ios::fixed, ios::floatfield ); 
    outFile << m_x_rf << '\t' << m_y_rf << '\t' << m_diam_rf << '\t' 
		<< m_lateral << '\t'; 
//	int type_index = (int) m_corner.m_type; 
//	outFile << featureName[ type_index] << '\t'; 
    m_corner.write( outFile, names, verbosity ); 
    m_main.write( outFile, names, verbosity ); 
    if (verbosity > 0) 
    { 
        // Write perpendicular features 
        m_90.write( outFile, names, verbosity ); 
/*        // Write +45 degree features 
        m_45.write( outFile, names, verbosity ); 
        // Write -45 degree features 
        m_135.write( outFile, names, verbosity );  */ 
    } 
    outFile << endl; 
    outFile.precision(oldPrecision); // restore previous precision 
} 
/*------------------------------------------------------------------------*/ 
// Compute distance from feature to base position 
float CFeature::PositionDifference 
(   Data_1D& base,     // the answer; true position of feature  
    bool fromRfCenter  // if true, use center of receptive field. 
                       // if false, use location of found feature  
) 
{ 
    double distance = 0.0; 
    if (m_corner.m_type == eNO_FEATURE || m_corner.m_type == eNOT_AVAILABE) 
        return (float) distance; 
 
    // intersection of base line with perpendicular thru m_main.m_x,y 
    double x = base.m_x; 
    double y = base.m_y; 
    double degrees = base.m_degrees; 
    double xFrom = m_main.m_x; 
    double yFrom = m_main.m_y; 
    if (fromRfCenter) 
    { 
        xFrom = m_x_rf; 
        yFrom = m_y_rf; 
    } 
    switch (base.m_type) 
    { 
    case eNO_FEATURE: 
    case eNOT_AVAILABE: 
        // distance is meaningless 
        break; 
 
    case eEDGE: 
    case eWHITE_LINE: 
    case eBLACK_LINE: 
        // TO DO: case of mismatched edges and bars 
        while (degrees > 180) 
            degrees -= 180; 
        while (degrees < 0) 
            degrees += 180; 
        if (degrees > 90 - EPSILON && degrees < 90 + EPSILON) 
            x = xFrom; 
        else if (degrees > 180 - EPSILON || degrees < EPSILON) 
            y = yFrom; 
        else 
        { 
            double slope = tan( degrees * PI / 180.0); 
            x = (slope / (slope * slope + 1.0)) * ( yFrom - base.m_y + 
                xFrom / slope + base.m_x * slope); 
            y = slope * (x - base.m_x) + base.m_y; 
        } 
        distance = sqrt( (x - xFrom) * (x - xFrom) + (y - yFrom) * (y - yFrom)); 
        if (xFrom - x + yFrom - y < 0) 
            distance = -distance; 
 
        // if matched an edge to a bar, report the distance to the edge of the bar 
        if (m_main.m_type == eEDGE &&  
           (base.m_type == eBLACK_LINE || base.m_type == eWHITE_LINE)) 
        { 
            double half_width =  base.m_width * 0.5; 
            distance += fabs( distance - half_width) < fabs( distance + half_width)? 
                -half_width : half_width; 
        } 
        if (base.m_type == eEDGE &&  
           (m_main.m_type == eBLACK_LINE || m_main.m_type == eWHITE_LINE)) 
        { 
            double cosine =  
                fabs( cos( (base.m_degrees - m_main.m_degrees) * PI / 180.0)); 
            if (cosine > EPSILON) 
            { 
                double half_width =  base.m_width * 0.5 / cosine; 
                distance += fabs( distance - half_width) <  
                    fabs( distance + half_width)? 
                    -half_width : half_width; 
            } 
        } 
        break; 
 
    case eCORNER: 
    case eDARK_CORNER: 
    case eLIGHT_CORNER: 
    case eBLOB: 
    case eDARK_BLOB: 
    case eLIGHT_BLOB: 
    default: 
        distance = sqrt( (x - xFrom) * (x - xFrom) + (y - yFrom) * (y - yFrom)); 
        if (xFrom - x + yFrom - y < 0) 
            distance = -distance; 
        break; 
    } 
    return (float) distance; 
 
} 
/*------------------------------------------------------------------------*/ 
// Compute error in type identification  
// Returned: difference in type values. 
int CFeature::TypeDifference( Data_1D& base ) 
{ 
    enum feature_type answer = base.m_type; 
    float half_width =  base.m_width * 0.5; 
	bool edge1seen = base.DoesIntersect(m_x_rf, m_y_rf, m_diam_rf/2, half_width); 
	bool edge2seen = edge1seen;			 
    if (answer == eBLACK_LINE || answer == eWHITE_LINE) 
		edge2seen = base.DoesIntersect(m_x_rf, m_y_rf, m_diam_rf/2, -half_width);			 
    // if the feature is not in the small receptive field, no feature. 
	if (!edge1seen && !edge2seen)			 
        answer = eNO_FEATURE; 
    // if the feature is a bar, can we see both edges? 
    if ((answer == eBLACK_LINE || answer == eWHITE_LINE) && 
		(!edge1seen || !edge2seen)) 
        answer = eEDGE;  // can't see far bar edge 
    return (int) (m_corner.m_type - answer); 
} 
/*------------------------------------------------------------------------*/ 
void CFeature::FindPlotInfo( Outline& shape, bool useMain ) 
{   // replaces CadData 
	Data_1D& oneD = useMain? m_main : m_90; 
    float radius = m_diam_rf / 2; /* effective radius of smaller filter */ 
    if (oneD.m_type == eNO_FEATURE || oneD.m_type == eNOT_AVAILABE) 
		return; 
    shape.color = EDGE_COLOR;  
	float half_width =  oneD.m_width * 0.5; 
	float p1x, p1y, p2x, p2y, p3x, p3y, p4x, p4y; 
	bool edge1seen = oneD.FindIntersect(m_x_rf, m_y_rf, radius, half_width, 
		&p1x, &p1y, &p2x, &p2y); 
	bool edge2seen = edge1seen;			 
    if (oneD.m_type == eBLACK_LINE || oneD.m_type == eWHITE_LINE) 
	{ 
		edge2seen = oneD.FindIntersect(m_x_rf, m_y_rf, radius, -half_width, 
		&p4x, &p4y, &p3x, &p3y); 
		shape.color = oneD.m_type == eWHITE_LINE? 
		LIGHT_BAR_COLOR: DARK_BAR_COLOR; 
	} 
	if (Is2D()) 
	{ 
        shape.color = CORNER_COLOR; 
		if (oneD.m_type == eEDGE) 
		{ 
			p2x = m_corner.m_x; 
			p2y = m_corner.m_y; 
		} 
	} 
	if (!edge1seen && !edge2seen)			 
        return;  // should not happen 
 
    // Make a quadrilateral describing the bar or line 
    POINT vertex;  
    vertex.x = Round(p1x);        
    vertex.y = Round(p1y); 
    shape.segment.Add( vertex ); 
    POINT firstPoint = vertex; 
    vertex.x = Round(p2x);        
    vertex.y = Round(p2y); 
    shape.segment.Add( vertex ); 
	if (oneD.m_type != eEDGE) 
	{ 
		vertex.x = Round(p3x);        
		vertex.y = Round(p3y); 
		shape.segment.Add( vertex ); 
		vertex.x = Round(p4x);        
		vertex.y = Round(p4y); 
		shape.segment.Add( vertex ); 
        shape.segment.Add( firstPoint ); 
	} 
    shape.thickness = 1; 
 
} 
/*------------------------------------------------------------------------*/ 
void CFeature::CadData( Outline& shape, bool useMain ) 
{ 
/* cad_info: write information that is useful for making a drawing of 
   the features that were located. */ 
/* This routine added 5/11/95 by Tyler C. Folsom */ 
 
//  float m_x_rf, m_y_rf,  /* center point of the filter */ 
//  float m_pos,           /* position of feature relative to center */ 
//  float m_angle,         /* angle of feature in radians */ 
//  int   m_type,          /* feature type */ 
//  float m_bar_width,     /* width of a bar */ 
 
	if (!useMain) 
		return;  // corner drawing is not yet working 
	Data_1D& oneD = useMain? m_main : m_90; 
    float radius = m_diam_rf / 2; /* effective radius of smaller filter */ 
 
    float sin_a, cos_a, half_bar; 
  /* intersection of circle with line from center to feature location. */ 
    float cx, cy; 
  /* end points of the feature in the smaller filter */ 
    float px, py, qx, qy; 
    /* centerline of bar */ 
    float pcx, pcy, qcx, qcy; 
  /* angle by which (cx, cy) is rotated into (px, py), (qx, qy). */ 
    double alpha; 
    POINT vertex;  
	float sin_th = oneD.m_sin_th; 
	float cos_th = oneD.m_cos_th; 
	/*  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 clockwise. 
		Thus for rotation, the origin is at (0,R) and the angle is positive. 
		x = -radius * sin_th; 
		y = radius * cos_th; 
		(cx, cy) is normal to (x,y) 
	*/ 
 
  /* intersection of circle with line from center to feature location. */ 
	if (oneD.m_degrees < 180) 
	{ 
		cx =  radius * cos_th; 
		cy =  radius * sin_th; 
	} 
	else 
	{ 
		cx = -radius * cos_th; 
		cy = -radius * sin_th; 
	} 
    if (oneD.m_type == eBLACK_LINE || oneD.m_type == eWHITE_LINE) 
    { 
        half_bar = (float) 0.5 * oneD.m_width; 
        shape.color = oneD.m_type == eWHITE_LINE? 
			LIGHT_BAR_COLOR: DARK_BAR_COLOR; 
        if (FL_ABS(oneD.m_pos) < half_bar) 
        {   /* center line is the longest segment */ 
            cos_a = oneD.m_pos / radius; 
            sin_a = (float) sqrt(1.0 - cos_a * cos_a); 
            /* centerline of bar: (cx, cy) rotated +/- alpha */ 
            pcx = m_x_rf + cx * cos_a + cy * sin_a; 
            pcy = m_y_rf - cx * sin_a + cy * cos_a; 
            qcx = m_x_rf + cx * cos_a - cy * sin_a; 
            qcy = m_y_rf + cx * sin_a + cy * cos_a;  
            /* inner edge of bar */ 
            px = pcx - half_bar * cos_th; 
            py = pcy - half_bar * sin_th; 
            qx = qcx - half_bar * cos_th; 
            qy = qcy - half_bar * sin_th; 
        } 
        else 
        {   /* inner edge is the longest segment */ 
            alpha = acos ((double) ((oneD.m_pos - half_bar) / radius)); 
            sin_a = (float) sin(alpha); 
            cos_a = (float) cos(alpha); 
            /* inner edge of bar */ 
            px = m_x_rf + cx * cos_a + cy * sin_a; 
            py = m_y_rf - cx * sin_a + cy * cos_a; 
            qx = m_x_rf + cx * cos_a - cy * sin_a; 
            qy = m_y_rf + cx * sin_a + cy * cos_a; 
            /* centerline of bar */ 
            pcx = px + half_bar * cos_th; 
            pcy = py + half_bar * sin_th; 
            qcx = qx + half_bar * cos_th; 
            qcy = qy + half_bar * sin_th; 
        } 
        // Make a quadrilateral describing the bar 
        vertex.x = Round(px);        
        vertex.y = Round(py); 
        POINT firstPoint = vertex; 
        shape.segment.Add( vertex ); 
        vertex.x = Round(qx);        
        vertex.y = Round(qy); 
        shape.segment.Add( vertex ); 
 
        /* outer edge of bar */ 
        px = pcx + half_bar * cos_th; 
        py = pcy + half_bar * sin_th; 
        qx = qcx + half_bar * cos_th; 
        qy = qcy + half_bar * sin_th; 
 
        vertex.x = Round(qx);        
        vertex.y = Round(qy); 
        shape.segment.Add( vertex ); 
        vertex.x = Round(px);        
        vertex.y = Round(py); 
        shape.segment.Add( vertex ); 
        shape.segment.Add( firstPoint ); 
 
    } 
    else   /* a step edge */ 
    {   /* compute the points where the edge meets the boundary of the smaller  
           receptive field */ 
        shape.color = EDGE_COLOR;  
        alpha = acos ((double) (oneD.m_pos / radius)); 
        sin_a = (float) sin(alpha); 
        cos_a = (float) cos(alpha); 
        // Points on the edge are (cx, cy) rotated by +/- alpha. 
        px = m_x_rf + cx * cos_a + cy * sin_a; 
        py = m_y_rf - cx * sin_a + cy * cos_a; 
        qx = m_x_rf + cx * cos_a - cy * sin_a; 
        qy = m_y_rf + cx * sin_a + cy * cos_a; 
 
        vertex.x = Round(px);        
        vertex.y = Round(py); 
        shape.segment.Add( vertex ); 
        vertex.x = Round(qx);        
        vertex.y = Round(qy); 
        shape.segment.Add( vertex ); 
    } 
/* 
	if (Is2D()) 
	{ 
        shape.color = CORNER_COLOR;  
	} 
*/ 
    if (oneD.m_strength < 0.5) 
    { 
        shape.thickness = 1; 
    } 
    else if (oneD.m_strength < 1.0) 
    { 
        shape.thickness = 2; 
    } 
    else if (oneD.m_strength < 1.5) 
    { 
        shape.thickness = 3; 
    } 
    else 
    { 
        shape.thickness = 4; 
    } 
} 
/*------------------------------------------------------------------------*/ 
// draw a circle corresponding to the small filter. 
BOOL CFeature::DrawCircle(HDC hDeviceContext) 
{ 
    return  
    (  
        Arc( hDeviceContext,  
            m_x_rf - m_diam_rf/2, // x Up Left 
            m_y_rf - m_diam_rf/2, // y Up Left 
            m_x_rf + m_diam_rf/2, // x Low Right 
            m_y_rf + m_diam_rf/2, // y Low Right 
            m_x_rf,                          // x start 
            m_y_rf + m_diam_rf/2, // y start 
            m_x_rf,                          // x end 
            m_y_rf + m_diam_rf/2  // y end 
            ) 
    ); 
 
}