processfeatures.cpp

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

{
	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
}