Rev Author Line No. Line
178 kaklik 1 //
2 // Generic image functions
3 //
4 // Bala Amavasai (bala@amavasai.org)
5 // Mon Oct 1 15:41:44 BST 2001
6 //
7 // $Header: /cvs/mimas2/include/image_funcs.h,v 1.1.1.1 2005/08/09 15:37:45 engmb Exp $
8 //
9  
10 #ifndef IMAGE_FUNCS_H
11 #define IMAGE_FUNCS_H
12  
13 #include <cassert>
14 #include <float.h>
15 #include <iostream>
16 #include <vector>
17 #include <algorithm>
18 #include <cstdlib>
19 #include <values.h>
20 #include <cmath>
21 #include "image.h"
22 #include "image_conv.h"
23 #include "image_op.h"
24 #include "rgba.h"
25  
26 namespace mimas {
27  
28 template< typename T, typename TPtr >
29 T min_val( const const_image_ref< T, TPtr > &image )
30 {
31 T retVal;
32 if ( image.initialised() ) {
33 int size = image.getWidth() * image.getHeight();
34 retVal = *std::min_element( image.rawData(),
35 image.rawData() + size );
36 } else
37 retVal = T();
38 return retVal;
39 }
40  
41 template< typename T, typename TPtr >
42 T max_val( const const_image_ref< T, TPtr > &image )
43 {
44 T retVal;
45 if ( image.initialised() ) {
46 int size = image.getWidth() * image.getHeight();
47 retVal = *std::max_element( image.rawData(),
48 image.rawData() + size );
49 } else
50 retVal = T();
51 return retVal;
52 }
53  
54 template< typename T >
55 image_ref< T > &normaliseIt( image_ref< T > &im,
56 const T &val1, const T &val2 )
57 {
58 _normalise< T > f( min_val( im ), max_val( im ), val1, val2 );
59 return image_apply( im, im,
60 _multi_help1< T, T, _normalise< T > >( f ) );
61 }
62  
63 template< typename T, typename TPtr >
64 image< T > normalise( const const_image_ref< T, TPtr > &im,
65 const T &val1, const T &val2 )
66 {
67 if ( im.initialised() ) {// Image must not be empty.
68 return image_func< T >( im,
69 _normalise< T >( min_val( im ),
70 max_val( im ), val1, val2 ) );
71 } else
72 return image< T >();
73 }
74  
75 template<typename T>
76 void equaliseIt( image_ref<T> &imagein )
77 {
78  
79 int H[256], TR[256];
80  
81 for (int count=0; count<256; count++)
82 H[count]=0;
83  
84 int intensity;
85 for (int j=0; j<imagein.getHeight(); j++)
86 for (int i=0; i<imagein.getWidth(); i++)
87 {
88 intensity=imagein.pixel(i,j);
89 assert( intensity <= 255 );
90 H[intensity]+=1;
91 }
92  
93 // convert to cumulative histogram
94 int imagesize=imagein.getWidth()*imagein.getHeight();
95 for (int count=1; count<256; count++)
96 {
97 H[count]=H[count]+H[count-1];
98 TR[count-1]=(255*H[count-1])/imagesize;
99 }
100 TR[255]=(255*H[255])/imagesize;
101  
102 for (int j=0; j<imagein.getHeight(); j++)
103 for (int i=0; i<imagein.getWidth(); i++)
104 imagein.pixel(i,j)=(T)TR[imagein.pixel(i,j)];
105 }
106  
107 template< typename T >
108 image_ref< T > &bilevelIt( image_ref< T > &im, T threshval,
109 T val1, T val2 )
110 {
111 std::binder2nd< _bilevel< T > > f( _bilevel< T >( val1, val2 ),
112 threshval );
113 return image_apply
114 ( im, im,
115 _multi_help1< T, T, std::binder2nd< _bilevel< T > > >( f ) );
116 }
117  
118 template< typename T, typename TPtr >
119 image< T > bilevel( const const_image_ref< T, TPtr > &im,
120 T threshval, T val1, T val2 )
121 {
122 return image_func< T >( im,
123 std::bind2nd( _bilevel< T >( val1, val2 ),
124 threshval ) );
125 }
126  
127 template< typename T >
128 image_ref< T > &thresholdIt( image_ref< T > &im, T threshval )
129 {
130 std::binder2nd< _threshold< T > > f( _threshold< T >(), threshval );
131 return image_apply
132 ( im, im,
133 _multi_help1< T, T, std::binder2nd< _threshold< T > > >( f ) );
134 }
135  
136 template< typename T, typename TPtr >
137 image< T > threshold( const const_image_ref< T, TPtr > &im, T threshval )
138 {
139 return image_func< T >( im, std::bind2nd( _threshold< T >(),
140 threshval ) );
141 }
142  
143 template< typename T >
144 image_ref< T > &bilevel_doubleIt( image_ref< T > &im,
145 T min, T max, T val1, T val2 )
146 {
147 _bilevel_double< T > f( val1, val2, min, max );
148 return image_apply( im, im,
149 _multi_help1< T, T, _bilevel_double< T > >( f ) );
150 }
151  
152 template< typename T, typename TPtr >
153 image< T > bilevel_double( const const_image_ref< T, TPtr > &im,
154 T min, T max, T val1, T val2 )
155 {
156 return multi_func< T >( im.rawData(),
157 _bilevel_double< T >( val1, val2, min, max ) );
158 }
159  
160 template<typename T, typename TPtr, typename U>
161 void despeckleKuwahara( const const_image_ref< T, TPtr > &imagein,
162 image< U > &imageout )
163 {
164 if ( &imageout != &imagein )
165 imageout.init(imagein.getWidth(),imagein.getHeight());
166  
167 U intensity;
168 float m_ul, m_ur, m_ll, m_lr;
169 float v_ul, v_ur, v_ll, v_lr;
170  
171 for (int y=0; y<imagein.getHeight(); y++)
172 for (int x=0; x<imagein.getWidth(); x++)
173 {
174 if (y<2 || y>imagein.getHeight()-3 || x<2 || x>imagein.getWidth()-3)
175 {
176 imageout.pixel(x,y)=imagein.pixel(x,y);
177 continue;
178  
179 }
180  
181  
182 // mean
183 m_ul=0.0; m_ur=0.0; m_ll=0.0; m_lr=0.0;
184  
185 for (int j=0; j<3; j++)
186 for (int i=0; i<3; i++)
187 {
188 m_ul+=(float)imagein.pixel(x-2+i,y-2+j);
189 m_ur+=(float)imagein.pixel(x+i,y-2+j);
190 m_ll+=(float)imagein.pixel(x-2+i,y+j);
191 m_lr+=(float)imagein.pixel(x+i,y+j);
192 }
193  
194 m_ul=m_ul/9.0;
195 m_ur=m_ur/9.0;
196 m_ll=m_ll/9.0;
197 m_lr=m_lr/9.0;
198  
199  
200 // estimated variance
201 v_ul=0.0; v_ur=0.0; v_ll=0.0; v_lr=0.0;
202  
203 for (int j=0; j<3; j++)
204 for (int i=0; i<3; i++)
205 {
206 v_ul+=(((float)imagein.pixel(x-2+i,y-2+j)-m_ul)*((float)imagein.pixel(x-2+i,y-2+j)-m_ul));
207 v_ur+=(((float)imagein.pixel(x+i,y-2+j)-m_ur)*((float)imagein.pixel(x+i,y-2+j)-m_ur));
208 v_ll+=(((float)imagein.pixel(x-2+i,y+j)-m_ll)*((float)imagein.pixel(x-2+i,y+j)-m_ll));
209 v_lr+=(((float)imagein.pixel(x+i,y+j)-m_lr)*((float)imagein.pixel(x+i,y+j)-m_lr));
210 }
211  
212 intensity=(U)m_ul;
213 intensity=v_ur<v_ul?(U)m_ur:intensity;
214 intensity=v_ll<v_ur?(U)m_ll:intensity;
215 intensity=v_lr<v_ll?(U)m_lr:intensity;
216  
217 imageout.pixel(x,y)=intensity;
218 }
219 }
220  
221  
222 template<typename T>
223 void despeckleKuwaharaIt( image_ref<T> &imagein )
224 {
225 image<T> tempimage;
226  
227 despeckleKuwahara(imagein, tempimage);
228  
229 imagein = tempimage;
230 }
231  
232 template< typename T, typename TPtr >
233 image< T > gradSobelX( const const_image_ref< T, TPtr > &im )
234 {
235 const T sobelx1[1][3]={ {-1, 0, 1} };
236 const T sobelx2[3][1]={ { 1 },
237 { 2 },
238 { 1 } };
239 const_image_ref< T >
240 filter1( &sobelx1[0][0], 3, 1 ),
241 filter2( &sobelx2[0][0], 1, 3 );
242 return correlate( correlate( im, filter1 ), filter2 );
243 }
244  
245  
246 template< typename T, typename TPtr >
247 image< T > gradSobelY( const const_image_ref< T, TPtr > &im )
248 {
249 const T sobely1[3][1]={ { -1 },
250 { 0 },
251 { 1 } };
252 const T sobely2[1][3]={ { 1, 2, 1} };
253 const_image_ref< T >
254 filter1( &sobely1[0][0], 1, 3 ),
255 filter2( &sobely2[0][0], 3, 1 );
256 return correlate( correlate( im, filter1 ), filter2 );
257 }
258  
259 template< typename T, typename TPtr >
260 image< T > edgeSobelSqr( const const_image_ref< T, TPtr > &im )
261 {
262 image< T >
263 gradX( gradSobelX( im ) ),
264 gradY( gradSobelY( im ) );
265 return gradX * gradX + gradY * gradY;
266 }
267  
268 template< typename T, typename TPtr >
269 image< T > edgeSobelNorm( const const_image_ref< T, TPtr > &im )
270 {
271 return squareRoot( edgeSobelSqr( im ) );
272 }
273  
274 // template<typename T, typename U>
275 // void edgeSobelOrientation(const image<T> &imagein, image<U> &imageout)
276 // {
277 // int x,y, sx,sy;
278 // int height, width;
279 // image<U> tempimage;
280 // float sxIntensity, syIntensity;
281  
282  
283 // static float sobelx[3][3]={ { -1.0, -2.0,-1.0},
284 // { 0.0, 0.0, 0.0},
285 // { 1.0, 2.0, 1.0}};
286  
287 // static float sobely[3][3]={ {-1.0, 0.0, 1.0},
288 // {-2.0, 0.0, 2.0},
289 // {-1.0, 0.0, 1.0}};
290  
291 // height=imagein.getHeight();
292 // width=imagein.getWidth();
293 // tempimage.init(width,height);
294  
295 // for(y=1;y<height-1;y++)
296 // for (x=1;x<width-1;x++)
297 // {
298 // sxIntensity=0.0; syIntensity=0.0;
299 // for(sx=0;sx<3;sx++)
300 // for(sy=0;sy<3;sy++)
301 // {
302 // sxIntensity+=sobelx[sy][sx]*(float)imagein.pixel(x+sx-1,y+sy-1);
303 // syIntensity+=sobely[sy][sx]*(float)imagein.pixel(x+sx-1,y+sy-1);
304 // }
305 // if (sxIntensity==0.0)
306 // tempimage.pixel(x,y)=0.0;
307 // else
308 // tempimage.pixel(x,y)=atan(syIntensity/sxIntensity);
309  
310  
311 // }
312  
313 // imageout = tempimage;
314 // }
315  
316 template< typename T, typename TPtr >
317 image< T > edgeLaplacian( const const_image_ref< T, TPtr > &im )
318 {
319 const T laplacian[3][3]={ { 0, -1, 0 },
320 { -1, 4, -1 },
321 { 0, -1, 0 } };
322 const_image_ref< T > f( &laplacian[0][0], 3, 3 );
323 return correlate( im, f );
324 }
325  
326 /** Laplacian of Gaussian filter.
327 @todo Determine zero-crossings to have a complete edge-detector. */
328 template< typename T, typename TPtr >
329 image< T > edgeLoG( const const_image_ref< T, TPtr > &imagein )
330 {
331 const T LoG[9][9]= { { 0, 0, 3, 2, 2, 2, 3, 0, 0},
332 { 0, 2, 3, 5, 5, 5, 3, 2, 0},
333 { 3, 3, 5, 3, 0, 3, 5, 3, 3},
334 { 2, 5, 3,-12,-23,-12, 3, 5, 2},
335 { 2, 5, 0,-23,-40,-23, 0, 5, 2},
336 { 2, 5, 3,-12,-23,-12, 3, 5, 2},
337 { 3, 3, 5, 3, 0, 3, 5, 3, 3},
338 { 0, 2, 3, 5, 5, 5, 3, 2, 0},
339 { 0, 0, 3, 2, 2, 2, 3, 0, 0} };
340  
341 const_image_ref< T > filter( &LoG[0][0], 9, 9 );
342 image< T > tempimage( correlate( imagein, filter ) );
343  
344 /*static float LoG[5][5]={{ 0.0, 0.0,-1.0, 0.0, 0.0},
345 { 0.0,-1.0,-2.0,-1.0, 0.0},
346 { -1.0,-2.0,16.0,-2.0,-1.0},
347 { 0.0,-1.0,-2.0,-1.0, 0.0},
348 { 0.0, 0.0,-1.0, 0.0, 0.0}};
349 */
350  
351 bilevelIt( tempimage, (T)0, (T)255, (T)0 );
352  
353 image< T > imageout;
354 imageout.init( imagein.getWidth(), imagein.getHeight() );
355  
356 // Shouldn't actually threshold - must check zero crossings
357 for ( int y=4; y<imagein.getHeight()-4; y++ )
358 for ( int x=4; x<imagein.getWidth()-4; x++ ) {
359 if (tempimage.pixel(x,y)<=127) continue;
360 if ( ( tempimage.pixel( x-1, y-1 ) > 127 ) ||
361 ( tempimage.pixel( x-1, y ) > 127 ) ||
362 ( tempimage.pixel( x-1, y+1 ) > 127 ) ||
363 ( tempimage.pixel( x , y-1 ) > 127 ) ||
364 ( tempimage.pixel( x , y+1 ) > 127 ) ||
365 ( tempimage.pixel( x+1, y-1 ) > 127 ) ||
366 ( tempimage.pixel( x+1, y ) > 127 ) ||
367 ( tempimage.pixel( x+1, y+1 ) > 127 ) )
368 imageout.pixel( x , y ) = 255;
369 }
370 return imageout;
371 }
372  
373 // IMPORTANT NOTE
374 // imagein1 is image at time t-1
375 // imagein2 is image at time t
376 template< typename T, typename T1Ptr, typename T2Ptr >
377 image< T > edgeHaynesJain( const const_image_ref< T, T1Ptr > &imagein1,
378 const const_image_ref< T, T2Ptr > &imagein2 )
379 {
380 MMERROR( imagein2.getWidth () == imagein1.getWidth () &&
381 imagein2.getHeight() == imagein1.getHeight(),
382 mimasexception, , "image_funcs::edgeHaynesJain - "
383 "size of input images do not match" );
384  
385 return absolute<T>( edgeSobelNorm<T>( imagein2 ) * ( imagein2 -
386 imagein1 ) );
387 }
388  
389  
390 template<typename T, typename TPtr, typename U>
391 void halfResolution( const const_image_ref< T, TPtr > &imagein,
392 image_ref<U> &imageout)
393 {
394 // note that this is only an estimate!
395 static float gauss[3][3]={ { 1.0, 2.0, 1.0},
396 { 2.0, 4.0, 2.0},
397 { 1.0, 2.0, 1.0}};
398 static float gaussfac=(1.0/16.0);
399  
400 assert( &imageout != &imagein );
401 imageout.init(imagein.getWidth()/2, imagein.getHeight()/2);
402  
403 for (int j=1; j<imagein.getHeight(); j+=2)
404 for (int i=1; i<imagein.getWidth(); i+=2)
405 {
406 float tempval=0.0;
407 for (int y=-1; y<=1; y++)
408 for (int x=-1; x<=1; x++)
409 tempval+=(imagein.pixel(i+x,j+y)*gauss[y+1][x+1]);
410 imageout.pixel(i/2,j/2)=(U)(tempval*gaussfac);
411 }
412 }
413  
414  
415 template<typename T>
416 void halfResolutionIt( image_ref<T> &imagein )
417 {
418 image<T> tempimage;
419  
420 halfResolution(imagein, tempimage);
421  
422 imagein = tempimage;
423 }
424  
425 // Note this assumes black edges on white background unlike halfResolution
426 // above
427 template<typename T, typename TPtr, typename U>
428 void halfResolutionEdgePyramid( const const_image_ref< T, TPtr > &imagein,
429 image<U> &imageout)
430 {
431 // set background colour of object to be white
432 assert( &imageout != &imagein );
433 imageout.init(imagein.getWidth()/2, imagein.getHeight()/2);
434 imageout.fill(255);
435  
436 for (int j=1; j<imagein.getHeight(); j+=2)
437 for (int i=1; i<imagein.getWidth(); i+=2)
438 {
439 int tempval=0;
440 for (int y=-1; y<=1; y++)
441 for (int x=-1; x<=1; x++)
442 if (imagein.getPixel(i+x,j+y)>127)
443 tempval=255;
444  
445 imageout.pixel(i/2,j/2)=(U)tempval;
446 }
447 }
448  
449  
450 template<typename T>
451 void halfResolutionEdgePyramidIt( image_ref<T> &imagein )
452 {
453 image<T> tempimage;
454  
455 halfResolutionEdgePyramid(imagein, tempimage);
456  
457 imagein = tempimage;
458 }
459  
460 // 3x3 median despecking filter
461 template< typename T, typename TPtr >
462 void despeckleMedian( const const_image_ref<T,TPtr> &imagein,
463 image<T> &imageout)
464 {
465 std::vector<T> pixvals;
466 assert( &imageout != &imagein );
467 imageout.init(imagein.getWidth(),imagein.getHeight());
468  
469 for (int j=1; j<imagein.getHeight()-1; j++)
470 for (int i=1; i<imagein.getWidth()-1; i++)
471 {
472 if (!pixvals.empty()) pixvals.clear();
473 for (int y=-1; y<=1; y++)
474 for (int x=-1; x<=1; x++)
475 pixvals.push_back(imagein.pixel(i+x,j+y));
476  
477 sort(pixvals.begin(), pixvals.end());
478 imageout.pixel(i,j)=pixvals[pixvals.size()/2];
479 }
480  
481 for (int j=0; j<imagein.getHeight(); j++)
482 {
483 imageout.pixel(0,j) = imagein.pixel(0,j);
484 imageout.pixel(imagein.getWidth()-1,j)=imagein.pixel(imagein.getWidth()-1,j);
485 }
486  
487 for (int i=0; i<imagein.getWidth(); i++)
488 {
489 imageout.pixel(i,0)=imagein.pixel(i,0);
490 imageout.pixel(i,imagein.getHeight()-1)=imagein.pixel(i,imagein.getHeight()-1);
491 }
492 }
493  
494 template<typename T>
495 void despeckleMedianIt( image_ref<T> &imagein )
496 {
497 image<T> tempimage;
498  
499 despeckleMedian(imagein, tempimage);
500  
501 imagein = tempimage;
502 }
503  
504 template<typename T>
505 void chamfer(image_ref<T> &inimage)
506 {
507 int nw,n,ne,w,e,sw,s,se;
508 int fp;
509  
510 for (int j=1; j<inimage.getHeight()-1; j++)
511 for (int i=1; i<inimage.getWidth(); i++)
512 {
513 fp=inimage.pixel(i,j);
514 nw=4+inimage.pixel(i-1,j-1);
515 n=3+inimage.pixel(i,j-1);
516 w=3+inimage.pixel(i-1,j);
517 sw=4+inimage.pixel(i-1,j+1);
518  
519 fp=fp<nw?fp:nw;
520 fp=fp<n?fp:n;
521 fp=fp<w?fp:w;
522 fp=fp<sw?fp:sw;
523 inimage.pixel(i,j)=fp;
524  
525 }
526  
527 for (int j=inimage.getHeight()-2; j>0; j--)
528 for (int i=inimage.getWidth()-2; i>=0; i--)
529 {
530 fp=inimage.pixel(i,j);
531 ne=4+inimage.pixel(i+1,j-1);
532 e=3+inimage.pixel(i+1,j);
533 s=3+inimage.pixel(i,j+1);
534 se=4+inimage.pixel(i+1,j+1);
535  
536 fp=fp<ne?fp:ne;
537 fp=fp<e?fp:e;
538 fp=fp<s?fp:s;
539 fp=fp<se?fp:se;
540 inimage.pixel(i,j)=fp;
541  
542 }
543  
544 for (int i=0; i<inimage.getWidth(); i++)
545 {
546 inimage.pixel(i,0)=inimage.pixel(i,1);
547 inimage.pixel(i,inimage.getHeight()-1)=inimage.pixel(i,inimage.getHeight()-2);
548 }
549  
550 for (int j=0; j<inimage.getHeight(); j++)
551 {
552 inimage.pixel(0,j)=inimage.pixel(1,j);
553 inimage.pixel(inimage.getWidth()-1,j)=inimage.pixel(inimage.getWidth()-2,j);
554 }
555 }
556  
557  
558 template<typename T,typename TPtr>
559 void paethRotateRadian(const const_image_ref<T,TPtr> &imagein,
560 image<T> &imageout,
561 float radian)
562 {
563 assert( &imageout != &imagein );
564 imageout.init(imagein.getWidth(),imagein.getHeight());
565 imageout.fill(255);
566  
567 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
568 float s1=-tan(-radian/2.0);
569 float s2=sin(-radian);
570  
571 float xr, yr;
572 int x,y;
573 for (int j=0; j<imageout.getHeight(); j++)
574 for (int i=0; i<imageout.getWidth(); i++)
575 {
576 x=i-imageout.getWidth()/2;
577 y=j-imageout.getHeight()/2;
578 //x=i;
579 //y=j;
580  
581 xr = (float)x + s1 * (float)y;
582 yr = (float)y + s2 * xr;
583 xr = xr + s1 * yr;
584  
585 xr+=(float)imageout.getWidth()/2;
586 yr+=(float)imageout.getHeight()/2;
587 imageout.pixel(i, j)=imagein.getPixel((int)xr,(int)yr);
588 //imageout.setPixel(i,j, imagein.getPixel(i,j));
589 }
590 }
591  
592 template< typename T, typename T1Ptr, typename T2Ptr >
593 void multiLevelCorrelation( const const_image_ref< T, T1Ptr > &edgemap,
594 const const_image_ref< T, T2Ptr > &objectedgemap,
595 int levels, int &x, int &y, float &corrval,
596 bool chamferObject)
597 {
598 image<T> object[levels], image[levels];
599  
600 image[0] = edgemap;
601 // invert(image[0]);
602  
603 object[0] = objectedgemap;
604 // invert(object[0]);
605  
606 for (int l=1; l<levels; l++)
607 {
608 halfResolutionEdgePyramid(image[l-1],image[l]);
609 halfResolutionEdgePyramid(object[l-1],object[l]);
610 }
611  
612 for (int l=0; l<levels; l++)
613 {
614 image[l] = (T)255 - image[l];
615 object[l] = (T)255 - object[l];
616 //image[l].display();
617 //object[l].display();
618 chamfer(image[l]);
619 if (chamferObject) chamfer(object[l]);
620 }
621  
622 int xmin = INT_MIN, ymin = INT_MIN;
623 int xminsearch, xmaxsearch, yminsearch, ymaxsearch;
624 float corr, mincorr;
625  
626 x=0; y=0;
627 xminsearch=0;
628 yminsearch=0;
629 xmaxsearch=image[levels-1].getWidth();
630 ymaxsearch=image[levels-1].getHeight();
631  
632 for (int l=levels-1; l>=0; l--)
633 {
634 //object[l].display();
635 //image[l].display();
636  
637 mincorr=FLT_MAX;
638 for (int j=yminsearch; j<ymaxsearch-object[l].getHeight(); j++)
639 for (int i=xminsearch; i<xmaxsearch-object[l].getWidth(); i++)
640 {
641 corr=0.0;
642 int count=0;
643 for (int b=0; b<object[l].getHeight(); b++)
644 for (int a=0; a<object[l].getWidth(); a++)
645 {
646 if (!chamferObject) if (object[l].pixel(a,b)>127) continue;
647 corr+=(float)abs((int)image[l].pixel(i+a,j+b)-(int)object[l].pixel(a,b));
648 count++;
649 }
650  
651 corr/=count;
652 if (corr<mincorr)
653 {
654 xmin=i; ymin=j;
655 mincorr=corr;
656 }
657 }
658  
659 if (l!=0)
660 {
661 xminsearch=(xmin*2)-1-20;
662 yminsearch=(ymin*2)-1-20;
663  
664 xmaxsearch=(xmin*2)-1+object[l-1].getWidth()+20;
665 ymaxsearch=(ymin*2)-1+object[l-1].getHeight()+20;
666  
667 if (xminsearch<0) xminsearch=0;
668 if (yminsearch<0) yminsearch=0;
669  
670 if (xmaxsearch>image[l-1].getWidth()) xmaxsearch=image[l-1].getWidth();
671 if (ymaxsearch>image[l-1].getHeight()) ymaxsearch=image[l-1].getHeight();
672 }
673  
674  
675 }
676  
677 assert( xmin != INT_MIN && ymin != INT_MIN );
678 x=xmin;
679 y=ymin;
680 corrval=mincorr;
681  
682 }
683  
684 template< typename T, typename T1Ptr, typename T2Ptr >
685 void correlation( const const_image_ref<T> &imagein,
686 const const_image_ref<T> &templatein,
687 int &x, int &y, float &corrval)
688 {
689  
690 int xmin = INT_MIN, ymin = INT_MIN;
691 int xminsearch, xmaxsearch, yminsearch, ymaxsearch;
692 float corr, mincorr;
693  
694 x=0; y=0;
695 xminsearch=0;
696 yminsearch=0;
697 xmaxsearch=imagein.getWidth();
698 ymaxsearch=imagein.getHeight();
699  
700 mincorr=FLT_MAX;
701 for (int j=yminsearch; j<ymaxsearch-templatein.getHeight(); j++)
702 for (int i=xminsearch; i<xmaxsearch-templatein.getWidth(); i++)
703 {
704 corr=0.0;
705 for (int b=0; b<templatein.getHeight(); b++)
706 for (int a=0; a<templatein.getWidth(); a++)
707 {
708 //if (object[l].getPixel(a,b)>127) continue;
709 corr+=(float)abs((int)imagein.pixel(i+a,j+b)-(int)templatein.pixel(a,b));
710 }
711  
712 if (corr<mincorr)
713 {
714 xmin=i; ymin=j;
715 mincorr=corr;
716 }
717 }
718  
719 assert( xmin != INT_MIN && ymin != INT_MIN );
720 x=xmin;
721 y=ymin;
722 corrval=mincorr;
723  
724 }
725  
726 template<typename T, typename T1Ptr, typename T2Ptr >
727 void correlationChamfer( const const_image_ref< T, T1Ptr > &imagein,
728 const const_image_ref< T, T2Ptr > &templatein,
729 int &x, int &y, float &corrval)
730 {
731  
732 int xmin = INT_MIN, ymin = INT_MIN;
733 int xminsearch, xmaxsearch, yminsearch, ymaxsearch;
734 float corr, mincorr;
735  
736 x=0; y=0;
737 xminsearch=0;
738 yminsearch=0;
739 xmaxsearch=imagein.getWidth();
740 ymaxsearch=imagein.getHeight();
741  
742 mincorr=FLT_MAX;
743 for (int j=yminsearch; j<ymaxsearch-templatein.getHeight(); j++)
744 for (int i=xminsearch; i<xmaxsearch-templatein.getWidth(); i++)
745 {
746 corr=0.0;
747 for (int b=0; b<templatein.getHeight(); b++)
748 for (int a=0; a<templatein.getWidth(); a++)
749 {
750 if (templatein.pixel(a,b)>127) continue;
751 corr+=(float)abs((int)imagein.pixel(i+a,j+b)-(int)templatein.pixel(a,b));
752 }
753  
754 if (corr<mincorr)
755 {
756 xmin=i; ymin=j;
757 mincorr=corr;
758 }
759 }
760  
761 assert( xmin != INT_MIN && ymin != INT_MIN );
762 x=xmin;
763 y=ymin;
764 corrval=mincorr;
765  
766 }
767  
768 template<typename T>
769 void rotationChamferCorrelationEdgeMap(image<T> &edgemap, image<T>
770 &objectedgemap, int angle1, int angle2, int anglestep, int &x, int &y, int &outangle, float &corrval)
771 {
772 image<T> chamferedImage, rotatedObject, chamferedObject;
773  
774 chamferedImage = (T)255 - edgemap;
775 chamfer(chamferedImage);
776  
777  
778 objectedgemap = (T)255 - objectedgemap;
779 // Make sure that object boundaries are white:
780 drawLine(objectedgemap, 0,0,0,objectedgemap.getHeight()-1,255);
781 drawLine(objectedgemap, objectedgemap.getWidth()-1,0,objectedgemap.getWidth()-1,objectedgemap.getHeight()-1,255);
782 drawLine(objectedgemap, 0,0,objectedgemap.getWidth()-1,objectedgemap.getHeight()-1,255);
783 drawLine(objectedgemap, 0,objectedgemap.getHeight()-1,objectedgemap.getWidth()-1,objectedgemap.getHeight()-1,255);
784 //chamfer(objectedgemap);
785  
786 float cval;
787 int tx,ty;
788 corrval=FLT_MAX;
789 for (int angle=angle1; angle<=angle2; angle+=anglestep)
790 {
791 rotatedObject = objectedgemap;
792 //invert(rotatedObject);
793 paethRotateRadian(rotatedObject,chamferedObject,(float)angle*M_PI/180.0);
794 //chamfer(chamferedObject);
795 correlationChamfer(chamferedImage, chamferedObject, tx, ty, cval);
796  
797 if (cval<corrval)
798 {
799 corrval=cval;
800 x=tx;
801 y=ty;
802 outangle=angle;
803 }
804 }
805 }
806  
807 template< typename T, typename TPtr >
808 image< T > focusEnhance( const const_image_ref< T, TPtr > &im )
809 {
810 const T enhance[3][3]={ { -1, 0, -1 },
811 { 0, 7, 0 },
812 { -1, 0, -1 } };
813 const_image_ref< T > filter( &enhance[0][0], 3, 3 );
814 return normalise( correlate( im, filter ), (T)0, (T)255 );
815 }
816  
817 template< typename T, typename TPtr >
818 image< T > featureEnhance( const const_image_ref< T, TPtr > &im )
819 {
820 const T enhance[3][3]={ { 0, -1, 0 },
821 { -1, 10, -1 },
822 { 0, -1, 0 } };
823 const_image_ref< T > filter( &enhance[0][0], 3, 3 );
824 return normalise( correlate( im, filter ), (T)0, (T)255 );
825 }
826  
827 template< typename T, typename TPtr >
828 image< T > emboss( const const_image_ref< T, TPtr > &im )
829 {
830 const T diff[3][3]={ { -1, 0, 0 },
831 { 0, 0, 0 },
832 { 0, 0, 1 } };
833 const_image_ref< T > filter( &diff[0][0], 3, 3 );
834 return normalise( correlate( im, filter ), (T)0, (T)255 );
835 }
836  
837 template< typename T, typename TPtr >
838 image< T > softenHeavy( const const_image_ref< T, TPtr > &im )
839 {
840 const T median[3]={ 1, 1, 1 };
841 const_image_ref< T >
842 filterx( &median[0], 3, 1 ),
843 filtery( &median[0], 1, 3 );
844 return
845 normalise( correlate( correlate( im, filterx ), filtery ),
846 (T)0, (T)255 );
847 }
848  
849 template< typename T, typename TPtr >
850 image< T > softenMedium( const const_image_ref< T, TPtr > &im )
851 {
852 const T soften[3][3]={ { 1, 1, 1 },
853 { 1, 2, 1 },
854 { 1, 1, 1 } };
855 const_image_ref< T > filter( &soften[0][0], 3, 3 );
856 return normalise( correlate( im, filter ), (T)0, (T)255 );
857 }
858  
859 template< typename T, typename TPtr >
860 image< T > softenLight( const const_image_ref< T, TPtr > &im )
861 {
862 const T soften[3][3]={ { 6, 12, 6 },
863 { 12, 25, 12 },
864 { 6, 12, 6 } };
865 const_image_ref< T > filter( &soften[0][0], 3, 3 );
866 return normalise( correlate( im, filter ), (T)0, (T)255 );
867 }
868  
869 template< typename T, typename TPtr >
870 void oilPainting( const const_image_ref< T, TPtr > &imagein,
871 image<T> &imageout,
872 int regiondim)
873 {
874 const int fy=regiondim;
875 const int fx=regiondim;
876  
877 assert( &imageout != &imagein );
878 imageout.init(imagein.getWidth(),imagein.getHeight());
879  
880 for (int y=0; y<imagein.getHeight(); y++)
881 for (int x=0; x<imagein.getWidth(); x++)
882 {
883 if (y-fy/2<0 || y+fy/2>imagein.getHeight()-1 || x-fx/2<0 || x+fx/2>imagein.getWidth()-1)
884 {
885 imageout.pixel(x,y)=imagein.pixel(x,y);
886 continue;
887 }
888  
889 std::vector<int> pixcolor(regiondim*regiondim);
890 std::vector<int> graylevel(256);
891  
892 int count=0;
893 for (int j=-fy/2; j<=fy/2; j++)
894 for (int i=-fx/2; i<=fx/2; i++)
895 {
896 pixcolor[count]=imagein.pixel(x+i,y+j);
897 graylevel[pixcolor[count]]+=1;
898  
899 count++;
900 }
901  
902 count=0;
903 int intensity=0, maxuse=0;
904 for (int j=-fy/2; j<=fy/2; j++)
905 for (int i=-fx/2; i<=fx/2; i++)
906 {
907 if (graylevel[pixcolor[count]]>maxuse)
908 {
909 maxuse=graylevel[pixcolor[count]];
910 intensity=pixcolor[count];
911 }
912  
913 count++;
914 }
915  
916 imageout.pixel(x,y)=(T)(intensity);
917 }
918  
919  
920 }
921  
922 //some common features here with vector class ...
923 //a cleverer man could probably use the same underlying base class
924 //but for now it doesnt matter
925 template< typename T, typename T1Ptr, typename T2Ptr >
926 double dotProduct( const const_image_ref<T,T1Ptr> &im1,
927 const const_image_ref<T,T2Ptr> &im2 ) throw (mimasexception)
928 {
929 MMERROR( im1.getWidth() == im2.getWidth() &&
930 im1.getHeight() == im2.getHeight(),
931 mimasexception, ,
932 "Can't (won't) compute the dot product of two images of "
933 "differing sizes." );
934  
935 double result = 0.0;
936  
937 const int size = (signed)im1.getSize();
938 const T
939 *p = im1.rawData(),
940 *q = im2.rawData();
941  
942 for ( int i=0; i<size; i++ )
943 result += *p++ * *q++;
944  
945 return result;
946 }
947  
948 /** Fast non-maxima suppression.
949 This function provides fast non-maxima suppression.
950 The function is provided separately, so that one can use the
951 Sobel-operator as well as the Gauss-gradient together with this
952 method.
953 @author Jan Wedekind (jan@wedesoft.de)
954 @date Mon Apr 10 13:53:00 2006 */
955 template<
956 typename T1, typename T2,
957 typename TAPtr, typename TBPtr, typename TCPtr
958 >
959 image< bool > nonMaximaSuppression
960 ( const const_image_ref< T1, TAPtr > &gradientX,
961 const const_image_ref< T1, TBPtr > &gradientY,
962 const const_image_ref< T2, TCPtr > &gradientSqr,
963 T2 threshold );
964  
965 }
966  
967 #include "image_funcs.tcc"
968  
969 #endif