8 #ifndef VIGRA_CORRELATION_HXX
9 #define VIGRA_CORRELATION_HXX
11 #include "stdimage.hxx"
12 #include "inspectimage.hxx"
13 #include "functorexpression.hxx"
14 #include "multi_math.hxx"
15 #include "multi_fft.hxx"
16 #include "integral_image.hxx"
19 #include "applywindowfunction.hxx"
86 template <
class T1,
class S1,
91 MultiArrayView<N, T2, S2>
const & mask,
92 MultiArrayView<N, T3, S3> out,
93 bool clearBorders=
true)
95 vigra_precondition(in.shape() == out.shape(),
96 "vigra::fastCrossCorrelation(): shape mismatch between input and output.");
101 typedef typename MultiArrayShape<N>::type Shape;
102 Shape maskRadius(
floor(mask.shape()/2));
174 template<
class T1,
class S1>
175 inline double integralMultiArrayWindowMean(MultiArrayView<1, T1, S1>
const & in,
176 typename MultiArrayView<1, T1, S1>::difference_type
const & left,
177 typename MultiArrayView<1, T1, S1>::difference_type
const & right)
179 return in[right]-in[left];
182 template<
class T1,
class S1>
183 inline double integralMultiArrayWindowMean(MultiArrayView<2, T1, S1>
const & in,
184 typename MultiArrayView<2, T1, S1>::difference_type
const & ul,
185 typename MultiArrayView<2, T1, S1>::difference_type
const & lr)
187 return in[lr] - in(lr[0],ul[1]) - in(ul[0],lr[1]) + in[ul];
190 template<
class T1,
class S1>
191 inline double integralMultiArrayWindowMean(MultiArrayView<3, T1, S1>
const & in,
192 typename MultiArrayView<3, T1, S1>::difference_type
const & ul,
193 typename MultiArrayView<3, T1, S1>::difference_type
const & lr)
195 return (in[lr] - in(lr[0],ul[1],lr[2]) - in(ul[0],lr[1],lr[2]) + in(ul[0],ul[1],lr[2]))
196 - (in(lr[0],lr[1],ul[2]) - in(lr[0],ul[1],ul[2]) - in(ul[0],lr[1],ul[2]) + in[ul]);
200 template <
class T1,
class S1,
204 inline void fastNormalizedCrossCorrelation(MultiArrayView<N, T1, S1>
const & in,
205 MultiArrayView<N, T2, S2>
const & mask,
206 MultiArrayView<N, T3, S3> out,
207 bool clearBorders=
true)
209 using namespace vigra::multi_math;
210 typedef typename MultiArrayShape<N>::type Shape;
212 vigra_precondition(in.shape() == out.shape(),
213 "vigra::fastNormalizedCrossCorrelation(): shape mismatch between input and output.");
215 vigra_precondition(N>0 && N<=3,
216 "vigra::fastNormalizedCrossCorrelation(): only implemented for arrays of 1, 2 or 3 dimensions.");
218 for(
unsigned int dim=0; dim<N; dim++)
220 vigra_precondition(mask.shape()[dim] % 2 == 1,
"vigra::fastNormalizedCrossCorrelation(): Mask width has to be of odd size!");
221 vigra_precondition(in.shape()[dim] >= mask.shape()[dim] ,
"vigra::fastNormalizedCrossCorrelation(): Mask is larger than image!");
225 double mask_mean = 0.0,
227 mask_size =
prod(mask.shape());
228 mask.meanVariance(&mask_mean, &mask_var);
231 double fix_denumerator = mask_size*
sqrt(mask_var);
233 if(fix_denumerator == 0)
240 MultiArray<N, double> norm_mask(mask.shape());
242 norm_mask -= mask_mean;
245 MultiArray<N, double> corr_result(in.shape());
252 MultiArray<N, double> sum_table(in.shape()+1),
253 sum_table2(in.shape()+1);
255 typename MultiArray<N, double>::difference_type zero_diff;
259 integralMultiArray(in,sum_table.subarray(zero_diff+1, in.shape()+1));
260 integralMultiArraySquared(in, sum_table2.subarray(zero_diff+1, in.shape()+1));
262 MultiCoordinateIterator<N> i(in.shape()-mask.shape()+1), end = i.getEndIterator();
264 Shape maskRadius(
floor(mask.shape()/2));
268 double window_mean = detail::integralMultiArrayWindowMean(sum_table, *i, *i+mask.shape()),
269 window_squared_mean = detail::integralMultiArrayWindowMean(sum_table2, *i, *i+mask.shape()),
270 var_denumerator =
sqrt(mask_size*window_squared_mean - window_mean*window_mean);
273 if(var_denumerator == 0)
275 out[*i+maskRadius] = 0;
279 out[*i+maskRadius] = mask_size*corr_result[*i+maskRadius]/(var_denumerator*fix_denumerator);
368 template<
class MaskIterator,
class MaskAccessor>
380 : m_mask_ul(m.first),
386 template <
class SrcIterator,
class SrcAccessor,
class DestIterator,
class DestAccessor>
387 void operator()(SrcIterator s, SrcAccessor s_acc, DestIterator d, DestAccessor d_acc)
const
389 using namespace vigra;
391 SrcIterator s_ul = s - windowShape()/2,
392 s_lr = s + windowShape()/2+
Diff2D(1,1);
398 SrcIterator ys = s_ul;
401 MaskIterator ym = m_mask_ul;
402 MaskIterator xm = ym;
406 for( ; ys.y != s_lr.y; ys.y++, ym.y++)
408 for(xs = ys, xm = ym; xs.x != s_lr.x; xs.x++, xm.x++)
410 res += m_mask_acc(xm)*s_acc(xs);
416 Diff2D windowShape()
const
418 return m_mask_lr - m_mask_ul;
422 MaskIterator m_mask_ul;
423 MaskIterator m_mask_lr;
424 MaskAccessor m_mask_acc;
428 template <
class SrcIterator,
class SrcAccessor,
429 class MaskIterator,
class MaskAccessor,
430 class DestIterator,
class DestAccessor>
431 inline void crossCorrelation(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
432 MaskIterator m_ul, MaskIterator m_lr, MaskAccessor m_acc,
433 DestIterator d_ul, DestAccessor d_acc,
434 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
440 template <
class SrcIterator,
class SrcAccessor,
441 class MaskIterator,
class MaskAccessor,
442 class DestIterator,
class DestAccessor>
443 inline void crossCorrelation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
444 triple<MaskIterator, MaskIterator, MaskAccessor> mask,
445 pair<DestIterator, DestAccessor> dest,
446 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
448 crossCorrelation(src.first, src.second, src.third,
449 mask.first, mask.second, mask.third,
450 dest.first, dest.second,
454 template <
class T1,
class S1,
462 "vigra::crossCorrelation(): shape mismatch between input and output.");
463 crossCorrelation(srcImageRange(in),
554 template<
class MaskIterator,
class MaskAccessor>
570 : m_mask_ul(mask.first),
571 m_mask_lr(mask.second),
572 m_mask_acc(mask.third),
579 template <
class SrcIterator,
class SrcAccessor,
class DestIterator,
class DestAccessor>
580 void operator()(SrcIterator s, SrcAccessor s_acc, DestIterator d, DestAccessor d_acc)
const
582 using namespace vigra;
584 SrcIterator s_ul = s - windowShape()/2,
585 s_lr = s + windowShape()/2+
Diff2D(1,1);
597 SrcIterator ys = s_ul;
600 MaskIterator ym = m_mask_ul;
601 MaskIterator xm = ym;
603 double s1=0,s2=0, s12=0, s22=0;
605 for( ; ys.y != s_lr.y; ys.y++, ym.y++)
607 for(xs = ys, xm = ym; xs.x != s_lr.x; xs.x++, xm.x++)
611 s12 += (s1-m_avg1)*(s2-average());
612 s22 += (s2-average())*(s2-average());
621 d_acc.set(s12/
sqrt(m_s11*s22),d);
626 Diff2D windowShape()
const
628 return m_mask_lr - m_mask_ul;
636 inspectImage(srcIterRange(m_mask_ul, m_mask_lr, m_mask_acc), average);
638 MaskIterator ym = m_mask_ul;
639 MaskIterator xm = ym;
643 for( ; ym.y != m_mask_lr.y; ym.y++)
645 for(xm = ym; xm.x != m_mask_lr.x; xm.x++)
647 m_s11 += (m_mask_acc(xm)-m_avg1)*(m_mask_acc(xm)-m_avg1);
652 MaskIterator m_mask_ul;
653 MaskIterator m_mask_lr;
654 MaskAccessor m_mask_acc;
661 template <
class SrcIterator,
class SrcAccessor,
662 class MaskIterator,
class MaskAccessor,
663 class DestIterator,
class DestAccessor>
664 inline void normalizedCrossCorrelation(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
665 MaskIterator m_ul, MaskIterator m_lr, MaskAccessor m_acc,
666 DestIterator d_ul, DestAccessor d_acc,
667 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
673 template <
class SrcIterator,
class SrcAccessor,
674 class MaskIterator,
class MaskAccessor,
675 class DestIterator,
class DestAccessor>
676 inline void normalizedCrossCorrelation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
677 triple<MaskIterator, MaskIterator, MaskAccessor> mask,
678 pair<DestIterator, DestAccessor> dest,
679 BorderTreatmentMode border = BORDER_TREATMENT_AVOID)
681 normalizedCrossCorrelation(src.first, src.second, src.third,
682 mask.first, mask.second, mask.third,
683 dest.first, dest.second,
687 template <
class T1,
class S1,
695 "vigra::normalizedCrossCorrelation(): shape mismatch between input and output.");
696 normalizedCrossCorrelation(srcImageRange(in),
705 #endif //VIGRA_CORRELATION_HXX