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"
116 template <
class T1,
class S1,
121 MultiArrayView<N, T2, S2>
const & mask,
122 MultiArrayView<N, T3, S3> out,
123 bool clearBorders=
true)
125 vigra_precondition(in.shape() == out.shape(),
126 "vigra::fastCrossCorrelation(): shape mismatch between input and output.");
131 typedef typename MultiArrayShape<N>::type Shape;
132 Shape maskRadius(
floor(mask.shape()/2));
139 template<
class T1,
class S1>
140 inline double integralMultiArrayWindowMean(MultiArrayView<1, T1, S1>
const & in,
141 typename MultiArrayView<1, T1, S1>::difference_type
const & left,
142 typename MultiArrayView<1, T1, S1>::difference_type
const & right)
144 return in[right]-in[left];
147 template<
class T1,
class S1>
148 inline double integralMultiArrayWindowMean(MultiArrayView<2, T1, S1>
const & in,
149 typename MultiArrayView<2, T1, S1>::difference_type
const & ul,
150 typename MultiArrayView<2, T1, S1>::difference_type
const & lr)
152 return in[lr] - in(lr[0],ul[1]) - in(ul[0],lr[1]) + in[ul];
155 template<
class T1,
class S1>
156 inline double integralMultiArrayWindowMean(MultiArrayView<3, T1, S1>
const & in,
157 typename MultiArrayView<3, T1, S1>::difference_type
const & ul,
158 typename MultiArrayView<3, T1, S1>::difference_type
const & lr)
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]);
223 template <
class T1,
class S1,
228 MultiArrayView<N, T2, S2>
const & mask,
229 MultiArrayView<N, T3, S3> out,
230 bool clearBorders=
true)
232 using namespace vigra::multi_math;
233 typedef typename MultiArrayShape<N>::type Shape;
235 vigra_precondition(in.shape() == out.shape(),
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,
250 mask_size =
prod(mask.shape());
251 mask.meanVariance(&mask_mean, &mask_var);
254 double fix_denumerator = mask_size*
sqrt(mask_var);
256 if(fix_denumerator == 0)
263 MultiArray<N, double> norm_mask(mask.shape());
265 norm_mask -= mask_mean;
268 MultiArray<N, double> corr_result(in.shape());
275 MultiArray<N, double> sum_table(in.shape()+1),
276 sum_table2(in.shape()+1);
278 typename MultiArray<N, double>::difference_type zero_diff;
282 integralMultiArray(in,sum_table.subarray(zero_diff+1, in.shape()+1));
283 integralMultiArraySquared(in, sum_table2.subarray(zero_diff+1, in.shape()+1));
285 MultiCoordinateIterator<N> i(in.shape()-mask.shape()+1), end = i.getEndIterator();
287 Shape maskRadius(
floor(mask.shape()/2));
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);
308 template<
class MaskIterator,
class MaskAccessor>
309 class 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;
447 template <
class SrcIterator,
class SrcAccessor,
448 class MaskIterator,
class MaskAccessor,
449 class DestIterator,
class DestAccessor>
450 inline 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);
459 template <
class SrcIterator,
class SrcAccessor,
460 class MaskIterator,
class MaskAccessor,
461 class DestIterator,
class DestAccessor>
462 inline 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,
473 template <
class T1,
class S1,
481 "vigra::crossCorrelation(): shape mismatch between input and output.");
487 template<
class MaskIterator,
class MaskAccessor>
488 class 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;
674 template <
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);
686 template <
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,
700 template <
class T1,
class S1,
708 "vigra::normalizedCrossCorrelation(): shape mismatch between input and output.");
718 #endif //VIGRA_CORRELATION_HXX
void fastNormalizedCrossCorrelation(...)
This function performes a fast normalized cross-correlation.
void fastCrossCorrelation(...)
This function performes a fast cross-correlation.
void initMultiArrayBorder(...)
Write values to the specified border values in the array.
const difference_type & shape() const
Definition: multi_array.hxx:1648
Two dimensional difference vector.
Definition: diff2d.hxx:185
Find the average pixel value in an image or ROI.
Definition: inspectimage.hxx:1248
NumericTraits< V >::Promote prod(TinyVectorBase< V, SIZE, D1, D2 > const &l)
product of the vector's elements
Definition: tinyvector.hxx:2097
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void crossCorrelation(...)
This function performes a (slow) cross-correlation.
void correlateFFT(...)
Correlate an array with a kernel by means of the Fourier transform.
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:704
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
void normalizedCrossCorrelation(...)
This function performes a (slow) normalized cross-correlation.
void inspectImage(...)
Apply read-only functor to every pixel in the image.
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
void applyWindowFunction(...)
Apply a window function to each pixels of a given image.