//// 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 <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"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 );} elseretVal = 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 );} elseretVal = 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 ) );} elsereturn 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 histogramint 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;}// meanm_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 variancev_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 );}/** Laplacian of Gaussian filter.@todo Determine zero-crossings to have a complete edge-detector. */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 crossingsfor ( 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 ttemplate< 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// abovetemplate<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 whiteassert( &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 filtertemplate< 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 infinityfloat 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 );returnnormalise( 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);}}//some common features here with vector class ...//a cleverer man could probably use the same underlying base class//but for now it doesnt mattertemplate< 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;}/** Fast non-maxima suppression.This function provides fast non-maxima suppression.The function is provided separately, so that one can use theSobel-operator as well as the Gauss-gradient together with thismethod.@author Jan Wedekind (jan@wedesoft.de)@date Mon Apr 10 13:53:00 2006 */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