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


/* 
** 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. 
*/ 
/* 
 
ProcessFeatures.cpp 
 
 
The user opens an image, which is displayed in its own window.   
He / She then selects Process / Find Features.  This establishes 
a linkage between the input image (the opened image most  
recently active), the grid file (default name is grid.txt),  
the text output file (same name as the image file with an extension 
of ".txt") and the graphical output file (same name as the image  
file, but with a "$" inserted before the extension) 
If the grid file has not previously been read, it is opened and read. 
 
For each filter location 
{ 
	Create new filters of this size if needed. 
	Acquire a subimage for the bigger filter 
	Apply filters to the big subimage. 
	Acquire a subimage for the smaller filter 
	Apply filters to the small subimage. 
	Use the filter results to produce feature characteristics 
	Add the feature to the text output file. 
} 
Save and display the text output file. 
Generate information for graphics display. 
Display the graphics output. 
*/ 
 
#include "stdafx.h" 
#include "ProcessFeatures.h" 
 
#include "quad_dis.h"  // for FL_ABS() 
#include "typeinfo.h" 
 
#include "profile.h" 
 
// Initialize static data member 
CMap CProcessFeatures::m_filters; 
 
/*------------------------------------------------------------------------*/ 
// This does the main work. 
bool CProcessFeatures::Process( grid_type grid) 
{ 
	PROFILE("CProcessFeatures::Process"); 
 
#ifdef TRACE_ME 
	ofstream traceFile( "Trace.txt", ios::out ); 
	traceFile << "Starting Process" << endl; 
	traceFile.close(); 
#endif // TRACE_ME 
    bool success = false; 
    // Read a file and construct CLocations at which to filter 
    switch (grid) 
    { 
    case eTest1D: 
        success = Make1DGrid(); 
        break; 
    case eTest2D: 
        success = Make2DGrid(); 
        break; 
    default: 
    case eMakeGrid: 
        success = MakeGrid(); 
        break; 
    case eReadGrid: 
        success = ReadGrid(); 
    } 
    if (!success) 
        return false;   // can't get the sampling locations. 
 
    m_features.RemoveAll();  // get rid of any previous features. 
    // iterate on m_locs   
    for ( int n = 0; n < m_locs.GetSize(); n++)  
    { 
        CLocation location = m_locs[n]; 
        if (Filter( location )) // filter at this location 
        { 
            Interpret( location ); // find feature at this location 
        } 
    } 
 
	FindNeighbors(); 
 
    // iterate on m_features   
    for ( n = 0; n < m_features.GetSize(); n++)  
    {   // facilitate or inhibit based on neighbors. 
        m_features[n].LateralAntagonism(); 
    } 
    for ( n = 0; n < m_features.GetSize(); n++)  
    {   // remove features below threshold. 
		if (m_features[n].GetStrength()  
			+ m_lateralFactor * m_features[n].GetAntagonism() 
			< m_threshold) 
			m_features.RemoveAt(n); 
	} 
 
#ifdef TRACE_ME 
	traceFile << "Completed Process" << endl; 
	traceFile.close(); 
#endif // TRACE_ME 
 
    return true; 
} 
/*------------------------------------------------------------------------*/ 
// If using a templated version of the image, would declare as: 
// templateCProcessFeatures::CProcessFeatures( ) 
CProcessFeatures::CProcessFeatures( ) 
{ 
    m_gridFileName = "grid.txt"; 
    m_threshold = (float) 0.20; 
	m_fixedDirection = false; 
	// This is the minimum overlap that will cover the whole image 
	m_overlap = (float) (PI_2 / (3 * ROOT_3)); // 1.209 
	m_filterDiam = 8; 
	m_displayTestImages = false; 
	m_testImageContrast = 20; 
	m_testImageNoise = 20; 
    m_viewEdgesOnly = false; 
    m_viewGrid = false; 
	m_subSampleDiam = DEFAULT_SUBSAMPLE_DIAM; 
	m_lateralFactor = DEFAULT_LATERAL; 
    SetAnswer(); 
} 
//---------------------------------------------------------------------------------------- 
// constructor for creating the answer for test images  
void CProcessFeatures::SetAnswer 
(   feature_type type, 
    float x,    // center point of feature 
    float y,    // center point of feature 
    float degrees, 
    float strength, 
    float width 
) 
{ 
    // Originally a right pointing horizontal line was 0 degrees. 
    // It has been redefined so that 0 degrees points down.  7/1/03 
    // The 0 degree kernel had been horizontal; it is now vertical.  5/20/03. 
    // Test images use the old convention. 
    m_answer.m_degrees = degrees + 90; 
    while (m_answer.m_degrees >= 360) 
        m_answer.m_degrees -= 360; 
    m_answer.m_corrEven = 0; 
    m_answer.m_corrOdd = 0; 
    m_answer.m_corrBigEven = 0; 
    m_answer.m_corrBigOdd = 0; 
    m_answer.m_strength = strength; 
    m_answer.m_x = x; 
    m_answer.m_y = y; 
    m_answer.m_pos = 0;  // depends on where filter is 
    m_answer.m_type = type; 
    m_answer.m_width = width; 
	m_answer.m_cos_th = cos(m_answer.m_degrees * PI / 180.0); 
	m_answer.m_sin_th = sin(m_answer.m_degrees * PI / 180.0); 
} 
/*------------------------------------------------------------------------*/ 
// Construct CLocations at which to filter an image. 
// Return "true" if succesful. 
bool CProcessFeatures::MakeGrid() 
{ 
    if (m_filterDiam < 2) 
        return false;  // filter is too small 
    if (m_bounds.right < 2 * m_filterDiam || 
        m_bounds.bottom < 2 * m_filterDiam) 
        return false;   // filter is too big 
    int row = m_filterDiam - 1; // big filter just hits row 0 
    int col = m_filterDiam - 1; // big filter just hits column 0 
    double centerToCenter = m_filterDiam *  
        sqrt(PI_HALF / (ROOT_3 * m_overlap)); 
    int row_space = (int) centerToCenter;   
    if (row_space < 2) 
        row_space = 2; 
    int col_space = (int) (ROOT_3 *  centerToCenter); 
    if (col_space < 2) 
        col_space = 2; 
 
    m_locs.RemoveAll();  // get rid of any previous locations. 
    for (int format = 0; format < 2; format++) // overlapped grids 
    { 
        for (int c = col; c < m_bounds.right; c += col_space) 
        { 
            for (int r = row; r < m_bounds.bottom; r += row_space) 
            { 
                CLocation newLoc( c, r, m_filterDiam ); // create location 
                newLoc.SetSmall( false ); // look at bigger region 
                RECT big = newLoc.GetRect(); 
                if (big.left   >= m_bounds.left && 
                    big.top    >= m_bounds.top  && 
                    big.right  <= m_bounds.right && 
                    big.bottom <= m_bounds.bottom) 
                {   // add this location to the list. 
                    m_locs.Add( newLoc ); 
                } 
            } 
        }	// loop on columns 
        row += row_space/2; 
		col += col_space/2; 
    }	// loop on overlapped grids 
    return true; 
} 
/*------------------------------------------------------------------------*/ 
// Construct CLocations at which to filter a 1D test image. 
// Return "true" if succesful. 
bool CProcessFeatures::Make1DGrid() 
{ 
    if (m_filterDiam < 2) 
        return false;  // filter is too small 
    // assumes m_bounds.left = m_bounds.top = 0 
    if (m_bounds.right < 2 * m_filterDiam || 
        m_bounds.bottom < 2 * m_filterDiam) 
        return false;   // filter is too big 
    int row = m_bounds.bottom/2; // only the center row 
    int col = (m_bounds.right- m_filterDiam)/2 - 1; // small filter just hits center 
    int col_last = (m_bounds.right + m_filterDiam)/2; // small filter just hits center 
    int col_space = 1; 
 
    m_locs.RemoveAll();  // get rid of any previous locations. 
    for (int c = col; c < col_last; c += col_space) 
    { 
        CLocation newLoc( c, row, m_filterDiam ); // create location 
        newLoc.SetSmall( false ); // look at bigger region 
        RECT big = newLoc.GetRect(); 
        if (big.left   >= m_bounds.left && 
            big.top    >= m_bounds.top  && 
            big.right  <= m_bounds.right && 
            big.bottom <= m_bounds.bottom) 
        {   // add this location to the list. 
            m_locs.Add( newLoc ); 
        } 
    }	// loop on columns 
    return true; 
} 
/*------------------------------------------------------------------------*/ 
// Construct CLocations at which to filter a 2D test image. 
// Return "true" if succesful. 
bool CProcessFeatures::Make2DGrid() 
{ 
    if (m_filterDiam < 2) 
        return false;  // filter is too small 
    if (m_bounds.right < 2 * m_filterDiam || 
        m_bounds.bottom < 2 * m_filterDiam) 
        return false;   // filter is too big 
    int row = (m_bounds.bottom- m_filterDiam)/2 - 1; // small filter just hits center 
    int col = (m_bounds.right- m_filterDiam)/2 - 1; // small filter just hits center 
    int row_space = 1; 
    int col_space = 1; 
    float radius_squared = m_filterDiam * m_filterDiam * (float) 0.25; 
    float row_center = (m_bounds.bottom - 1 )* (float) 0.5; 
    float col_center = (m_bounds.right - 1 )* (float) 0.5; 
    int row_last = (int) (row_center + m_filterDiam * (float) 0.5); 
    int col_last = (int) (col_center + m_filterDiam * (float) 0.5); 
 
    m_locs.RemoveAll();  // get rid of any previous locations. 
    for (int c = col; c < col_last; c += col_space) 
    { 
        float c_sqr = (c - col_center) * (c - col_center); 
        for (int r = row; r < row_last; r += row_space) 
        {    
            if ((r - row_center) * (r - row_center) + c_sqr > radius_squared) 
                continue;  // filter would not hot the center point 
            CLocation newLoc( c, r, m_filterDiam ); // create location 
            newLoc.SetSmall( false ); // look at bigger region 
            RECT big = newLoc.GetRect(); 
            if (big.left   >= m_bounds.left && 
                big.top    >= m_bounds.top  && 
                big.right  <= m_bounds.right && 
                big.bottom <= m_bounds.bottom) 
            {   // add this location to the list. 
                m_locs.Add( newLoc ); 
            } 
        } 
    }	// loop on columns 
    return true; 
} 
/*------------------------------------------------------------------------*/ 
// Read a file and construct CLocations at which to filter 
// If m_locs is empty, fill it by reading from m_gridFileName. 
// Return "true" if succesful. 
bool CProcessFeatures::ReadGrid() 
{ 
	PROFILE("ReadGrid"); 
 
    ifstream inGrid( m_gridFileName, ios::in ); 
//  if (!inGrid) 
//		return false; // can't open file 
 
//  If the file read does not work, provide these hardcoded values. 
    int gridMethod = 0; 
    float threshold;  // ignored; m_threshold is hard-coded. 
    int diameter = 32; // 8; 
    int row = 35;  // 3 or 7; 
    int col = 32;  // 3 or 7; 
    int row_space = 27; // 8; 
    int col_space = 54; // 8; 
 
    int image_rows = 64; // only used for old style gridMethod 1. 
    int image_cols = 64; 
    int format = 0; 
    bool got_some = false; 
 
    m_locs.RemoveAll();  // get rid of any previous locations. 
 
//  if (inGrid >> gridMethod) 
//      return false;  // can't read 
 
    switch (gridMethod) 
    { 
    case 1: // old style; everything given explicitly 
        if (inGrid >>threshold) 
            return got_some; // can't read threshold. 
        if (inGrid >> diameter) 
            return got_some; // can't read. 
        if (diameter < 3) 
            return got_some; 
        if (inGrid >> format >> image_rows >> image_cols) // ignore 
            return got_some; // can't read. 
        while (inGrid >> row >> col) 
        { 
            CLocation newLoc( col, row, diameter ); // create location 
            newLoc.SetSmall( false ); // look at bigger region 
            RECT big = newLoc.GetRect(); 
            if (big.left   >= m_bounds.left && 
                big.top    >= m_bounds.top  && 
                big.right  <= m_bounds.right && 
                big.bottom <= m_bounds.bottom) 
            {   // add this location to the list. 
                m_locs.Add( newLoc ); 
            } 
        } 
        break; 
 
    case 0: // read failure; use hardcoded values. 
        for (format = 0; format < 2; format++) // overlapped grids 
        { 
            for (int c = col; c < m_bounds.right; c += col_space) 
            { 
                for (int r = row; r < m_bounds.bottom; r += row_space) 
                { 
                    CLocation newLoc( c, r, diameter ); // create location 
                    newLoc.SetSmall( false ); // look at bigger region 
                    RECT big = newLoc.GetRect(); 
                    if (big.left   >= m_bounds.left && 
                        big.top    >= m_bounds.top  && 
                        big.right  <= m_bounds.right && 
                        big.bottom <= m_bounds.bottom) 
                    {   // add this location to the list. 
                        m_locs.Add( newLoc ); 
                    } 
                } 
            }	// loop on columns 
            row += row_space/2; 
			col += col_space/2; 
        }	// loop on overlapped grids 
        got_some = true; 
        break; 
 
    case 2: // new style; use repeat block 
        if (inGrid >>threshold) 
            return got_some; // can't read threshold. 
        while (inGrid >> diameter >> row >> col >> row_space >> col_space) 
        { 
            if (diameter < 3) 
                return got_some; // diameter too small 
            for (int c = col; c < m_bounds.right; c += col_space) 
            { 
                for (int r = row; r < m_bounds.bottom; r += row_space) 
                { 
                    CLocation newLoc( c, r, diameter ); // create location 
                    newLoc.SetSmall( false ); // look at bigger region 
                    RECT big = newLoc.GetRect(); 
                    if (big.left   >= m_bounds.left && 
                        big.top    >= m_bounds.top  && 
                        big.right  <= m_bounds.right && 
                        big.bottom <= m_bounds.bottom) 
                    {   // add this location to the list. 
                        m_locs.Add( newLoc ); 
                    } 
                } 
            } 
        } // exit loop when no more to read. 
        break; 
 
    case 97: // read failure; use hardcoded values. 
        row = 5;  // larger grid with 12 pixel filters. 
        col = 5; 
        row_space = 12; 
        col_space = 12; 
        for (format = 0; format < 2; format++) 
        { 
            for (int c = col; c < m_bounds.right; c += col_space) 
            { 
                for (int r = row; r < m_bounds.bottom; r += row_space) 
                { 
                    CLocation newLoc( c, r, diameter ); // create location 
                    newLoc.SetSmall( false ); // look at bigger region 
                    RECT big = newLoc.GetRect(); 
                    if (big.left   >= m_bounds.left && 
                        big.top    >= m_bounds.top  && 
                        big.right  <= m_bounds.right && 
                        big.bottom <= m_bounds.bottom) 
                    {   // add this location to the list. 
                        m_locs.Add( newLoc ); 
                    } 
                } 
            } 
            row = col = 11; 
        } 
        got_some = true; 
        break; 
 
 
    case 98: // Special test values for debug on cube 
        { 
        CLocation newLoc( 16, 31, diameter ); // white to black 
        m_locs.Add( newLoc ); 
        newLoc.Initialize( 24, 16, diameter ); // 211 deg 
        m_locs.Add( newLoc ); 
        } 
        got_some = true; 
        break; 
 
    case 99: // Special test values for 64 x 64 square 
        { 
        CLocation newLoc( 8, 8, diameter ); // all gray 
        m_locs.Add( newLoc ); 
        newLoc.Initialize( 32, 32, diameter ); // all dark 
        m_locs.Add( newLoc ); 
        newLoc.Initialize( 32, 16, diameter );  
        m_locs.Add( newLoc ); // top light, lower dark 
        newLoc.Initialize( 32, 48, diameter ); 
        m_locs.Add( newLoc ); // top dark, lower light 
        newLoc.Initialize( 48, 32, diameter ); 
        m_locs.Add( newLoc ); // left dark, right light 
        newLoc.Initialize( 15, 32, diameter ); 
        m_locs.Add( newLoc ); // left light, right dark 
        newLoc.Initialize( 16, 32, diameter ); 
        m_locs.Add( newLoc ); // left light, right dark 
        newLoc.Initialize( 17, 32, diameter ); 
        m_locs.Add( newLoc ); // center at 16.5 
        newLoc.Initialize( 18, 32, diameter ); 
        m_locs.Add( newLoc ); // center at 17.5 
        newLoc.Initialize( 19, 32, diameter ); 
        m_locs.Add( newLoc ); // center at 18.5 
        } 
        got_some = true; 
        break; 
 
    default: 
        return got_some;  // bad style 
 
    } 
    return got_some; 
} 
/*------------------------------------------------------------------------*/ 
// Set the grid file name to the new value. 
// If this is the same as before, do nothing. 
// Otherwise, read and validate the new grid file. 
// If the new grid is not valid, retain the old one  
// and return "false".  Otherwise return "true". 
bool CProcessFeatures::SetGrid( CString gridFileName ) 
{ 
    if (gridFileName == m_gridFileName) 
        return true; 
    CString oldName = m_gridFileName; 
    if (ReadGrid()) 
        return true; 
    // else failed to read new file 
    m_gridFileName = oldName; 
    ReadGrid(); 
    return false; 
} 
/*------------------------------------------------------------------------*/ 
// Filter the image at a particular location 
// return true if successful. 
bool CProcessFeatures::Filter( CLocation& location ) 
{ 
	PROFILE("Filter"); 
 
    location.SetSmall(true); 
    location.SetOdd(false); 
    CFilter even, odd; 
    bool newFilter = false; 
	float magSqr, evenMagSqr, oddMagSqr; 
	float evenResults[4]; 
	float oddResults[3]; 
 
	int xcent, ycent, orn; 
    location.GetLocation(&xcent, &ycent); 
#ifdef TRACE_ME 
	ofstream traceFile( "Trace.txt", ios::app ); 
	traceFile << "Starting Filter at " << xcent << ", " << ycent << endl; 
	traceFile.close(); 
#endif // TRACE_ME 
    if (xcent == 62 && ycent == 146) 
//  if (xcent == 149 && ycent == 197) 
    {   // trouble spot on cube image. 
        int i = 0;  // do nothing statement for a breakpoint. 
    } 
 
	int bestBand = 1; 
    // Small diameter even filters 
    location.SetOdd(false); 
    location.SetSmall(true); 
    int diam = location.GetDiam(); 
    if (!m_filters.Lookup( diam, even )) 
    { 
        // create filters of this size 
        CFilter newFilters( diam, diam/m_subSampleDiam ); 
        newFilter = true; 
        // put filters into the CMap 
        m_filters.SetAt( diam, newFilters ); 
        if (!m_filters.Lookup( diam, even )) 
            return false;  // failed to put filter into CMap. 
    } 
#ifdef TRACE_ME 
	traceFile << "Starting 1st correlation at " << xcent << ", " << ycent << endl; 
	traceFile.close(); 
#endif // TRACE_ME 
    // Correlate small diameter even filters with subimage 
    evenMagSqr = correlate( location, even, bestBand); 
	// save the results 
	for (orn = 0; orn < even.m_orientations; orn ++) 
		evenResults[orn] = location.GetResult(orn); 
 
    // Small diameter odd filters 
    location.SetOdd(true); 
    diam = -diam;   // negative key for odd filters 
    if (!m_filters.Lookup( diam, odd )) 
    { 
        // create filters of this size 
        CFilter newFilters( diam,  diam/m_subSampleDiam ); 
        newFilter = true; 
        // put filters into the CMap 
        m_filters.SetAt( diam, newFilters ); 
        if (!m_filters.Lookup( diam, odd )) 
            return false;  // failed to put filter into CMap. 
    } 
    // Correlate small diameter odd filters with subimage 
    oddMagSqr = correlate( location, odd, bestBand); 
	// save the results 
	for (orn = 0; orn < odd.m_orientations; orn ++) 
		oddResults[orn] = location.GetResult(orn); 
 
    if (newFilter) 
    {   // compute cross correlations. 
        even.fl_correlate( odd ); 
        odd.fl_correlate( even ); 
        // put filters into the CMap 
        m_filters.SetAt( diam, odd ); // negative diam for odd 
        diam = -diam;   // positive key for even filters 
        m_filters.SetAt( diam, even ); 
		newFilter = false; 
    } 
 
	// Now see if there are other color bands that give more response 
	// Ideally, the image is preprocessed into a luminance and 2 chrominance bands. 
	// More often, we are looking at red, green and blue bands. 
	// The edge found will be the maximum over all bands for the small diameter. 
	// Possible undesired effect: it might find edges in alpha 
	magSqr = evenMagSqr + oddMagSqr; 
	for (int band = 2; band <= m_pImage->NBands(); band++) 
	{ 
		// Correlate small diameter even filters with subimage 
		location.SetOdd(false); 
		evenMagSqr = correlate( location, even, band); 
		// Correlate small diameter odd filters with subimage 
	    location.SetOdd(true); 
		oddMagSqr = correlate( location, odd, band); 
		if (evenMagSqr + oddMagSqr > magSqr) 
		{ 
			magSqr = evenMagSqr + oddMagSqr; 
			bestBand = band; 
			for (orn = 0; orn < odd.m_orientations; orn ++) 
				oddResults[orn] = location.GetResult(orn); 
			location.SetOdd(false); 
			for (orn = 0; orn < even.m_orientations; orn ++) 
				evenResults[orn] = location.GetResult(orn); 
		} 
	} 
	// Restore the results from the best band 
	for (orn = 0; orn < odd.m_orientations; orn ++) 
		location.SetResult( orn, oddResults[orn]); 
	location.SetOdd(false); 
	for (orn = 0; orn < even.m_orientations; orn ++) 
		location.SetResult( orn, evenResults[orn]); 
#ifdef TRACE_ME 
	traceFile << "Done with small correlations at " << xcent << ", " << ycent << endl; 
	traceFile.close(); 
#endif // TRACE_ME 
 
 
    // Next do the even large diameter 
    location.SetOdd(false); 
    location.SetSmall(false); 
    diam = location.GetDiam(); 
    if (!m_filters.Lookup( diam, even )) 
    { 
        // create filters of this size 
        CFilter newFilters( diam, diam/m_subSampleDiam ); 
        // put filters into the CMap 
        newFilter = true; 
        m_filters.SetAt( diam, newFilters ); 
        if (!m_filters.Lookup( diam, even )) 
            return false;  // failed to put filter into CMap. 
    } 
    // Correlate large diameter even filters with subimage 
    correlate( location, even, bestBand); 
 
    // Next is the odd large diameter 
    diam = -diam;   // negative key for odd filters 
    location.SetOdd(true); 
    if (!m_filters.Lookup( diam, odd )) 
    { 
        // create filters of this size 
        CFilter newFilters( diam, diam/m_subSampleDiam ); 
        newFilter = true; 
        // put filters into the CMap 
        m_filters.SetAt( diam, newFilters ); 
        if (!m_filters.Lookup( diam, odd )) 
            return false;  // failed to put filter into CMap. 
    } 
    // Correlate large diameter odd filters with subimage 
    correlate( location, odd, bestBand); 
 
    if (newFilter) 
    {   // compute cross correlations. 
        even.fl_correlate( odd ); 
        odd.fl_correlate( even ); 
        // put filters into the CMap 
        m_filters.SetAt( diam, odd ); // negative diam for odd 
        diam = -diam;   // positive key for even filters 
        m_filters.SetAt( diam, even ); 
    } 
#ifdef TRACE_ME 
	traceFile << "Finished Filter at " << xcent << ", " << ycent << endl; 
	traceFile.close(); 
#endif // TRACE_ME 
     
    return true;   
} 
/*------------------------------------------------------------------------*/ 
// Find the features at a particular location 
bool CProcessFeatures::Interpret( CLocation& loc ) 
{ 
	PROFILE("Interpret"); 
 
    float corrEven;  // not to be confused with the array of same name in loc 
    float corrOdd; 
    float magn_sqr; // squared magnitude of response at angles 
 
    if (loc.GetStrengthBound() < m_threshold * PRE_THRESH) 
        return true;  // Nothing is there 
 
    // This is for debugging behavior at a particular place. 
    int xcent, ycent; 
    loc.GetLocation(&xcent, &ycent); 
    if (xcent == 62 && ycent == 146) 
//  if (xcent == 149 && ycent == 197) 
    {   // trouble spot on cube image. 
        magn_sqr = 0;  // do nothing statement for a breakpoint. 
#ifdef TRACE_ME 
	ofstream traceFile( "Trace.txt", ios::app ); 
	traceFile << "Starting Interpret" << endl; 
	traceFile.close(); 
#endif // TRACE_ME 
    } 
 
    // Find the feature orientation 
    float angle = loc.GetAngle( &magn_sqr,  
     &corrEven, &corrOdd, m_fixedDirection); 
 
//  Retreive the filter used 
    // Small diameter odd filters 
    loc.SetSmall(true); 
    loc.SetOdd(true); 
    int diam = loc.GetDiam(); 
    diam = -diam;   // negative key for odd filters 
    CFilter filter; 
    if (!m_filters.Lookup( diam, filter )) 
        return false;  // can't find filter 
    float base_strength = filter.m_maxResponse / 2; 
//  float base_strength = STRENGTH_SCALE * filter.m_self_ovlap[0]; 
 
    loc.SetSmall(false); 
    int big_diam = loc.GetDiam(); 
    big_diam = -big_diam;   // negative key for odd filters 
    if (!m_filters.Lookup( big_diam, filter )) 
        return false;  // can't find filter 
    float big_base_strength = filter.m_maxResponse / 2; 
//  float big_base_strength = STRENGTH_SCALE * filter.m_self_ovlap[0]; 
    float strength = (float) sqrt((double)(magn_sqr)) / base_strength; 
 
    if (m_threshold < (float) 0.0) 
        m_threshold *= -strength; 
    if (strength > m_threshold * PRE_THRESH) 
    {   /* Strength is significant within the residual image */ 
 
 
        int x, y, diam; 
        diam = loc.GetLocation( &x, &y); 
 
        //  Retreive the filter used 
        // Small diameter even filters 
        loc.SetSmall(true); 
        loc.SetOdd(false); 
        diam = loc.GetDiam(); // positive key for even filters 
        if (!m_filters.Lookup( diam, filter )) 
            return false;  // can't find filter 
        float center_lobe = filter.m_center_lobe; 
 
        // Large diameter even filters 
        loc.SetSmall(false); 
        big_diam = loc.GetDiam(); // positive key for even filters 
        if (!m_filters.Lookup( big_diam, filter )) 
            return false;  // can't find filter 
        float big_center_lobe = filter.m_center_lobe; 
 
        float width_tuning = (float) 0.5 * (center_lobe + 
         big_center_lobe); 
 
        // Steer big filters to this angle 
        loc.SetSmall(false); 
        loc.SetOdd( false); 
        float corrBigEven = loc.Steer( angle); 
        loc.SetOdd( true); 
        float corrBigOdd  = loc.Steer( angle); 
 
        // Create a CFeature 
        CFeature feature( x, y, diam, angle, 
         corrEven, corrOdd, corrBigEven, corrBigOdd, 
         base_strength, big_base_strength, big_diam,  
         big_center_lobe, width_tuning); 
 
        feature.multi_find(0, m_viewEdgesOnly); 
        // steer to other orientations. 
//      for (int i = 1; i < 4; i++) 
		int i = 2; 
        { 
            angle += PI_1/2;  // PI_1/4 to use m_45 
            // Steer small filters to this angle 
            loc.SetSmall(true); 
            loc.SetOdd( false); 
            corrEven = loc.Steer( angle); 
            loc.SetOdd( true); 
            corrOdd  = loc.Steer( angle); 
            // Steer big filters to this angle 
            loc.SetSmall(false); 
            loc.SetOdd( false); 
            corrBigEven = loc.Steer( angle); 
            loc.SetOdd( true); 
            corrBigOdd  = loc.Steer( angle); 
            feature.AddSecondary( i, angle, corrEven, corrOdd, corrBigEven,  
             corrBigOdd); 
            feature.multi_find(i, m_viewEdgesOnly); 
        } 
		feature.corners(); 
 
        // Add the feature to m_features if it is above threshold. 
		// We do not have any results from lateral antagonism at this point. 
		// GetAnagonism is used here for future expansion to the case where 
		// we are processing multiple frames and have previous results available. 
        if (feature.GetStrength() /*  + m_lateralFactor * feature.GetAntagonism() */ 
			>= m_threshold * PRE_THRESH) 
        { 
            m_features.Add( feature ); 
        } 
    } 
    return true;   
} 
/*------------------------------------------------------------------------*/ 
/* This is used for correlating an image with a kernel. */ 
 
