testimage.cpp

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

// TestImage.cpp
/*
** Copyright (C) 2000 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 "TestImage.h"
#include "quad_dis.h"

/*------------------------------------------------------------------------*/
void CTestImage::MakeEdge
// Fill in the test image to make an edge with the given contrast.
( int contrast  /* in; gray-scale difference across the edge. */
)
{
    /* dark and light are the two gray levels to use*/
    int ndark = MEDIUM_GRAY - contrast/2;
    PIXEL dark = (PIXEL) (ndark < MIN_PIXEL?
        MIN_PIXEL: ndark);
    int nlight = ndark + contrast;
    PIXEL light = (PIXEL) (nlight > MAX_PIXEL?
        MAX_PIXEL: nlight);

    CRect full = Rect();  // The full test image
    int center = Height() / 2 - 1;
    CRect half = full;
    half.bottom = center; // Upper half of image
    SetRect(half);        // Activate only the upper half.
    FillPixels( light);   // upper half is light
    half = full;        
    half.top = center;    // Lower half of image
    SetRect(half);        // Activate lower half
    FillPixels( dark);    // bottom half is dark
    SetRect(full);        // whole image is active
}
/*------------------------------------------------------------------------*/
void CTestImage::MakeBar( int contrast, int barWidth)
{
    int ndark = MEDIUM_GRAY - contrast/2;
    PIXEL dark = (PIXEL) (ndark < MIN_PIXEL?
        MIN_PIXEL: ndark);
    int nlight = ndark + contrast;
    PIXEL light = (PIXEL) (nlight > MAX_PIXEL?
        MAX_PIXEL: nlight);

    CRect full = Rect();
    FillPixels( dark);  // whole image is dark
    int center = Height() / 2;
    CRect bar = full;
    bar.top = center - barWidth/2;
    bar.bottom = center + (barWidth+1)/2;
    SetRect(bar);
    FillPixels( light);   // add a light bar
    SetRect(full);       // whole image is active
}
/*------------------------------------------------------------------------*/
void CTestImage::MakeCorner
// Fill in the test image to make a corner with the given contrast.
( int contrast  /* in; gray-scale difference across the corner. */
)
{
    /* dark and light are the two gray levels to use*/
    int ndark = MEDIUM_GRAY - contrast/2;
    PIXEL dark = (PIXEL) (ndark < MIN_PIXEL?
        MIN_PIXEL: ndark);
    int nlight = ndark + contrast;
    PIXEL light = (PIXEL) (nlight > MAX_PIXEL?
        MAX_PIXEL: nlight);

    CRect full = Rect();  // The full test image
    int centerUp = Height() / 2 + 1;
    int centerHoriz = Width() / 2 + 1;
    CRect quarter = full;
    quarter.bottom = centerUp; // Upper half of image
    SetRect(quarter);        // Activate only the upper half.
    FillPixels( light);   // upper half is light

    quarter = full;        
    quarter.top = centerUp;    // Lower half of image
	quarter.right = centerHoriz;  // lower left
    SetRect(quarter);        // Activate lower left quarter
    FillPixels( light);    //  Make it light

    quarter = full;        
    quarter.top = centerUp;    // Lower half of image
	quarter.left = centerHoriz;  // lower right
    SetRect(quarter);        // Activate lower right quarter
    FillPixels( dark);    // make it dark
    SetRect(full);        // whole image is active
}
/*------------------------------------------------------------------------*/
void CTestImage::Rotate(int degrees)
{
    CRect bounds = Rect();
    // Retain the old image
    const CTestImage* pNewImage = this;
    CTestImage NewImage = *pNewImage;
    // copy the current image to a new image
    // Copy will make the current image be NewImage
    Copy(NewImage);
    // Do the rotation
    Warp( NewImage, (double) (degrees - 180));
    // Copy new image to current
    Copy(NewImage);
    // delete new image
//  delete pNewImage;
}
/*------------------------------------------------------------------------*/
/* Add Gaussian noise to an image */
void CTestImage::AddNoise
( float sigma )  /* in; Standard deviation for the noise. */
{
   /* Seed the random-number generator with current time so that
    * the numbers will be different every time we run.    */
    srand( (unsigned)time( NULL ) );
    CRect bounds = Rect();
    for (int row = bounds.top; row < bounds.bottom; row++)
    {
        for (int col = bounds.left; col < bounds.right; col++)
        {
            float x = (float) Pixel( col, row);
            float noisy = normal_rv (x, sigma); /* add gaussian noise */
            if (noisy < MIN_PIXEL)
                noisy = MIN_PIXEL;
            if (noisy > MAX_PIXEL)
                noisy = MAX_PIXEL;
            Pixel( col, row) = (PIXEL) (noisy + 0.5);
        }
   }
}
/*------------------------------------------------------------------------*/
/* Returns a normally distributed random variable. */
/* Tyler Folsom  10/28/92 */
float CTestImage::normal_rv 
(   float mean,  /* in; Mean value for the random variable */
    float sigma  /* in: Standard deviation */
)
{
    double r_num;
    float rand_noise;
    static struct density_function
    { double x; double area;}  ndf[] =  /* normal distribution function */
    {0.0, 0.5000,  0.1, 0.5398,  0.2, 0.5793,  0.3, 0.6179,  0.4, 0.6554,
     0.5, 0.6915,  0.6, 0.7257,  0.7, 0.7580,  0.8, 0.7881,  0.9, 0.8159,
     1.0, 0.8413,  1.1, 0.8643,  1.2, 0.8849,  1.3, 0.9032,  1.4, 0.9192,
     1.6, 0.9452,  1.8, 0.9641,  2.0, 0.9772,  2.2, 0.9861,  2.4, 0.9918,
     2.7, 0.9965,  3.0, 0.9987,  3.3, 0.9995,  3.6, 0.9998,  4.0, 1.0000};
    int i;

    r_num = rand() / double (RAND_MAX + 1);
    if (r_num < 0.5)
    {   /* ndf is symmetric about 0.5; reflect it. */
        sigma = -sigma;
        r_num = r_num + 0.5;
    }
    /* make it gaussian */
    for (i = 1; r_num > ndf[i].area; i++)
        ;
    /* Interpolate using ndf[i].x to get rand_noise. */
    rand_noise = (float) (ndf[i-1].x + (r_num - ndf[i-1].area) * 
     (ndf[i].x - ndf[i-1].x) / (ndf[i].area - ndf[i-1].area));

    return (mean + sigma * rand_noise);
}
/*------------------------------------------------------------------------*/
void CTestImage::Warp
(   CTestImage& newImage,
    double rot, // Rotation (degrees)
    double h1, double h2, // Translation (pixels from image center)
    double scale,  // scale factor: new height / old height
    double tilt_ax, // Tilt axis (degrees)
    double tilt    // Tilt angle (degrees)
)
{
    /* multiply coordinates of image by this matrix to get newImage coords */
    double old_new[4];
    /* multiply coordinates of newImage by this matrix to get original coords */
    double new_old[4];

    angles_matrix( scale, rot, tilt_ax, tilt, old_new);
    invert( old_new, new_old);
    warp_it( newImage, new_old, h1, h2);

}
/*----------------------------------------------------------------------*/
/* Convert the warp scale and angles to a linear transform matrix */
void CTestImage::angles_matrix( double scale, double rotat, double tilt_axis, 
    double tilt_angle, double *matrix)
{
    double scale_rot[4], tilts[4];
    double rot, tilt_ax, tilt;
    double cos_axis, sin_axis, cos_beta;

    /* Convert to radians */
    rot = rotat * PI / 180.0;
    tilt_ax = tilt_axis * PI / 180.0;
    tilt = tilt_angle * PI / 180.0;
    /* set up scaling and rotation */
    scale_rot[0] = scale * cos(rot);
    scale_rot[3] = scale_rot[0];
    scale_rot[2] = scale * sin(rot);
    scale_rot[1] = -scale_rot[2];
    /* set up tilt matrix */
    cos_axis = cos(tilt_ax);
    sin_axis = sin(tilt_ax);
    cos_beta = cos(tilt);
    tilts[0] = cos_beta + (1.0 - cos_beta) * cos_axis * cos_axis;
    tilts[1] =
    tilts[2] = (1.0 - cos_beta) * cos_axis * sin_axis;
    tilts[3] = cos_beta + (1.0 - cos_beta) * sin_axis * sin_axis;
    /* multiply these matrices to get transformation */
    matrix[0] = scale_rot[0] * tilts[0] + scale_rot[1] * tilts[2];
    matrix[1] = scale_rot[0] * tilts[1] + scale_rot[1] * tilts[3];
    matrix[2] = scale_rot[2] * tilts[0] + scale_rot[3] * tilts[2];
    matrix[3] = scale_rot[2] * tilts[1] + scale_rot[3] * tilts[3];
}
/*----------------------------------------------------------------------*/
/* Invert a 2x2 matrix */
void CTestImage::invert( double *matrix, double *inverse)
{
    double determinant;

    determinant = matrix[0] * matrix[3] - matrix[1] * matrix[2];
    if (fabs(determinant) <= 1.0e-20)
        determinant = 1.0;
    inverse[0] =  matrix[3] / determinant;
    inverse[1] = -matrix[1] / determinant;
    inverse[2] = -matrix[2] / determinant;
    inverse[3] =  matrix[0] / determinant;
}
/*----------------------------------------------------------------------*/
/* Transform iold to im_new by mapping coordinates:
  (x, y) = (matrix) (x - h1, y - h2)   */
