36#ifndef VIGRA_CORRELATION_HXX
37#define VIGRA_CORRELATION_HXX
39#include "stdimage.hxx"
40#include "inspectimage.hxx"
41#include "functorexpression.hxx"
42#include "multi_math.hxx"
43#include "multi_fft.hxx"
44#include "integral_image.hxx"
47#include "applywindowfunction.hxx"
116template <
class T1,
class S1,
123 bool clearBorders=
true)
126 "vigra::fastCrossCorrelation(): shape mismatch between input and output.");
139template<
class T1,
class S1>
140inline double integralMultiArrayWindowMean(MultiArrayView<1, T1, S1>
const & in,
144 return in[right]-in[left];
147template<
class T1,
class S1>
148inline double integralMultiArrayWindowMean(MultiArrayView<2, T1, S1>
const & in,
152 return in[lr] - in(lr[0],ul[1]) - in(ul[0],lr[1]) + in[ul];
155template<
class T1,
class S1>
156inline double integralMultiArrayWindowMean(MultiArrayView<3, T1, S1>
const & in,
160 return (in[lr] - in(lr[0],ul[1],lr[2]) - in(ul[0],lr[1],lr[2]) + in(ul[0],ul[1],lr[2]))
161 - (in(lr[0],lr[1],ul[2]) - in(lr[0],ul[1],ul[2]) - in(ul[0],lr[1],ul[2]) + in[ul]);
223template <
class T1,
class S1,
230 bool clearBorders=
true)
236 "vigra::fastNormalizedCrossCorrelation(): shape mismatch between input and output.");
238 vigra_precondition(N>0 && N<=3,
239 "vigra::fastNormalizedCrossCorrelation(): only implemented for arrays of 1, 2 or 3 dimensions.");
241 for(
unsigned int dim=0; dim<N; dim++)
243 vigra_precondition(mask.
shape()[dim] % 2 == 1,
"vigra::fastNormalizedCrossCorrelation(): Mask width has to be of odd size!");
244 vigra_precondition(in.
shape()[dim] >= mask.
shape()[dim] ,
"vigra::fastNormalizedCrossCorrelation(): Mask is larger than image!");
248 double mask_mean = 0.0,
254 double fix_denumerator = mask_size*
sqrt(mask_var);
256 if(fix_denumerator == 0)
265 norm_mask -= mask_mean;
276 sum_table2(in.
shape()+1);
282 integralMultiArray(in,sum_table.
subarray(zero_diff+1, in.
shape()+1));
283 integralMultiArraySquared(in, sum_table2.subarray(zero_diff+1, in.
shape()+1));
291 double window_mean = detail::integralMultiArrayWindowMean(sum_table, *i, *i+mask.
shape()),
292 window_squared_mean = detail::integralMultiArrayWindowMean(sum_table2, *i, *i+mask.
shape()),
293 var_denumerator =
sqrt(mask_size*window_squared_mean - window_mean*window_mean);
296 if(var_denumerator == 0)
298 out[*i+maskRadius] = 0;
302 out[*i+maskRadius] = mask_size*corr_result[*i+maskRadius]/(var_denumerator*fix_denumerator);
308template<
class MaskIterator,
class MaskAccessor>
309class CorrelationFunctor
312 CorrelationFunctor(MaskIterator m_ul, MaskIterator m_lr, MaskAccessor m_acc)
319 CorrelationFunctor(triple<MaskIterator,MaskIterator,MaskAccessor> m)
320 : m_mask_ul(m.first),
326 template <
class SrcIterator,
class SrcAccessor,
class DestIterator,
class DestAccessor>
327 void operator()(SrcIterator s, SrcAccessor s_acc, DestIterator d, DestAccessor d_acc)
const
329 using namespace vigra;
331 SrcIterator s_ul = s - windowShape()/2,
332 s_lr = s + windowShape()/2+
Diff2D(1,1);
338 SrcIterator ys = s_ul;
341 MaskIterator ym = m_mask_ul;
342 MaskIterator xm = ym;
346 for( ; ys.y != s_lr.y; ys.y++, ym.y++)
348 for(xs = ys, xm = ym; xs.x != s_lr.x; xs.x++, xm.x++)
350 res += m_mask_acc(xm)*s_acc(xs);
356 Diff2D windowShape()
const
358 return m_mask_lr - m_mask_ul;
362 MaskIterator m_mask_ul;
363 MaskIterator m_mask_lr;
364 MaskAccessor m_mask_acc;
447template <
class SrcIterator,
class SrcAccessor,
448 class MaskIterator,
class MaskAccessor,
449 class DestIterator,
class DestAccessor>
450inline void crossCorrelation(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
451 MaskIterator m_ul, MaskIterator m_lr, MaskAccessor m_acc,
452 DestIterator d_ul, DestAccessor d_acc,
453 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
455 CorrelationFunctor<MaskIterator, MaskAccessor> func(m_ul, m_lr, m_acc);
459template <
class SrcIterator,
class SrcAccessor,
460 class MaskIterator,
class MaskAccessor,
461 class DestIterator,
class DestAccessor>
462inline void crossCorrelation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
463 triple<MaskIterator, MaskIterator, MaskAccessor> mask,
464 pair<DestIterator, DestAccessor> dest,
465 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
468 mask.first, mask.second, mask.third,
469 dest.first, dest.second,
473template <
class T1,
class S1,
481 "vigra::crossCorrelation(): shape mismatch between input and output.");
487template<
class MaskIterator,
class MaskAccessor>
488class NormalizedCorrelationFunctor
492 NormalizedCorrelationFunctor(MaskIterator m_ul, MaskIterator m_lr, MaskAccessor m_acc)
502 NormalizedCorrelationFunctor(triple<MaskIterator,MaskIterator,MaskAccessor> mask)
503 : m_mask_ul(mask.first),
504 m_mask_lr(mask.second),
505 m_mask_acc(mask.third),
512 template <
class SrcIterator,
class SrcAccessor,
class DestIterator,
class DestAccessor>
513 void operator()(SrcIterator s, SrcAccessor s_acc, DestIterator d, DestAccessor d_acc)
const
515 using namespace vigra;
517 SrcIterator s_ul = s - windowShape()/2,
518 s_lr = s + windowShape()/2+
Diff2D(1,1);
530 SrcIterator ys = s_ul;
533 MaskIterator ym = m_mask_ul;
534 MaskIterator xm = ym;
536 double s1=0,s2=0, s12=0, s22=0;
538 for( ; ys.y != s_lr.y; ys.y++, ym.y++)
540 for(xs = ys, xm = ym; xs.x != s_lr.x; xs.x++, xm.x++)
544 s12 += (s1-m_avg1)*(s2-average());
545 s22 += (s2-average())*(s2-average());
554 d_acc.set(s12/
sqrt(m_s11*s22),d);
559 Diff2D windowShape()
const
561 return m_mask_lr - m_mask_ul;
569 inspectImage(srcIterRange(m_mask_ul, m_mask_lr, m_mask_acc), average);
571 MaskIterator ym = m_mask_ul;
572 MaskIterator xm = ym;
576 for( ; ym.y != m_mask_lr.y; ym.y++)
578 for(xm = ym; xm.x != m_mask_lr.x; xm.x++)
580 m_s11 += (m_mask_acc(xm)-m_avg1)*(m_mask_acc(xm)-m_avg1);
585 MaskIterator m_mask_ul;
586 MaskIterator m_mask_lr;
587 MaskAccessor m_mask_acc;
674template <
class SrcIterator,
class SrcAccessor,
675 class MaskIterator,
class MaskAccessor,
676 class DestIterator,
class DestAccessor>
678 MaskIterator m_ul, MaskIterator m_lr, MaskAccessor m_acc,
679 DestIterator d_ul, DestAccessor d_acc,
680 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
682 NormalizedCorrelationFunctor<MaskIterator, MaskAccessor> func(m_ul, m_lr, m_acc);
686template <
class SrcIterator,
class SrcAccessor,
687 class MaskIterator,
class MaskAccessor,
688 class DestIterator,
class DestAccessor>
690 triple<MaskIterator, MaskIterator, MaskAccessor> mask,
691 pair<DestIterator, DestAccessor> dest,
692 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
695 mask.first, mask.second, mask.third,
696 dest.first, dest.second,
700template <
class T1,
class S1,
708 "vigra::normalizedCrossCorrelation(): shape mismatch between input and output.");
Two dimensional difference vector.
Definition diff2d.hxx:186
Find the average pixel value in an image or ROI.
Definition inspectimage.hxx:1249
Base class for, and view to, vigra::MultiArray.
Definition multi_fwd.hxx:127
MultiArrayShape< actual_dimension >::type difference_type
Definition multi_array.hxx:739
MultiArrayView subarray(difference_type p, difference_type q) const
Definition multi_array.hxx:1530
const difference_type & shape() const
Definition multi_array.hxx:1650
void meanVariance(U *mean, U *variance) const
Definition multi_array.hxx:1782
Main MultiArray class containing the memory management.
Definition multi_fwd.hxx:131
Iterate over a virtual array where each element contains its coordinate.
Definition multi_iterator.hxx:89
Class for fixed size vectors.
Definition tinyvector.hxx:1008
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition fixedpoint.hxx:667
void fastNormalizedCrossCorrelation(...)
This function performes a fast normalized cross-correlation.
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition fixedpoint.hxx:616
void crossCorrelation(...)
This function performes a (slow) cross-correlation.
void correlateFFT(...)
Correlate an array with a kernel by means of the Fourier transform.
void inspectImage(...)
Apply read-only functor to every pixel in the image.
void fastCrossCorrelation(...)
This function performes a fast cross-correlation.
void applyWindowFunction(...)
Apply a window function to each pixels of a given image.
NumericTraits< V >::Promote prod(TinyVectorBase< V, SIZE, D1, D2 > const &l)
product of the vector's elements
Definition tinyvector.hxx:2097
void normalizedCrossCorrelation(...)
This function performes a (slow) normalized cross-correlation.
void initMultiArrayBorder(...)
Write values to the specified border values in the array.