image_funcs.h

Go to the documentation of this file.

//
// Generic image functions
//
// Bala Amavasai (bala@amavasai.org)
// Mon Oct  1 15:41:44 BST 2001
//
// $Header: /cvs/mimas2/include/image_funcs.h,v 1.1.1.1 2005/08/09 15:37:45 engmb Exp $
//

#ifndef IMAGE_FUNCS_H
#define IMAGE_FUNCS_H

#include "mimasconfig.h"
#include <cassert>
#include <float.h>
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstdlib>
#include <values.h>
#include <cmath>
#include "image.h"
#include "image_conv.h"
#include "image_op.h"
#include "rgba.h"
#include "prop.h"
#include "property_image.h"
#include "image_draw.h"

namespace mimas {

template< typename T, typename TPtr >
T min_val( const const_image_ref< T, TPtr > &image )
{
  T retVal;
  if ( image.initialised() ) {
    int size = image.getWidth() * image.getHeight();
    retVal = *std::min_element( image.rawData(),
                                image.rawData() + size );
  } else
    retVal = T();
  return retVal;
}

template< typename T, typename TPtr >
T max_val( const const_image_ref< T, TPtr > &image )
{
  T retVal;
  if ( image.initialised() ) {
    int size = image.getWidth() * image.getHeight();
    retVal = *std::max_element( image.rawData(),
                                image.rawData() + size );
  } else
    retVal = T();
  return retVal;
}

template< typename T >
image_ref< T > &normaliseIt( image_ref< T > &im,
                             const T &val1, const T &val2 )
{
  _normalise< T > f( min_val( im ), max_val( im ), val1, val2 );
  return image_apply( im, im,
                      _multi_help1< T, T, _normalise< T > >( f ) );
}

template< typename T, typename TPtr >
image< T > normalise( const const_image_ref< T, TPtr > &im,
                      const T &val1, const T &val2 )
{
  if ( im.initialised() ) {// Image must not be empty.
    return image_func< T >( im,
                            _normalise< T >( min_val( im ),
                                             max_val( im ), val1, val2 ) );
  } else
    return image< T >();
}

template<typename T>
void equaliseIt( image_ref<T> &imagein )
{

    int H[256], TR[256];

    for (int count=0; count<256; count++)
  H[count]=0;

    int intensity;
    for (int j=0; j<imagein.getHeight(); j++)
        for (int i=0; i<imagein.getWidth(); i++)
  {
      intensity=imagein.pixel(i,j);
       assert( intensity <= 255 );
      H[intensity]+=1;
  }
    
   // convert to cumulative histogram 
    int imagesize=imagein.getWidth()*imagein.getHeight();
    for (int count=1; count<256; count++)
    {
  H[count]=H[count]+H[count-1];
  TR[count-1]=(255*H[count-1])/imagesize;
    }
    TR[255]=(255*H[255])/imagesize;
     
    for (int j=0; j<imagein.getHeight(); j++)
        for (int i=0; i<imagein.getWidth(); i++)
          imagein.pixel(i,j)=(T)TR[imagein.pixel(i,j)];
}

template< typename T >
image_ref< T > &bilevelIt( image_ref< T > &im, T threshval,
                           T val1, T val2 )
{
  std::binder2nd< _bilevel< T > > f( _bilevel< T >( val1, val2 ),
                                     threshval );
  return image_apply
    ( im, im,
      _multi_help1< T, T, std::binder2nd< _bilevel< T > > >( f ) );
}

template< typename T, typename TPtr >
image< T > bilevel( const const_image_ref< T, TPtr > &im,
                    T threshval, T val1, T val2 )
{
  return image_func< T >( im,
                          std::bind2nd( _bilevel< T >( val1, val2 ),
                                        threshval ) );
}

template< typename T >
image_ref< T > &thresholdIt( image_ref< T > &im, T threshval )
{
  std::binder2nd< _threshold< T > > f( _threshold< T >(), threshval );
  return image_apply
    ( im, im,
      _multi_help1< T, T, std::binder2nd< _threshold< T > > >( f ) );
}

template< typename T, typename TPtr >
image< T > threshold( const const_image_ref< T, TPtr > &im, T threshval )
{
  return image_func< T >( im, std::bind2nd( _threshold< T >(),
                                            threshval ) );
}

template< typename T >
image_ref< T > &bilevel_doubleIt( image_ref< T > &im,
                                  T min, T max, T val1, T val2 )
{
  _bilevel_double< T > f( val1, val2, min, max );
  return image_apply( im, im,
                      _multi_help1< T, T, _bilevel_double< T > >( f ) );
}

template< typename T, typename TPtr >
image< T > bilevel_double( const const_image_ref< T, TPtr > &im,
                           T min, T max, T val1, T val2 )
{
  return multi_func< T >( im.rawData(),
                          _bilevel_double< T >( val1, val2, min, max ) );
}

template<typename T, typename TPtr, typename U>
void despeckleKuwahara( const const_image_ref< T, TPtr > &imagein,
                        image< U > &imageout )
{
  if ( &imageout != &imagein )
    imageout.init(imagein.getWidth(),imagein.getHeight());
    
    U intensity;
    float m_ul, m_ur, m_ll, m_lr;
    float v_ul, v_ur, v_ll, v_lr;
    
    for (int y=0; y<imagein.getHeight(); y++)
  for (int x=0; x<imagein.getWidth(); x++)
  {
      if (y<2 || y>imagein.getHeight()-3 || x<2 || x>imagein.getWidth()-3)
      {
         imageout.pixel(x,y)=imagein.pixel(x,y);
    continue;

      }


      // mean
      m_ul=0.0; m_ur=0.0; m_ll=0.0; m_lr=0.0;

      for (int j=0; j<3; j++)
    for (int i=0; i<3; i++)
    {
        m_ul+=(float)imagein.pixel(x-2+i,y-2+j);
        m_ur+=(float)imagein.pixel(x+i,y-2+j);
        m_ll+=(float)imagein.pixel(x-2+i,y+j);
        m_lr+=(float)imagein.pixel(x+i,y+j);
    }

      m_ul=m_ul/9.0;
      m_ur=m_ur/9.0;
      m_ll=m_ll/9.0;
      m_lr=m_lr/9.0;

      
      // estimated variance
      v_ul=0.0; v_ur=0.0; v_ll=0.0; v_lr=0.0;
      
      for (int j=0; j<3; j++)
    for (int i=0; i<3; i++)
    {
        v_ul+=(((float)imagein.pixel(x-2+i,y-2+j)-m_ul)*((float)imagein.pixel(x-2+i,y-2+j)-m_ul));
        v_ur+=(((float)imagein.pixel(x+i,y-2+j)-m_ur)*((float)imagein.pixel(x+i,y-2+j)-m_ur));
        v_ll+=(((float)imagein.pixel(x-2+i,y+j)-m_ll)*((float)imagein.pixel(x-2+i,y+j)-m_ll));
        v_lr+=(((float)imagein.pixel(x+i,y+j)-m_lr)*((float)imagein.pixel(x+i,y+j)-m_lr));
    }

      intensity=(U)m_ul;
      intensity=v_ur<v_ul?(U)m_ur:intensity;
      intensity=v_ll<v_ur?(U)m_ll:intensity;
      intensity=v_lr<v_ll?(U)m_lr:intensity;

      imageout.pixel(x,y)=intensity;
  }  
}


template<typename T>
void despeckleKuwaharaIt( image_ref<T> &imagein )
{
    image<T> tempimage;

    despeckleKuwahara(imagein, tempimage);

    imagein = tempimage;
}

template< typename T, typename TPtr >
image< T > gradSobelX( const const_image_ref< T, TPtr > &im )
{
  const T sobelx1[1][3]={ {-1, 0, 1} };
  const T sobelx2[3][1]={ { 1 },
                          { 2 },
                          { 1 } };
  const_image_ref< T >
    filter1( &sobelx1[0][0], 3, 1 ),
    filter2( &sobelx2[0][0], 1, 3 );
  return correlate( correlate( im, filter1 ), filter2 );
}


template< typename T, typename TPtr >
image< T > gradSobelY( const const_image_ref< T, TPtr > &im )
{
  const T sobely1[3][1]={ { -1 },
                          {  0 },
                          {  1 } };
  const T sobely2[1][3]={ { 1, 2, 1} };
  const_image_ref< T >
    filter1( &sobely1[0][0], 1, 3 ),
    filter2( &sobely2[0][0], 3, 1 );
  return correlate( correlate( im, filter1 ), filter2 );
}

template< typename T, typename TPtr >
image< T > edgeSobelSqr( const const_image_ref< T, TPtr > &im )
{
  image< T >
    gradX( gradSobelX( im ) ),
    gradY( gradSobelY( im ) );
  return gradX * gradX + gradY * gradY;
}

template< typename T, typename TPtr >
image< T > edgeSobelNorm( const const_image_ref< T, TPtr > &im )
{
  return squareRoot( edgeSobelSqr( im ) );
}

// template<typename T, typename U>
// void edgeSobelOrientation(const image<T> &imagein, image<U> &imageout)
// {
//     int  x,y, sx,sy;
//     int  height, width;
//     image<U> tempimage;
//     float  sxIntensity, syIntensity;


//     static float sobelx[3][3]={  { -1.0, -2.0,-1.0},
//                                 {  0.0,  0.0, 0.0},
//                                 {  1.0,  2.0, 1.0}};

//     static float sobely[3][3]={  {-1.0, 0.0, 1.0},
//                                  {-2.0, 0.0, 2.0},
//                                  {-1.0, 0.0, 1.0}};

//     height=imagein.getHeight();
//     width=imagein.getWidth();
//     tempimage.init(width,height);

//     for(y=1;y<height-1;y++)
//         for (x=1;x<width-1;x++)
//         {
//             sxIntensity=0.0; syIntensity=0.0;
//             for(sx=0;sx<3;sx++)
//                 for(sy=0;sy<3;sy++)
//                 {
//                     sxIntensity+=sobelx[sy][sx]*(float)imagein.pixel(x+sx-1,y+sy-1);
//                     syIntensity+=sobely[sy][sx]*(float)imagein.pixel(x+sx-1,y+sy-1);
//                 }
//             if (sxIntensity==0.0)
//               tempimage.pixel(x,y)=0.0;
//             else
//               tempimage.pixel(x,y)=atan(syIntensity/sxIntensity);


//         }

//     imageout = tempimage;
// }

template< typename T, typename TPtr >
image< T > edgeLaplacian( const const_image_ref< T, TPtr > &im )
{
  const T laplacian[3][3]={ {  0, -1,  0 },
                            { -1,  4, -1 },
                            {  0, -1,  0 } };
  const_image_ref< T > f( &laplacian[0][0], 3, 3 );
  return correlate( im, f );
}

template< typename T, typename TPtr >
image< T > edgeLoG( const const_image_ref< T, TPtr > &imagein )
{
  const T LoG[9][9]= { {   0,  0,  3,  2,  2,  2,  3,  0,  0},
                       {   0,  2,  3,  5,  5,  5,  3,  2,  0},
                       {   3,  3,  5,  3,  0,  3,  5,  3,  3},
                       {   2,  5,  3,-12,-23,-12,  3,  5,  2},
                       {   2,  5,  0,-23,-40,-23,  0,  5,  2},
                       {   2,  5,  3,-12,-23,-12,  3,  5,  2},
                       {   3,  3,  5,  3,  0,  3,  5,  3,  3},
                       {   0,  2,  3,  5,  5,  5,  3,  2,  0},
                       {   0,  0,  3,  2,  2,  2,  3,  0,  0} };

  const_image_ref< T > filter( &LoG[0][0], 9, 9 );
  image< T > tempimage( correlate( imagein, filter ) );
    
  /*static float LoG[5][5]={{  0.0, 0.0,-1.0, 0.0, 0.0},
                            {  0.0,-1.0,-2.0,-1.0, 0.0},
                            { -1.0,-2.0,16.0,-2.0,-1.0},
                            {  0.0,-1.0,-2.0,-1.0, 0.0},
                            {  0.0, 0.0,-1.0, 0.0, 0.0}};
  */

  bilevelIt( tempimage, (T)0, (T)255, (T)0 );
  
  image< T > imageout;
  imageout.init( imagein.getWidth(), imagein.getHeight() );

  // Shouldn't actually threshold - must check zero crossings
  for ( int y=4; y<imagein.getHeight()-4; y++ )
    for ( int x=4; x<imagein.getWidth()-4; x++ ) {
      if (tempimage.pixel(x,y)<=127) continue;
      if ( ( tempimage.pixel( x-1, y-1 ) > 127 ) ||
           ( tempimage.pixel( x-1, y   ) > 127 ) ||
           ( tempimage.pixel( x-1, y+1 ) > 127 ) ||
           ( tempimage.pixel( x  , y-1 ) > 127 ) ||
           ( tempimage.pixel( x  , y+1 ) > 127 ) ||
           ( tempimage.pixel( x+1, y-1 ) > 127 ) ||
           ( tempimage.pixel( x+1, y   ) > 127 ) ||
           ( tempimage.pixel( x+1, y+1 ) > 127 ) )
        imageout.pixel( x  , y   ) = 255;
    }
  return imageout;
}

// IMPORTANT NOTE
// imagein1 is image at time t-1
// imagein2 is image at time t
template< typename T, typename T1Ptr, typename T2Ptr >
image< T > edgeHaynesJain( const const_image_ref< T, T1Ptr > &imagein1,
                           const const_image_ref< T, T2Ptr > &imagein2 )
{
    MMERROR( imagein2.getWidth () == imagein1.getWidth () &&
             imagein2.getHeight() == imagein1.getHeight(),
             mimasexception, , "image_funcs::edgeHaynesJain - "
             "size of input images do not match" );

    return absolute<T>( edgeSobelNorm<T>( imagein2 ) * ( imagein2 -
                                                         imagein1 ) );
}


template<typename T, typename TPtr, typename U>
void halfResolution( const const_image_ref< T, TPtr > &imagein,
                     image_ref<U> &imageout)
{
    // note that this is only an estimate!
    static float gauss[3][3]={  {  1.0,  2.0, 1.0},
                               {  2.0,  4.0, 2.0},
                               {  1.0,  2.0, 1.0}};
    static float gaussfac=(1.0/16.0);

    assert( &imageout != &imagein );
    imageout.init(imagein.getWidth()/2, imagein.getHeight()/2);

    for (int j=1; j<imagein.getHeight(); j+=2)
        for (int i=1; i<imagein.getWidth(); i+=2)
        {
            float tempval=0.0;
            for (int y=-1; y<=1; y++)
                for (int x=-1; x<=1; x++)
                    tempval+=(imagein.pixel(i+x,j+y)*gauss[y+1][x+1]);
            imageout.pixel(i/2,j/2)=(U)(tempval*gaussfac);
        }
}


template<typename T>
void halfResolutionIt( image_ref<T> &imagein )
{
    image<T> tempimage;

    halfResolution(imagein, tempimage);

    imagein = tempimage;
}

// Note this assumes black edges on white background unlike halfResolution
// above
template<typename T, typename TPtr, typename U>
void halfResolutionEdgePyramid( const const_image_ref< T, TPtr > &imagein,
                                image<U> &imageout)
{
  // set background colour of object to be white
  assert( &imageout != &imagein );
  imageout.init(imagein.getWidth()/2, imagein.getHeight()/2);
  imageout.fill(255);

    for (int j=1; j<imagein.getHeight(); j+=2)
        for (int i=1; i<imagein.getWidth(); i+=2)
        {
            int tempval=0;
      for (int y=-1; y<=1; y++)
                for (int x=-1; x<=1; x++)
                    if (imagein.getPixel(i+x,j+y)>127)
      tempval=255;
      
      imageout.pixel(i/2,j/2)=(U)tempval;
        }
}


template<typename T>
void halfResolutionEdgePyramidIt( image_ref<T> &imagein )
{
    image<T> tempimage;

    halfResolutionEdgePyramid(imagein, tempimage);

    imagein = tempimage;
}

// 3x3 median despecking filter
template< typename T, typename TPtr >
void despeckleMedian( const const_image_ref<T,TPtr> &imagein,
                      image<T> &imageout)
{
  std::vector<T> pixvals;
  assert( &imageout != &imagein );
  imageout.init(imagein.getWidth(),imagein.getHeight());
 
  for (int j=1; j<imagein.getHeight()-1; j++)
  for (int i=1; i<imagein.getWidth()-1; i++)
  {
      if (!pixvals.empty()) pixvals.clear();
      for (int y=-1; y<=1; y++)
    for (int x=-1; x<=1; x++)
        pixvals.push_back(imagein.pixel(i+x,j+y));
      
      sort(pixvals.begin(), pixvals.end());
      imageout.pixel(i,j)=pixvals[pixvals.size()/2];
  }
 
    for (int j=0; j<imagein.getHeight(); j++)
    {
      imageout.pixel(0,j) = imagein.pixel(0,j);
      imageout.pixel(imagein.getWidth()-1,j)=imagein.pixel(imagein.getWidth()-1,j);
    }
    
    for (int i=0; i<imagein.getWidth(); i++)
    {
      imageout.pixel(i,0)=imagein.pixel(i,0);
      imageout.pixel(i,imagein.getHeight()-1)=imagein.pixel(i,imagein.getHeight()-1);
    }
}
 
template<typename T>
void despeckleMedianIt( image_ref<T> &imagein )
{
    image<T> tempimage;

    despeckleMedian(imagein, tempimage);

    imagein = tempimage;
}

template<typename T>
void chamfer(image_ref<T> &inimage)
{
    int nw,n,ne,w,e,sw,s,se;
    int fp;

    for (int j=1; j<inimage.getHeight()-1; j++)
        for (int i=1; i<inimage.getWidth(); i++)
        {
            fp=inimage.pixel(i,j);
            nw=4+inimage.pixel(i-1,j-1);
            n=3+inimage.pixel(i,j-1);
            w=3+inimage.pixel(i-1,j);
            sw=4+inimage.pixel(i-1,j+1);

            fp=fp<nw?fp:nw;
            fp=fp<n?fp:n;
            fp=fp<w?fp:w;
            fp=fp<sw?fp:sw;
            inimage.pixel(i,j)=fp;

        }

    for (int j=inimage.getHeight()-2; j>0; j--)
        for (int i=inimage.getWidth()-2; i>=0; i--)
        {
            fp=inimage.pixel(i,j);
            ne=4+inimage.pixel(i+1,j-1);
            e=3+inimage.pixel(i+1,j);
            s=3+inimage.pixel(i,j+1);
            se=4+inimage.pixel(i+1,j+1);

            fp=fp<ne?fp:ne;
            fp=fp<e?fp:e;
            fp=fp<s?fp:s;
            fp=fp<se?fp:se;
            inimage.pixel(i,j)=fp;

        }
    
    for (int i=0; i<inimage.getWidth(); i++)
    {
      inimage.pixel(i,0)=inimage.pixel(i,1);
      inimage.pixel(i,inimage.getHeight()-1)=inimage.pixel(i,inimage.getHeight()-2);
    }
    
    for (int j=0; j<inimage.getHeight(); j++)
    {
      inimage.pixel(0,j)=inimage.pixel(1,j);
      inimage.pixel(inimage.getWidth()-1,j)=inimage.pixel(inimage.getWidth()-2,j);
    }
 }


template<typename T,typename TPtr>
void paethRotateRadian(const const_image_ref<T,TPtr> &imagein,
                       image<T> &imageout,
                       float radian)
{
  assert( &imageout != &imagein );
  imageout.init(imagein.getWidth(),imagein.getHeight());
  imageout.fill(255);

    if (fabs(radian-M_PI)<1e-3) radian=radian<M_PI?M_PI-1e-3:M_PI+1e-3; // Quick hack to fix the problem when tan goes to infinity
    float s1=-tan(-radian/2.0);
    float s2=sin(-radian);

    float xr, yr;
    int x,y;
    for (int j=0; j<imageout.getHeight(); j++)
        for (int i=0; i<imageout.getWidth(); i++)
        {
            x=i-imageout.getWidth()/2;
            y=j-imageout.getHeight()/2;
            //x=i;
            //y=j;

            xr = (float)x + s1 * (float)y;
            yr = (float)y + s2 * xr;
            xr = xr + s1 * yr;

            xr+=(float)imageout.getWidth()/2;
            yr+=(float)imageout.getHeight()/2;
            imageout.pixel(i, j)=imagein.getPixel((int)xr,(int)yr);
            //imageout.setPixel(i,j, imagein.getPixel(i,j));
        }
}

template< typename T, typename T1Ptr, typename T2Ptr >
void multiLevelCorrelation( const const_image_ref< T, T1Ptr > &edgemap,
                            const const_image_ref< T, T2Ptr > &objectedgemap,
                            int levels, int &x, int &y, float &corrval,
                            bool chamferObject)
{
    image<T> object[levels], image[levels];
    
    image[0] = edgemap;
   // invert(image[0]);

    object[0] = objectedgemap;
   // invert(object[0]);
  
     for (int l=1; l<levels; l++)
    {
        halfResolutionEdgePyramid(image[l-1],image[l]);
        halfResolutionEdgePyramid(object[l-1],object[l]);
    }
    
    for (int l=0; l<levels; l++)
    {
      image[l] = (T)255 - image[l];
      object[l] = (T)255 - object[l];
  //image[l].display();
  //object[l].display();
  chamfer(image[l]);
  if (chamferObject) chamfer(object[l]);
    }
 
     int xmin = INT_MIN, ymin = INT_MIN;
    int xminsearch, xmaxsearch, yminsearch, ymaxsearch;
    float corr, mincorr;

    x=0; y=0;
    xminsearch=0;
    yminsearch=0;
    xmaxsearch=image[levels-1].getWidth();
    ymaxsearch=image[levels-1].getHeight();

    for (int l=levels-1; l>=0; l--)
    {
  //object[l].display();
  //image[l].display();

        mincorr=FLT_MAX;
        for (int j=yminsearch; j<ymaxsearch-object[l].getHeight(); j++)
            for (int i=xminsearch; i<xmaxsearch-object[l].getWidth(); i++)
            {
                corr=0.0;
    int count=0;
                for (int b=0; b<object[l].getHeight(); b++)
                    for (int a=0; a<object[l].getWidth(); a++)
                    {
                        if (!chamferObject) if (object[l].pixel(a,b)>127) continue;
                        corr+=(float)abs((int)image[l].pixel(i+a,j+b)-(int)object[l].pixel(a,b));
      count++;
                    }

    corr/=count;
                if (corr<mincorr)
                {
                    xmin=i; ymin=j;
                    mincorr=corr;
                }
            }

        if (l!=0)
        {
            xminsearch=(xmin*2)-1-20;
            yminsearch=(ymin*2)-1-20;

            xmaxsearch=(xmin*2)-1+object[l-1].getWidth()+20;
            ymaxsearch=(ymin*2)-1+object[l-1].getHeight()+20;

            if (xminsearch<0) xminsearch=0;
            if (yminsearch<0) yminsearch=0;

            if (xmaxsearch>image[l-1].getWidth()) xmaxsearch=image[l-1].getWidth();
            if (ymaxsearch>image[l-1].getHeight()) ymaxsearch=image[l-1].getHeight();
        }


    }

    assert( xmin != INT_MIN && ymin != INT_MIN );
    x=xmin;
    y=ymin;
    corrval=mincorr;

}

template< typename T, typename T1Ptr, typename T2Ptr >
void correlation( const const_image_ref<T> &imagein,
                  const const_image_ref<T> &templatein,
                  int &x, int &y, float &corrval)
{

    int xmin = INT_MIN, ymin = INT_MIN;
    int xminsearch, xmaxsearch, yminsearch, ymaxsearch;
    float corr, mincorr;

    x=0; y=0;
    xminsearch=0;
    yminsearch=0;
    xmaxsearch=imagein.getWidth();
    ymaxsearch=imagein.getHeight();

    mincorr=FLT_MAX;
    for (int j=yminsearch; j<ymaxsearch-templatein.getHeight(); j++)
        for (int i=xminsearch; i<xmaxsearch-templatein.getWidth(); i++)
        {
            corr=0.0;
            for (int b=0; b<templatein.getHeight(); b++)
                for (int a=0; a<templatein.getWidth(); a++)
                {
                    //if (object[l].getPixel(a,b)>127) continue;
                    corr+=(float)abs((int)imagein.pixel(i+a,j+b)-(int)templatein.pixel(a,b));
                }

            if (corr<mincorr)
            {
                xmin=i; ymin=j;
                mincorr=corr;
            }
        }

    assert( xmin != INT_MIN && ymin != INT_MIN );
    x=xmin;
    y=ymin;
    corrval=mincorr;

}

template<typename T, typename T1Ptr, typename T2Ptr >
void correlationChamfer( const const_image_ref< T, T1Ptr > &imagein,
                         const const_image_ref< T, T2Ptr > &templatein,
                        int &x, int &y, float &corrval)
{

    int xmin = INT_MIN, ymin = INT_MIN;
    int xminsearch, xmaxsearch, yminsearch, ymaxsearch;
    float corr, mincorr;

    x=0; y=0;
    xminsearch=0;
    yminsearch=0;
    xmaxsearch=imagein.getWidth();
    ymaxsearch=imagein.getHeight();

    mincorr=FLT_MAX;
    for (int j=yminsearch; j<ymaxsearch-templatein.getHeight(); j++)
        for (int i=xminsearch; i<xmaxsearch-templatein.getWidth(); i++)
        {
            corr=0.0;
            for (int b=0; b<templatein.getHeight(); b++)
                for (int a=0; a<templatein.getWidth(); a++)
                {
        if (templatein.pixel(a,b)>127) continue;
                    corr+=(float)abs((int)imagein.pixel(i+a,j+b)-(int)templatein.pixel(a,b));
                }

            if (corr<mincorr)
            {
                xmin=i; ymin=j;
                mincorr=corr;
            }
        }

    assert( xmin != INT_MIN && ymin != INT_MIN );
    x=xmin;
    y=ymin;
    corrval=mincorr;

}

template<typename T>
void rotationChamferCorrelationEdgeMap(image<T> &edgemap, image<T>
  &objectedgemap, int angle1, int angle2, int anglestep, int &x, int &y, int &outangle, float &corrval)
{
    image<T> chamferedImage, rotatedObject, chamferedObject;

    chamferedImage = (T)255 - edgemap;
    chamfer(chamferedImage);

  
    objectedgemap = (T)255 - objectedgemap;
    // Make sure that object boundaries are white:
    drawLine(objectedgemap, 0,0,0,objectedgemap.getHeight()-1,255);
    drawLine(objectedgemap, objectedgemap.getWidth()-1,0,objectedgemap.getWidth()-1,objectedgemap.getHeight()-1,255);
    drawLine(objectedgemap, 0,0,objectedgemap.getWidth()-1,objectedgemap.getHeight()-1,255);
    drawLine(objectedgemap, 0,objectedgemap.getHeight()-1,objectedgemap.getWidth()-1,objectedgemap.getHeight()-1,255);
    //chamfer(objectedgemap);

    float cval;
    int tx,ty;
    corrval=FLT_MAX;
    for (int angle=angle1; angle<=angle2; angle+=anglestep)
    {
   rotatedObject = objectedgemap;
  //invert(rotatedObject);
  paethRotateRadian(rotatedObject,chamferedObject,(float)angle*M_PI/180.0);
  //chamfer(chamferedObject);
  correlationChamfer(chamferedImage, chamferedObject, tx, ty, cval);

  if (cval<corrval)
  {
      corrval=cval;
      x=tx;
      y=ty;
      outangle=angle;
  }
    }
}

template< typename T, typename TPtr >
image< T > focusEnhance( const const_image_ref< T, TPtr > &im )
{
  const T enhance[3][3]={ { -1,  0, -1 },
                          {  0,  7,  0 },
                          { -1,  0, -1 } };
  const_image_ref< T > filter( &enhance[0][0], 3, 3 );
  return normalise( correlate( im, filter ), (T)0, (T)255 );
}

template< typename T, typename TPtr >
image< T > featureEnhance( const const_image_ref< T, TPtr > &im )
{
  const T enhance[3][3]={ {  0, -1,  0 },
                          { -1, 10, -1 },
                          {  0, -1,  0 } };
  const_image_ref< T > filter( &enhance[0][0], 3, 3 );
  return normalise( correlate( im, filter ), (T)0, (T)255 );
}

template< typename T, typename TPtr >
image< T > emboss( const const_image_ref< T, TPtr > &im )
{
  const T diff[3][3]={ { -1,  0,  0 },
                       {  0,  0,  0 },
                       {  0,  0,  1 } };
  const_image_ref< T > filter( &diff[0][0], 3, 3 );
  return normalise( correlate( im, filter ), (T)0, (T)255 );
}

template< typename T, typename TPtr >
image< T > softenHeavy( const const_image_ref< T, TPtr > &im )
{
  const T median[3]={ 1, 1, 1 };
  const_image_ref< T >
    filterx( &median[0], 3, 1 ),
    filtery( &median[0], 1, 3 );
 return
   normalise( correlate( correlate( im, filterx ), filtery ),
              (T)0, (T)255 );
}

template< typename T, typename TPtr >
image< T > softenMedium( const const_image_ref< T, TPtr > &im )
{
  const T soften[3][3]={ {  1,  1,  1 },
                         {  1,  2,  1 },
                         {  1,  1,  1 } };
  const_image_ref< T > filter( &soften[0][0], 3, 3 );
  return normalise( correlate( im, filter ), (T)0, (T)255 );
}

template< typename T, typename TPtr >
image< T > softenLight( const const_image_ref< T, TPtr > &im )
{
  const T soften[3][3]={ {  6,  12,  6 },
                         { 12,  25, 12 },
                         {  6,  12,  6 } };
  const_image_ref< T > filter( &soften[0][0], 3, 3 );
  return normalise( correlate( im, filter ), (T)0, (T)255 );
}

template< typename T, typename TPtr >
void oilPainting( const const_image_ref< T, TPtr > &imagein,
                  image<T> &imageout,
                  int regiondim)
{
   const int fy=regiondim;
    const int fx=regiondim;

    assert( &imageout != &imagein );
    imageout.init(imagein.getWidth(),imagein.getHeight());
    
    for (int y=0; y<imagein.getHeight(); y++)
  for (int x=0; x<imagein.getWidth(); x++)
  {
      if (y-fy/2<0 || y+fy/2>imagein.getHeight()-1 || x-fx/2<0 || x+fx/2>imagein.getWidth()-1)
      {
         imageout.pixel(x,y)=imagein.pixel(x,y);
    continue;
      }

       std::vector<int> pixcolor(regiondim*regiondim);
       std::vector<int> graylevel(256);
    
      int count=0;
      for (int j=-fy/2; j<=fy/2; j++)
    for (int i=-fx/2; i<=fx/2; i++)
    {
        pixcolor[count]=imagein.pixel(x+i,y+j);
        graylevel[pixcolor[count]]+=1;
        
        count++;
    }
      
      count=0;
      int intensity=0, maxuse=0;
      for (int j=-fy/2; j<=fy/2; j++)
    for (int i=-fx/2; i<=fx/2; i++)
    {
        if (graylevel[pixcolor[count]]>maxuse)
        {
      maxuse=graylevel[pixcolor[count]];
      intensity=pixcolor[count];
        }

        count++;
    }
      
      imageout.pixel(x,y)=(T)(intensity);
  }

    
}

#ifdef HAVE_OPENGL
template<typename T,typename TPtr>
void displayWithOverlay( image_output< rgba< unsigned char > > &display,
                         const const_image_ref<T,TPtr>& im,
                         const property_image& pim)
{
  image< rgba< unsigned char > > output;
  output.init( im.getWidth(), im.getHeight() );
    
  for( int j=0;j<im.getHeight();j++)
    for( int i=0;i<im.getWidth();i++) {
      prop_ptr p( pim.pixel(i,j) );
      if( p )
  p->draw( output, i, j, false);
    }
    
  for( int j=0;j<im.getHeight();j++)
    for( int i=0;i<im.getWidth();i++)
      if ( (int)output.pixel( i, j ) ==  0 ) {
  unsigned char value = (unsigned char)im.pixel( i, j );
  output.pixel( i, j ) =
    rgba< unsigned char >( value, value, value );
      };
          
  display << output;
}
#endif


//some common features here with vector class ...
//a cleverer man could probably use the same underlying base class
//but for now it doesnt matter
 template< typename T, typename T1Ptr, typename T2Ptr >
double dotProduct( const const_image_ref<T,T1Ptr> &im1,
                   const const_image_ref<T,T2Ptr> &im2 ) throw (mimasexception)
{
  MMERROR( im1.getWidth() == im2.getWidth() &&
           im1.getHeight() == im2.getHeight(),
           mimasexception, ,
           "Can't (won't) compute the dot product of two images of "
           "differing sizes." );

  double result = 0.0;

  const int size = (signed)im1.getSize();
  const T
    *p = im1.rawData(),
    *q = im2.rawData();

  for ( int i=0; i<size; i++ )
    result += *p++ * *q++;

  return result;
}

template<
  typename T1, typename T2,
  typename TAPtr, typename TBPtr, typename TCPtr
>
image< bool > nonMaximaSuppression
  ( const const_image_ref< T1, TAPtr > &gradientX,
    const const_image_ref< T1, TBPtr > &gradientY,
    const const_image_ref< T2, TCPtr > &gradientSqr,
    T2 threshold );

}

#include "image_funcs.tcc"

#endif

---