// Correlate the image at the given location with the filters. 
// returned value is the squared sum of magnitudes at all orientations. 
float CProcessFeatures::correlate(  
        CLocation& location,  
        CFilter& filters, 
		int band) 
{ 
	PROFILE("correlate"); 
 
    CDefaultImage region = m_pImage->SubImage(location.GetRect()); 
    int orn, r, c; 
    float result; 
 
    int first_col = region.Left(); 
    int first_row = region.Top(); 
    int last_col = region.Right(); 
    int last_row = region.Bottom(); 
 
    float *pKern = filters.kern; 
	int nbands = region.NBands(); 
 
	float magSqr = 0; 
#ifdef DO_COLORS 
// Add support for color images   5/31/03 
	for (orn = 0; orn < filters.m_orientations; orn ++) 
	{ 
		result = (float) 0.0; 
		for (r = first_row; r < last_row; r+=filters.m_sampleSpacing) 
		{ 
			for (c = first_col; c < last_col; c+=filters.m_sampleSpacing) 
			{ 
				register unsigned char pel = *region.PbPixel(c,r,band); 
				result += (float) pel * *pKern++; 
			} 
		} 
		location.SetResult( orn, result ); 
		magSqr += result * result; 
	} 
#else  // Gray level byte image 
	/*  Could check for gray scale with: 
	type_info pixel_type = typeid( region.Pixel(first_col, first_row)); 
    type_info byte = typeid( char); 
    if (region.NBands() == 1 && pixel_type == byte) 
	*/ 
//#define DEBUG_ME 1 
#ifdef DEBUG_ME 
	float kernels[60]; 
	float results[60]; 
	int pixels[60]; 
    int xcent, ycent; 
    location.GetLocation(&xcent, &ycent); 
    if (xcent == 77 && ycent == 70) 
	{ 
		ofstream outFile( "Convolve.txt", ios::out ); 
		// set to show 4 decimal places 
		int oldPrecision = outFile.precision(4); 
		outFile.setf( ios::fixed, ios::floatfield ); 
		 
 
		for (orn = 0; orn < filters.m_orientations; orn ++) 
		{ 
			result = (float) 0.0; 
			for (r = first_row; r < last_row; r+=filters.m_sampleSpacing) 
			{ 
				outFile << endl << "Orientation:\t" << orn << "\tRow:\t" << r << endl; 
				unsigned char *pPixel = region.RowPointer(r); 
				for (c = first_col; c < last_col; c+=filters.m_sampleSpacing) 
				{ 
					kernels[c-first_col] = *pKern; 
					pixels[c-first_col] = pPixel[c]; 
					result += pPixel[c] * *pKern++; 
					results[c-first_col] = result; 
				} 
				for (c = first_col; c < last_col; c+=filters.m_sampleSpacing) 
					outFile << pixels[c-first_col] << '\t'; 
				outFile << endl; 
				for (c = first_col; c < last_col; c+=filters.m_sampleSpacing) 
					outFile << kernels[c-first_col] << '\t'; 
				outFile << endl; 
				for (c = first_col; c < last_col; c+=filters.m_sampleSpacing) 
					outFile << results[c-first_col] << '\t'; 
				outFile << endl; 
			} 
		} 
		outFile.precision(oldPrecision); // restore previous precision 
		outFile.close(); 
		pKern = filters.kern; 
    } 
#endif  // DEBUG_ME 
	// This should be a bit faster than the more general method for color images. 
	if (filters.m_sampleSpacing == 1) 
	{ 
		for (orn = 0; orn < filters.m_orientations; orn ++) 
		{ 
			result = (float) 0.0; 
			for (r = first_row; r < last_row; r++) 
			{ 
				unsigned char *pPixel = region.RowPointer(r); 
				for (c = first_col; c < last_col; c++) 
				{ 
					result += pPixel[c] * *pKern++; 
				} 
			} 
			location.SetResult( orn, result ); 
			magSqr += result * result; 
		} 
	} 
	else 
	{ 
		for (orn = 0; orn < filters.m_orientations; orn ++) 
		{ 
			result = (float) 0.0; 
			for (r = first_row; r < last_row; r+=filters.m_sampleSpacing) 
			{ 
				unsigned char *pPixel = region.RowPointer(r); 
				for (c = first_col; c < last_col; c+=filters.m_sampleSpacing) 
				{ 
					result += pPixel[c] * *pKern++; 
				} 
			} 
			location.SetResult( orn, result ); 
			magSqr += result * result; 
		} 
	} 
#endif // DO_COLORS 
	return magSqr; 
} 
/*------------------------------------------------------------------------*/ 
// Fill in CFeature::m_pNbr with pointers to adjacent features 
void CProcessFeatures::FindNeighbors() 
{ 
	PROFILE("FindNeighbors"); 
 
    for ( int n = 0; n < m_features.GetSize(); n++)  
    { 
		for ( int m = n+1; m < m_features.GetSize(); m++)  
		{ 
			// Neighbor if the distance between centers is less than the sum of radii. 
			m_features[n].IsNeighbor(&m_features[m]); 
		} 
    } 
 
} 
/*------------------------------------------------------------------------*/ 
// Write text of features.  Return "true" if succesful. 
bool CProcessFeatures::WriteFeatureList 
(   int verbosity,  // 0 for terse, 1 or 2 for more verbose  
    int mode,       // ios::app to append to file; ios::out for new file 
    char* title,    // put this on first line 
    float noise     // noise level of test image 
) 
{ 
	PROFILE("WriteFeatureList"); 
    ofstream outFile( "Results.txt", mode ); 
    if (!outFile) 
    { 
        return false; // can't open file 
    } 
 
    int i; 
    if (mode == ios::out) 
    {  // put header on new file 
        outFile << endl << title << endl; 
        if (m_answer.m_type != eNOT_AVAILABE) 
        { 
            outFile << "EAngle"  << '\t'; 
            outFile << "EIntens" << '\t'; 
            outFile << "ErrPos" << '\t'; 
            outFile << "ErrType"   << '\t'; 
            outFile << "ErrWide"  << '\t'; 
            outFile << "Noise"  << '\t'; 
        } 
        outFile << "Region center" << '\t' << "Diam" << '\t' << "Lateral" << '\t' 
			<< "X" << '\t' << "Y" << '\t' << "Type2D" << '\t'; 
        int numSteers = verbosity == 0? 1: 2; 
        for (i = 0; i < numSteers; i++) 
        { 
            outFile << "Angle"  << '\t'; 
            outFile << "Intens" << '\t'; 
            outFile << "Column" << '\t'; 
            outFile << "Row"    << '\t'; 
            outFile << "Type"   << '\t'; 
            outFile << "Pos"    << '\t'; 
            outFile << "Width"  << '\t'; 
            if (verbosity > 1) 
            { 
                outFile << "Odd"      << '\t'; 
                outFile << "Even"     << '\t'; 
                outFile << "BigOdd"  << '\t'; 
                outFile << "BigEven" << '\t'; 
            } 
        } 
        outFile << endl; 
    } // end of header 
 
    // set to show 2 decimal places 
    int oldPrecision = outFile.precision(2); 
    outFile.setf( ios::fixed, ios::floatfield ); 
 
    if (m_answer.m_type != eNOT_AVAILABE) 
    { 
        /*  Write the answer * 
        outFile << "\t\t\t\t\t";  // skip errors 
        outFile << "\t\t\t";  // skip region center and diam 
        m_answer.write( outFile, true, 1 ); 
        outFile << endl; 
        */ 
        for (i = 0; i < m_features.GetSize(); i++) 
        { 
			float angleError = -m_answer.m_degrees + m_features[i].GetDegrees(); 
			if (angleError > 180) 
				angleError -= 360; 
			if (m_answer.m_type == eWHITE_LINE || 
				m_answer.m_type == eBLACK_LINE) 
			{ 
				while(angleError > 90) 
					angleError -= 180; 
				while(angleError < -90) 
					angleError += 180; 
			} 
            outFile << angleError << '\t'; 
			// 6/11/03  Strength output changed to pixels. 
			// Previously, it had been a percent. 
            outFile << (MAX_PIXEL-MIN_PIXEL)/2 * (-m_answer.m_strength + m_features[i].GetStrength()) << '\t'; 
            float distance =  m_features[i].PositionDifference( m_answer ); 
            outFile << distance << '\t'; 
            int typeError = m_features[i].TypeDifference( m_answer ); 
            outFile << typeError << '\t'; 
            outFile << -m_answer.m_width + m_features[i].GetWidth() << '\t'; 
            outFile << noise << '\t'; 
            m_features[i].write( outFile, true, verbosity ); 
        } 
    } 
    else 
    { 
        for (i = 0; i < m_features.GetSize(); i++) 
        { 
            m_features[i].write( outFile, true, verbosity ); 
        } 
    } 
//  outFile << endl; 
    outFile.precision(oldPrecision); // restore previous precision 
    outFile.close(); 
    return true; 
} 
/*------------------------------------------------------------------------*/ 
// Overlay graphics on the image 
void CProcessFeatures::MakeFeatureImage(CVisRGBAByteImage& featureImage) 
{ 
    // for debugging, show the biggest kernels. 
//  plotKernels( featureImage );  
    drawLines( featureImage );  
} 
/*------------------------------------------------------------------------*/ 
void CProcessFeatures::plotKernels(CVisRGBAByteImage& featureImage) 
{ 
    // for debugging, show the biggest kernels. 
    // find the biggest kernel 
    int max_diam = 0; 
    int diam; 
    CFilter filter; 
    CFilter biggest; 
    POSITION pos = m_filters.GetStartPosition( ); 
    while ( pos != NULL ) 
    { 
        m_filters.GetNextAssoc( pos, diam, filter ); 
        if (diam > max_diam) 
        { 
            max_diam = diam; 
            biggest = filter; 
        } 
    } 
    // find largest value in the filter 
    int size = (max_diam / biggest.m_sampleSpacing) * (max_diam / biggest.m_sampleSpacing); 
    float max = 0; 
    for (int i = 0; i < size; i++) 
    { 
        if (FL_ABS(biggest.kern[i]) > max ) 
            max = FL_ABS(biggest.kern[i]); 
    } 
    float scale = (float) .92 * 128 / max; // Assume 256 gray levels. 
 
    i = 0; 
    int c, r; 
    for (int orn = 0; orn < biggest.m_orientations; orn++) 
    { 
        for (c = 0; c < max_diam; c += biggest.m_sampleSpacing) 
        { 
            for (r = orn*max_diam; r < orn*max_diam+max_diam; r += biggest.m_sampleSpacing) 
            { 
                int value = 128 + biggest.kern[i++] * scale; 
                featureImage.Pixel(r,c) = value; 
            } 
        } 
    } 
 
    // Now show the quadrature phase filter. 
    diam = -max_diam; 
    m_filters.Lookup( diam, biggest ); 
    i = 0; 
    for (orn = 0; orn < biggest.m_orientations; orn++) 
    { 
        for (c = max_diam; c < 2*max_diam; c += biggest.m_sampleSpacing) 
        { 
            for (r = orn*max_diam; r < orn*max_diam+max_diam; r += biggest.m_sampleSpacing) 
            { 
                int value = 128 + biggest.kern[i++] * scale; 
                featureImage.Pixel(r,c) = value; 
            } 
        } 
    } 
 
} 
/*------------------------------------------------------------------------*/ 
// Overlay graphics on the image 
void CProcessFeatures::drawLines(CVisRGBAByteImage& image) 
{ 
	PROFILE("drawLines"); 
 
    HDC hdcImage = image.Hdc(); 
    // create a handle for a device context overlaying the image 
    if (hdcImage != 0) 
    { 
        int oldThickness = 2; 
        COLORREF oldColor = EDGE_COLOR; 
        HPEN hPen = CreatePen( PS_SOLID, 1, oldColor); 
        HPEN oldPen; 
        oldPen = (HPEN) SelectObject( hdcImage, hPen ); 
 
        Outline shape; 
		bool useMain = true; 
        for (int i = 0; i < m_features.GetSize(); ) 
        { 
            shape.segment.RemoveAll(); 
			// Find how to draw the feature 
			m_features[i].CadData( shape, useMain ); 
//			m_features[i].FindPlotInfo( shape, useMain ); 
 
            if (shape.segment.GetSize() > 1) 
            { 
                if (shape.color != oldColor || shape.thickness != oldThickness) 
                {   // select a pen based on strength and type of feature 
                    switch (3 /* shape.thickness */) 
                    { 
                    case 1: 
                        hPen = CreatePen( PS_DOT, 1, shape.color); 
                        break; 
                    case 2: 
                    default: 
                        hPen = CreatePen( PS_DASH, 1, shape.color); 
                        break; 
                    case 3: 
                        hPen = CreatePen( PS_SOLID, 1, shape.color); 
                        break; 
                    case 4: 
                        hPen = CreatePen( PS_SOLID, 2, shape.color); 
                        break; 
                    } 
                    SelectObject( hdcImage, hPen ); 
                    oldColor = shape.color; 
                    oldThickness = shape.thickness; 
                }  // end of selecting new pen 
                POINT vertex = shape.segment[0]; 
				if (!m_viewEdgesOnly || shape.color == EDGE_COLOR) 
				{	// If requested, only draw red lines (step edges) 
					// I.e. suppress drawing bars (blue or yellow) 
					MoveToEx( hdcImage, vertex.x, vertex.y, NULL ); 
					for (int j = 1; j < shape.segment.GetSize(); j++) 
					{   // draw it 
						vertex = shape.segment[j]; 
						LineTo( hdcImage, vertex.x, vertex.y); 
					} 
				}  
            } // end of checking that there is something to draw 
			if (m_features[i].Is2D() && useMain == true) 
			{   // do an extra iteration for perpendicular direction. 
				useMain = false; 
			} 
			else 
			{ 
				if (m_viewGrid) 
				{ 
					hPen = CreatePen( PS_SOLID, 1, GRID_COLOR); 
					oldColor = GRID_COLOR; 
					SelectObject( hdcImage, hPen ); 
					// Draw a circle for the grid 
					m_features[i].DrawCircle( hdcImage ); 
				} 
				useMain = true; 
				i++; 
			} 
        }  // end of loop on features 
 
        // restore original pen 
        SelectObject( hdcImage, oldPen ); 
        image.DestroyHdc(); 
    } // end of check if we created a handle for a device context 
}