void CTestImage::warp_it( CTestImage& newImage, double matrix[4], double h1, 
 double h2)
{
    CRect bounds = Rect();
    double x_old, y_old;
    int row_old, col_old;
    double x_new, y_new, y1, x1, new_value;

    double center_x = (bounds.right - bounds.left) / 2.0;
    double center_y = (bounds.bottom - bounds.top) / 2.0;
    CRect newSize = newImage.Rect();
    for (int row = newSize.top; row < newSize.bottom; row++)
    {
        y_new = (double) (row - h2);
        for (int col = newSize.left; col < newSize.right; col++)
        {
            x_new = (double) (col - h1);
            x_old = matrix[0] * (x_new - center_x) + matrix[1] * (y_new - 
             center_y) + center_x;
            y_old = matrix[2] * (x_new - center_x) + matrix[3] * (y_new -
             center_y) + center_y;
            row_old = (int) y_old;
            col_old = (int) x_old;

            if (col_old < bounds.left)
            {
                if (row_old < bounds.top)
                    newImage.Pixel(col, row) = Pixel(bounds.left, bounds.top);
                else if (row_old+1 >= bounds.bottom)
                    newImage.Pixel(col, row) = Pixel(bounds.left, bounds.bottom - 1);
                else
                {
                    y1 = y_old - (double) row_old;
                    new_value = y1 * Pixel(bounds.left, row_old+1) +
                     (1.0 - y1) * Pixel(bounds.left, row_old);
                    newImage.Pixel(col, row) = (PIXEL) (new_value + 0.5);
                }
            }
            else if (col_old >= bounds.right - 1)
            {
                if (row_old < bounds.top)
                    newImage.Pixel(col, row) = Pixel(bounds.right-1, bounds.top);
                else if (row_old >= bounds.bottom - 1)
                    newImage.Pixel(col, row) = Pixel(bounds.right-1, bounds.bottom - 1);
                else
                {
                    y1 = y_old - (double) row_old;
                    new_value = y1 * Pixel(bounds.right-1, row_old+1) +
                     (1.0 - y1) * Pixel(bounds.right-1, row_old);
                    newImage.Pixel(col, row) = (PIXEL) (new_value + 0.5);
                }
            }
            else   /* col_old is within bounds. */
            {
                x1 = x_old - (double) col_old;
                if (row_old < bounds.top)
                {
                    new_value = x1 * Pixel( col_old + 1, bounds.top)
                     + (1.0 - x1) * Pixel( col_old, bounds.top);
                }
                else if (row_old >= bounds.bottom - 1)
                {
                    new_value = x1 * Pixel( col_old + 1, bounds.bottom - 1)
                     + (1.0 - x1) * Pixel( col_old, bounds.bottom - 1);
                }
                else   /* Everything is within bounds */
                {
                    y1 = y_old - (double) row_old;
                    new_value = y1 * (x1 * Pixel( col_old + 1, row_old + 1)
                     + (1.0 - x1) * Pixel( col_old, row_old + 1))
                     + (1.0 - y1) * (x1 * Pixel( col_old + 1, row_old)
                     + (1.0 - x1) * Pixel( col_old, row_old));
                }   /* end of row_old if */
                newImage.Pixel(col, row) = (PIXEL) (new_value + 0.5);
            }   /* end of col_old if */
        }   /* end of loop on col */
    }   /* end of loop on row */
}