Correlation.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00024 #ifndef VIGRA_EXT_CORRELATION_H
00025 #define VIGRA_EXT_CORRELATION_H
00026 
00027 #include <hugin_shared.h>
00028 #include <vigra/stdimage.hxx>
00029 #include <vigra/inspectimage.hxx>
00030 #include <vigra/copyimage.hxx>
00031 #include <vigra/resizeimage.hxx>
00032 #include <vigra/transformimage.hxx>
00033 
00034 #include "hugin_utils/utils.h"
00035 #include "hugin_math/hugin_math.h"
00036 #include "vigra_ext/FitPolynom.h"
00037 #include "vigra_ext/utils.h"
00038 
00039 // hmm.. not really great.. should be part of vigra_ext
00040 #include "vigra_ext/ImageTransforms.h"
00041 #include "hugin_config.h"
00042 #ifdef HAVE_FFTW
00043 #include <vigra/fftw3.hxx>
00044 #include <vigra/functorexpression.hxx>
00045 #include <hugin_utils/openmp_lock.h>
00046 #include <vector>
00047 #else
00048 #define VIGRA_EXT_USE_FAST_CORR
00049 #endif
00050 
00051 //#define DEBUG_WRITE_FILES
00052 
00053 namespace vigra_ext{
00054 
00056 struct CorrelationResult
00057 {
00058     CorrelationResult()
00059         : maxi(-1), maxpos(0,0), curv(0,0), maxAngle(0)
00060         { }
00061     // value at correlation peak.
00062     double maxi;
00063     // position of maximum
00064     hugin_utils::FDiff2D maxpos;
00065     // position of correlated point, used only when using projection aware routine
00066     hugin_utils::FDiff2D corrPos;
00067     // curvature of the correlation peak
00068     hugin_utils::FDiff2D curv;
00069     double maxAngle;
00070 };
00071 
00072 #ifdef HAVE_FFTW
00073 
00075 static hugin_omp::Lock fftLock;
00076 
00077 // multiplication with conjugate number in Fourier space
00078 template <class Value>
00079 Value multiplyConjugate(Value const& v1, Value const& v2)
00080 {
00081     return v1 * vigra::conj(v2);
00082 };
00083 
00092 template <class SrcImage, class DestImage, class KernelImage>
00093 CorrelationResult correlateImageFastFFT(SrcImage & src, DestImage & dest, KernelImage & kernel, vigra::Diff2D kul, vigra::Diff2D klr)
00094 {
00095     vigra_precondition(kul.x <= 0 && kul.y <= 0,
00096         "correlateImageFastFFT(): coordinates of kernel's upper left must be <= 0.");
00097     vigra_precondition(klr.x >= 0 && klr.y >= 0,
00098         "correlateImageFastFFT(): coordinates of kernel's lower right must be >= 0.");
00099 
00100     // calculate width and height of the image
00101     const int kw = kernel.width();
00102     const int kh = kernel.height();
00103     const int sw = src.width();
00104     const int sh = src.height();
00105     vigra_precondition(sw >= kw && sh >= kh, "correlateImageFFT(): kernel larger than image.");
00106 
00107     CorrelationResult res;
00108     if (sw < kw || sh < kh)
00109     {
00110         return res;
00111     };
00112 
00113     // calculate mean and variance of kernel/template
00114     vigra::FindAverageAndVariance<typename KernelImage::PixelType> kMean;
00115     vigra::inspectImage(srcImageRange(kernel), kMean);
00116     // uniform patch, skip calculation
00117     if (kMean.variance(false) == 0)
00118     {
00119         return res;
00120     };
00121     // subtract mean from kernel/template
00122     vigra::transformImage(srcImageRange(kernel), destImage(kernel), vigra::functor::Arg1() - vigra::functor::Param(kMean.average()));
00123     // FFT: we are reusing the fftw_plan 
00124     vigra::FFTWComplexImage spatial(sw, sh);
00125     vigra::FFTWComplexImage fourier(sw, sh);
00126     vigra::FFTWComplexImage FKernel(sw, sh);
00127     // copy kernel to complex data structure
00128     vigra::copyImage(srcImageRange(kernel), destImage(spatial, vigra::FFTWWriteRealAccessor<typename KernelImage::PixelType>()));
00129     // now do FFT of kernel
00130     fftw_plan fftwplan;
00131     {
00132         // creation of plan is not thread safe, so use lock
00133         hugin_omp::ScopedLock sl(fftLock);
00134         fftwplan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(), FFTW_FORWARD, FFTW_ESTIMATE);
00135     };
00136     fftw_execute(fftwplan);
00137     vigra::copyImage(srcImageRange(fourier), destImage(FKernel));
00138     // now do FFT of search image, reuse fftw_plan
00139     vigra::copyImage(srcImageRange(src), destImage(spatial, vigra::FFTWWriteRealAccessor<typename KernelImage::PixelType>()));
00140     fftw_execute(fftwplan);
00141     // give now not anymore needed memory free
00142     spatial.resize(0, 0);
00143 
00144     // multiply SrcImage with conjugated kernel in frequency domain
00145     vigra::combineTwoImages(srcImageRange(fourier), srcImage(FKernel), destImage(fourier), &multiplyConjugate<vigra::FFTWComplex<typename KernelImage::PixelType> >);
00146     // FFT back into spatial domain (inplace)
00147     fftw_plan backwardPlan;
00148     {
00149         // creation of plan is not thread safe, so use lock
00150         hugin_omp::ScopedLock sl(fftLock);
00151         backwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)fourier.begin(), (fftw_complex *)fourier.begin(), FFTW_BACKWARD, FFTW_ESTIMATE);
00152     };
00153     fftw_execute(backwardPlan);
00154     {
00155         // delayed destroy to use only one lock
00156         hugin_omp::ScopedLock sl(fftLock);
00157         fftw_destroy_plan(fftwplan);
00158         fftw_destroy_plan(backwardPlan);
00159     };
00160 
00161     // calculate look up sum tables
00162     // use double instead of float!, otherwise there can be truncation errors
00163     vigra::DImage s(sw, sh);
00164     vigra::DImage s2(sw, sh);
00165     double val = src(0, 0);
00166     s(0, 0) = val;
00167     s2(0, 0) = val * val;
00168     // special treatment for first line
00169     for (int x = 1; x < sw; ++x)
00170     {
00171         val = src(x, 0);
00172         s(x, 0) = val + s(x - 1, 0);
00173         s2(x, 0) = val * val + s2(x - 1, 0);
00174     }
00175     // special treatment for first column
00176     for (int y = 1; y < sh; ++y)
00177     {
00178         val = src(0, y);
00179         s(0, y) = val + s(0, y - 1);
00180         s2(0, y) = val * val + s2(0, y - 1);
00181     }
00182     // final summation
00183     for (int y = 1; y < sh; ++y)
00184     {
00185         for (int x = 1; x < sw; ++x)
00186         {
00187             val = src(x, y);
00188             s(x, y) = val + s(x - 1, y) + s(x, y - 1) - s(x - 1, y - 1);
00189             s2(x, y) = val * val + s2(x - 1, y) + s2(x, y - 1) - s2(x - 1, y - 1);
00190         };
00191     };
00192     const int yend = sh - klr.y + kul.y;
00193     const int xend = sw - klr.x + kul.x;
00194     // calculate constant part
00195     const double normFactor = 1.0 / (sw * sh * sqrt(kMean.variance(false)));
00196     for (int yr = 0; yr < yend; ++yr)
00197     {
00198         for (int xr = 0; xr < xend; ++xr)
00199         {
00200             double value = fourier(xr, yr).re() * normFactor;
00201             // do final summation using the lookup tables
00202             double sumF = s(xr + kw - 1, yr + kh - 1);
00203             double sumF2 = s2(xr + kw - 1, yr + kh - 1);
00204             if (xr > 0)
00205             {
00206                 sumF -= s(xr - 1, yr + kh - 1);
00207                 sumF2 -= s2(xr - 1, yr + kh - 1);
00208             };
00209             if (yr > 0)
00210             {
00211                 sumF -= s(xr + kw - 1, yr - 1);
00212                 sumF2 -= s2(xr + kw - 1, yr - 1);
00213             };
00214             if (xr > 0 && yr > 0)
00215             {
00216                 sumF += s(xr - 1, yr - 1);
00217                 sumF2 += s2(xr - 1, yr - 1);
00218             };
00219 
00220             double den = sqrt((kw * kh * sumF2 - sumF * sumF));
00221             // prevent division through zero
00222             if (den != 0)
00223             {
00224                 value /= den;
00225                 if (value > res.maxi)
00226                 {
00227                     res.maxi = value;
00228                     res.maxpos.x = xr - kul.x;
00229                     res.maxpos.y = yr - kul.y;
00230                 }
00231                 dest(xr - kul.x, yr - kul.y) = vigra::NumericTraits<typename DestImage::value_type>::fromRealPromote(value);
00232             };
00233         };
00234     };
00235     return res;
00236 };
00237 
00238 #endif
00239 
00252 template <class SrcImage, class DestImage, class KernelImage>
00253 CorrelationResult correlateImageFast(SrcImage & src,
00254                                      DestImage & dest,
00255                                      KernelImage & kernel,
00256                                      vigra::Diff2D kul, vigra::Diff2D klr,
00257                                      double threshold = 0.7 )
00258 {
00259     vigra_precondition(kul.x <= 0 && kul.y <= 0,
00260                  "convolveImage(): coordinates of "
00261                  "kernel's upper left must be <= 0.");
00262     vigra_precondition(klr.x >= 0 && klr.y >= 0,
00263                  "convolveImage(): coordinates of "
00264                  "kernel's lower right must be >= 0.");
00265 
00266     // use traits to determine SumType as to prevent possible overflow
00267     typedef typename
00268         vigra::NumericTraits<typename SrcImage::value_type>::RealPromote SumType;
00269     typedef typename
00270         vigra::NumericTraits<typename KernelImage::value_type>::RealPromote KSumType;
00271     typedef
00272         vigra::NumericTraits<typename DestImage::value_type> DestTraits;
00273 
00274     // calculate width and height of the image
00275     int w = src.width();
00276     int h = src.height();
00277     int wk = kernel.width();
00278     int hk = kernel.height();
00279 
00280 //    DEBUG_DEBUG("correlate Image srcSize " << (slr - sul).x << "," << (slr - sul).y
00281 //                << " tmpl size: " << wk << "," << hk);
00282 
00283     vigra_precondition(w >= wk && h >= hk,
00284                  "convolveImage(): kernel larger than image.");
00285 
00286     int ystart = -kul.y;
00287     int yend   = h-klr.y;
00288     int xstart = -kul.x;
00289     int xend   = w-klr.x;
00290 
00291     // mean of kernel
00292     KSumType kmean=0;
00293     for(int y=0; y < hk; y++) {
00294         for(int x=0; x < wk; x++) {
00295             kmean += kernel(x,y);
00296         }
00297     }
00298     kmean = kmean / (hk*wk);
00299 
00300     CorrelationResult res;
00301 
00302 //    DEBUG_DEBUG("size: " << w << "," <<  h << " ystart: " << ystart <<", yend: " << yend);
00303     int unifwarned=0;
00304     for(int yr=ystart; yr < yend; ++yr)
00305     {
00306         for(int xr=xstart; xr < xend; ++xr)
00307         {
00308             if (dest(xr,yr) < threshold) {
00309                 continue;
00310             }
00311             // init the sum
00312             SumType numerator = 0;
00313             SumType div1 = 0;
00314             SumType div2 = 0;
00315             SumType spixel = 0;
00316             KSumType kpixel = 0;
00317 
00318             // mean of image patch
00319             KSumType mean=0;
00320                         int ym;
00321             for(ym=yr+kul.y; ym <= yr+klr.y; ym++) {
00322                 for(int xm=xr+kul.x; xm <= xr+klr.x; xm++) {
00323                     mean += src(xm,ym);
00324                 }
00325             }
00326             mean = mean / (hk*wk);
00327 
00328             // perform correlation (inner loop)
00329             ym=yr+kul.y;
00330             int yk;
00331             for(yk=0; yk < hk; yk++) {
00332                 int xm=xr+kul.x;
00333                 int xk;
00334                 for(xk=0; xk < wk; xk++) {
00335                     spixel = src(xm,ym) - mean;
00336                     kpixel = kernel(xk,yk) - kmean;
00337                     numerator += kpixel * spixel;
00338                     div1 += kpixel * kpixel;
00339                     div2 += spixel * spixel;
00340                     xm++;
00341                 }
00342                 ym++;
00343             }
00344             if (div1*div2 == 0) {
00345                 // This happens when one of the patches is perfectly uniform
00346                 numerator = 0;  // Set correlation to zero since this is uninteresting
00347                 if (!unifwarned) {
00348                     DEBUG_DEBUG("Uniform patch(es) during correlation computation");
00349                     unifwarned=1;
00350                 }
00351             } else
00352                 numerator = (numerator/sqrt(div1 * div2));
00353             
00354             if (numerator > res.maxi) {
00355                 res.maxi = numerator;
00356                 res.maxpos.x = xr;
00357                 res.maxpos.y = yr;
00358             }
00359             dest(xr,yr) = DestTraits::fromRealPromote(numerator);
00360         }
00361     }
00362     return res;
00363 }
00364 
00372 template <class Iterator, class Accessor>
00373 CorrelationResult subpixelMaxima(vigra::triple<Iterator, Iterator, Accessor> img,
00374                                  vigra::Diff2D max)
00375 {
00376     const int interpWidth = 1;
00377     CorrelationResult res;
00378     vigra_precondition(max.x > interpWidth && max.y > interpWidth,
00379                  "subpixelMaxima(): coordinates of "
00380                  "maxima must be > interpWidth, interpWidth.");
00381     vigra::Diff2D sz = img.second - img.first;
00382     vigra_precondition(sz.x - max.x >= interpWidth && sz.y - max.y >= interpWidth,
00383                  "subpixelMaxima(): coordinates of "
00384                  "maxima must interpWidth pixels from the border.");
00385     typedef typename Accessor::value_type T;
00386 
00387     T x[2*interpWidth+1];
00388     T zx[2*interpWidth+1];
00389     T zy[2*interpWidth+1];
00390 
00391 #ifdef DEBUG_CORRELATION
00392     exportImage(img,vigra::ImageExportInfo("test.tif"));
00393 #endif
00394     
00395     Accessor acc = img.third;
00396     Iterator begin=img.first;
00397     for (int i=-interpWidth; i<=interpWidth; i++) {
00398         // collect first row
00399         x[i+interpWidth] = i;
00400         zx[i+interpWidth] = acc(begin, max + vigra::Diff2D(i,0));
00401         zy[i+interpWidth] = acc(begin, max + vigra::Diff2D(0,i));
00402     }
00403 
00404     double a,b,c;
00405     FitPolynom(x, x + 2*interpWidth+1, zx, a,b,c);
00406     if (std::isnan(a) || std::isnan(b) || std::isnan(c)) {
00407 #ifdef DEBUG_CORRELATION
00408         exportImage(img,vigra::ImageExportInfo("test.tif"));
00409 #endif
00410         DEBUG_DEBUG("Bad polynomial fit results");
00411         res.maxpos.x=max.x;
00412         res.maxpos.y=max.y;
00413         return res;
00414     }
00415 
00416     // calculate extrema of x position by setting
00417     // the 1st derivate to zero
00418     // 2*c*x + b = 0
00419     if (c==0)
00420         res.maxpos.x=0;
00421     else
00422         res.maxpos.x = -b/(2*c);
00423 
00424     res.curv.x = 2*c;
00425 
00426     // calculate result at maxima
00427     double maxx = c*res.maxpos.x*res.maxpos.x + b*res.maxpos.x + a;
00428 
00429     FitPolynom(x, x + 2*interpWidth+1, zy, a,b,c);
00430     if (std::isnan(a) || std::isnan(b) || std::isnan(c))
00431     {
00432         DEBUG_DEBUG("Bad polynomial fit results");
00433         res.maxpos.x = max.x;
00434         res.maxpos.y = max.y;
00435         return res;
00436     }
00437     // calculate extrema of y position
00438     if (c==0)
00439         res.maxpos.y=0;
00440     else
00441         res.maxpos.y = -b/(2*c);
00442 
00443     res.curv.y = 2*c;
00444     double maxy = c*res.maxpos.y*res.maxpos.y + b*res.maxpos.y + a;
00445 
00446     // use mean of both maxima as new interpolation result
00447     res.maxi = (maxx + maxy) / 2;
00448     DEBUG_DEBUG("value at subpixel maxima (" << res.maxpos.x << " , "
00449                 <<res.maxpos.y << "): " << maxx << "," << maxy
00450                 << " mean:" << res.maxi << " coeff: a=" << a
00451                 << "; b=" << b << "; c=" << c);
00452 
00453     // test if the point has moved too much ( more than 1.5 pixel).
00454     // actually, I should also test the confidence of the fit.
00455     if (fabs(res.maxpos.x) > 1 || fabs(res.maxpos.y) > 1) {
00456         DEBUG_NOTICE("subpixel Maxima has moved to much, ignoring");
00457         res.maxpos.x = max.x;
00458         res.maxpos.y = max.y;
00459     } else {
00460         res.maxpos.x = res.maxpos.x + max.x;
00461         res.maxpos.y = res.maxpos.y + max.y;
00462     }
00463     return res;
00464 }
00465 
00468 class RotateTransform
00469 {
00470 public:
00471     RotateTransform(double phi, hugin_utils::FDiff2D origin, hugin_utils::FDiff2D transl)
00472         : m_phi(phi), m_origin(origin), m_transl(transl)
00473         { }
00474 
00475     bool transformImgCoord(double &destx, double &desty, double srcx, double srcy) const
00476     {
00477         srcx -= m_origin.x;
00478         srcy -= m_origin.y;
00479 
00480         destx = srcx * cos(m_phi) + srcy * sin(m_phi);
00481         desty = srcx * -sin(m_phi) + srcy * cos(m_phi);
00482 
00483         destx += m_transl.x;
00484         desty += m_transl.y;
00485         return true;
00486     }
00487 
00488     double m_phi;
00489     hugin_utils::FDiff2D m_origin;
00490     hugin_utils::FDiff2D m_transl;
00491 };
00492 
00503 template <class IMAGET, class ACCESSORT, class IMAGES, class ACCESSORS>
00504 CorrelationResult PointFineTune(const IMAGET & templImg, ACCESSORT access_t,
00505                                 vigra::Diff2D templPos,
00506                                 int templSize,
00507                                 const IMAGES & searchImg, ACCESSORS access_s,
00508                                 vigra::Diff2D searchPos,
00509                                 int sWidth)
00510 {
00511 //    DEBUG_TRACE("templPos: " vigra::<< templPos << " searchPos: " vigra::<< searchPos);
00512 
00513     // extract patch from template
00514 
00515     int templWidth = templSize/2;
00516     vigra::Diff2D tmplUL(templPos.x - templWidth, templPos.y-templWidth);
00517     // lower right iterators "are past the end"
00518     vigra::Diff2D tmplLR(templPos.x + templWidth + 1, templPos.y + templWidth + 1);
00519     // clip corners to ensure the template is inside the image.
00520     vigra::Diff2D tmplImgSize(templImg.size());
00521     tmplUL = hugin_utils::simpleClipPoint(tmplUL, vigra::Diff2D(0,0), tmplImgSize);
00522     tmplLR = hugin_utils::simpleClipPoint(tmplLR, vigra::Diff2D(0,0), tmplImgSize);
00523     vigra::Diff2D tmplSize = tmplLR - tmplUL;
00524     DEBUG_DEBUG("template position: " << templPos << "  tmplUL: " << tmplUL
00525                     << "  templLR:" << tmplLR << "  size:" << tmplSize);
00526 
00527     // extract patch from search region
00528     // make search region bigger, so that interpolation can always be done
00529     int swidth = sWidth/2 +(2+templWidth);
00530 //    DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
00531 //    Diff2D subjPoint(searchPos);
00532     // clip search window
00533     if (searchPos.x < 0) searchPos.x = 0;
00534     if (searchPos.x > (int) searchImg.width()) searchPos.x = searchImg.width()-1;
00535     if (searchPos.y < 0) searchPos.y = 0;
00536     if (searchPos.y > (int) searchImg.height()) searchPos.x = searchImg.height()-1;
00537 
00538     vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
00539     // point past the end
00540     vigra::Diff2D searchLR(searchPos.x + swidth+1, searchPos.y + swidth+1);
00541     // clip search window
00542     vigra::Diff2D srcImgSize(searchImg.size());
00543     searchUL = hugin_utils::simpleClipPoint(searchUL, vigra::Diff2D(0,0), srcImgSize);
00544     searchLR = hugin_utils::simpleClipPoint(searchLR, vigra::Diff2D(0,0), srcImgSize);
00545 //    DEBUG_DEBUG("search borders: " << searchLR.x << "," << searchLR.y);
00546 
00547     vigra::Diff2D searchSize = searchLR - searchUL;
00548     // create output image
00549 
00550     // source input
00551     vigra::FImage srcImage(searchLR-searchUL);
00552     vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
00553                                         searchImg.upperLeft() + searchLR,
00554                                         access_s),
00555                      destImage(srcImage) );
00556 
00557     vigra::FImage templateImage(tmplSize);
00558     vigra::copyImage(vigra::make_triple(templImg.upperLeft() + tmplUL,
00559                                         templImg.upperLeft() + tmplLR,
00560                                         access_t),
00561                      destImage(templateImage));
00562 #ifdef DEBUG_WRITE_FILES
00563     vigra::ImageExportInfo tmpli("hugin_templ.tif");
00564     vigra::exportImage(vigra::srcImageRange(templateImage), tmpli);
00565 
00566     vigra::ImageExportInfo srci("hugin_searchregion.tif");
00567     vigra::exportImage(vigra::srcImageRange(srcImage), srci);
00568 #endif
00569 
00570     vigra::FImage dest(searchSize);
00571     dest.init(-1);
00572     // we could use the multiresolution version as well.
00573     // but usually the region is quite small.
00574     CorrelationResult res;
00575 #ifdef HAVE_FFTW
00576     DEBUG_DEBUG("+++++ starting fast correlation");
00577     res = correlateImageFastFFT(srcImage, dest, templateImage,
00578         tmplUL - templPos, tmplLR - templPos - vigra::Diff2D(1, 1));
00579 #elif defined VIGRA_EXT_USE_FAST_CORR
00580     DEBUG_DEBUG("+++++ starting fast correlation");
00581     res = correlateImageFast(srcImage,
00582         dest,
00583         templateImage,
00584         tmplUL - templPos, tmplLR - templPos - vigra::Diff2D(1, 1),
00585         -1);
00586 #else
00587     DEBUG_DEBUG("+++++ starting normal correlation");
00588     res = correlateImage(srcImage.upperLeft(),
00589                          srcImage.lowerRight(),
00590                          srcImage.accessor(),
00591                          dest.upperLeft(),
00592                          dest.accessor(),
00593                          templateImage.upperLeft() + templPos,
00594                          templateImage.accessor(),
00595                          tmplUL, tmplLR, -1);
00596 #endif
00597     DEBUG_DEBUG("normal search finished, max:" << res.maxi
00598                 << " at " << res.maxpos);
00599     // do a subpixel maxima estimation
00600     // check if the max is inside the pixel boundaries,
00601     // and there are enought correlation values for the subpixel
00602     // estimation, (2 + templWidth)
00603     if (res.maxpos.x > 2 + templWidth && res.maxpos.x < 2*swidth+1-2-templWidth
00604         && res.maxpos.y > 2+templWidth && res.maxpos.y < 2*swidth+1-2-templWidth)
00605     {
00606         // subpixel estimation
00607         res = subpixelMaxima(vigra::srcImageRange(dest), res.maxpos.toDiff2D());
00608         DEBUG_DEBUG("subpixel position: max:" << res.maxi
00609                     << " at " << res.maxpos);
00610     } else {
00611         // not enough values for subpixel estimation.
00612         DEBUG_DEBUG("subpixel estimation not done, maxima too close to border");
00613     }
00614 
00615     res.maxpos = res.maxpos + searchUL;
00616     return res;
00617 }
00618 
00629 #ifndef HAVE_FFTW
00630 template <class IMAGET, class IMAGES>
00631 CorrelationResult PointFineTuneRotSearch(const IMAGET & templImg,
00632                                          vigra::Diff2D templPos,
00633                                          int templSize,
00634                                          const IMAGES & searchImg,
00635                                          vigra::Diff2D searchPos,
00636                                          int sWidth,
00637                                          double startAngle,
00638                                          double stopAngle,
00639                                          int angleSteps)
00640 {
00641     DEBUG_TRACE("templPos: " << templPos << " searchPos: " << searchPos);
00642 
00643     // extract patch from template
00644 
00645     // make template size user configurable as well?
00646     int templWidth = templSize/2;
00647     vigra::Diff2D tmplUL(-templWidth, -templWidth);
00648     vigra::Diff2D tmplLR(templWidth, templWidth);
00649     // clip template
00650     if (tmplUL.x + templPos.x < 0) tmplUL.x = -templPos.x;
00651     if (tmplUL.y + templPos.y < 0) tmplUL.y = -templPos.y;
00652     if (tmplLR.x + templPos.x> templImg.width())
00653         tmplLR.x = templImg.width() - templPos.x;
00654     if (tmplLR.y + templPos.y > templImg.height())
00655         tmplLR.y = templImg.height() - templPos.y;
00656     vigra::Diff2D tmplSize = tmplLR - tmplUL + vigra::Diff2D(1,1);
00657 
00658     // extract patch from search region
00659     // make search region bigger, so that interpolation can always be done
00660     int swidth = sWidth/2 +(2+templWidth);
00661     DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
00662 //    Diff2D subjPoint(searchPos);
00663     // clip search window
00664     if (searchPos.x < 0) searchPos.x = 0;
00665     if (searchPos.x > (int) searchImg.width()) searchPos.x = searchImg.width()-1;
00666     if (searchPos.y < 0) searchPos.y = 0;
00667     if (searchPos.y > (int) searchImg.height()) searchPos.x = searchImg.height()-1;
00668 
00669     vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
00670     vigra::Diff2D searchLR(searchPos.x + swidth, searchPos.y + swidth);
00671     // clip search window
00672     if (searchUL.x < 0) searchUL.x = 0;
00673     if (searchUL.x > searchImg.width()) searchUL.x = searchImg.width();
00674     if (searchUL.y < 0) searchUL.y = 0;
00675     if (searchUL.y > searchImg.height()) searchUL.y = searchImg.height();
00676     if (searchLR.x > searchImg.width()) searchLR.x = searchImg.width();
00677     if (searchLR.x < 0) searchLR.x = 0;
00678     if (searchLR.y > searchImg.height()) searchLR.y = searchImg.height();
00679     if (searchLR.y < 0) searchLR.y = 0;
00680     DEBUG_DEBUG("search borders: " << searchUL << "," << searchLR << " size: " << searchLR -searchUL);
00681 
00682 
00683 #ifdef DEBUG_WRITE_FILES
00684     DEBUG_DEBUG("template area: " << templPos + tmplUL << " -> " << templPos + tmplLR);
00685     vigra::exportImage(vigra::make_triple(templImg.upperLeft() + templPos + tmplUL,
00686                                           templImg.upperLeft() + templPos + tmplLR + vigra::Diff2D(1,1),
00687                                           templImg.accessor()),
00688                        vigra::ImageExportInfo("00_original_template.png"));
00689 
00690     vigra::exportImage(make_triple(searchImg.upperLeft() + searchUL,
00691                                    searchImg.upperLeft() + searchLR,
00692                                    searchImg.accessor()),
00693                        vigra::ImageExportInfo("00_searcharea.png"));
00694 #endif
00695 
00696     // rotated template
00697     vigra::FImage rotTemplate(tmplSize);
00698     vigra::FImage alpha(tmplSize);
00699 
00700     // correlation output
00701     vigra::Diff2D searchSize = searchLR - searchUL;
00702     vigra::FImage dest(searchSize);
00703     dest.init(1);
00704     vigra::FImage bestCorrelation(searchSize);
00705 
00706     // source input
00707     vigra::FImage srcImage(searchLR-searchUL);
00708     vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
00709                                         searchImg.upperLeft() + searchLR,
00710                                         vigra::RGBToGrayAccessor<typename IMAGES::value_type>() ),
00711                      destImage(srcImage) );
00712 
00713     CorrelationResult bestRes;
00714     bestRes.maxi = -1;
00715     double bestAngle = 0;
00716 
00717     AppBase::DummyProgressDisplay dummy;
00718     // test the image at rotation angles with 30 deg. steps.
00719     double step = (stopAngle - startAngle)/(angleSteps-1);
00720     double phi=startAngle;
00721     vigra_ext::PassThroughFunctor<float> nf;
00722     for (double i=0; phi <= stopAngle; i++, phi += step) {
00723         DEBUG_DEBUG("+++++ Rotating image, phi: " << RAD_TO_DEG(phi));
00724         RotateTransform t(phi, hugin_utils::FDiff2D(templWidth, templWidth), templPos);
00725         vigra_ext::transformImage(srcImageRange(templImg, vigra::RGBToGrayAccessor<typename IMAGET::value_type>()),
00726                            destImageRange(rotTemplate),
00727                            destImage(alpha),
00728                            vigra::Diff2D(0,0),
00729                            t,
00730                            nf,
00731                            false,
00732                            vigra_ext::INTERP_CUBIC,
00733                            &dummy);
00734         DEBUG_DEBUG("----- Image rotated");
00735 
00736         // force a search in at all points.
00737         dest.init(-1);
00738         DEBUG_DEBUG("+++++ starting correlation");
00739 
00740         CorrelationResult res;
00741         // we could use the multiresolution version as well.
00742         // but usually the region is quite small.
00743 #ifdef HAVE_FFTW
00744         res = correlateImageFastFFT(srcImage, dest, rotTemplate,
00745             vigra::Diff2D(-templWidth, -templWidth), vigra::Diff2D(templWidth, templWidth));
00746 #elif defined VIGRA_EXT_USE_FAST_CORR
00747         res = correlateImageFast(srcImage,
00748                                  dest,
00749                                  rotTemplate,
00750                                  vigra::Diff2D(-templWidth, -templWidth),
00751                                  vigra::Diff2D(templWidth, templWidth),
00752                                  -1);
00753 #else
00754         res = correlateImage(srcImage.upperLeft(),
00755                              srcImage.lowerRight(),
00756                              srcImage.accessor(),
00757                              dest.upperLeft(),
00758                              dest.accessor(),
00759                              rotTemplate.upperLeft() + vigra::Diff2D(templWidth, templWidth),
00760                              rotTemplate.accessor(),
00761                              vigra::Diff2D(-templWidth, -templWidth), vigra::Diff2D(templWidth, templWidth),
00762                              -1);
00763 #endif
00764         DEBUG_DEBUG("---- correlation finished");
00765 
00766 #ifdef DEBUG_WRITE_FILES
00767         char fname[256];
00768         vigra::BImage tmpImg(rotTemplate.size());
00769         vigra::copyImage(srcImageRange(rotTemplate),destImage(tmpImg));
00770         sprintf(fname, "%3.0f_rotated_template.png", phi/M_PI*180);
00771         vigra::exportImage(srcImageRange(tmpImg), vigra::ImageExportInfo(fname));
00772 
00773         vigra::transformImage(vigra::srcImageRange(dest), vigra::destImage(dest),
00774                               vigra::linearRangeMapping(
00775                                   -1, 1,               // src range
00776                                   (unsigned char)0, (unsigned char)255) // dest range
00777             );
00778         tmpImg.resize(dest.size());
00779         vigra::copyImage(srcImageRange(dest),destImage(tmpImg));
00780         sprintf(fname, "%3.0f_corr_result.png", phi/M_PI*180);
00781         vigra::exportImage(srcImageRange(tmpImg), vigra::ImageExportInfo(fname));
00782 #endif
00783 
00784 
00785         DEBUG_DEBUG("normal search finished, max:" << res.maxi
00786                     << " at " << res.maxpos << " angle:" << phi/M_PI*180 << "");
00787 
00788         if (res.maxi > bestRes.maxi) {
00789             // remember best correlation.
00790             bestCorrelation = dest;
00791             bestRes = res;
00792             bestAngle = phi;
00793         }
00794 
00795     }
00796 
00797     DEBUG_DEBUG("rotation search finished, max:" << bestRes.maxi
00798                 << " at " << bestRes.maxpos << " angle:" << bestAngle/M_PI*180 << "");
00799 
00800     // do a subpixel maxima estimation
00801     // check if the max is inside the pixel boundaries,
00802     // and there are enought correlation values for the subpixel
00803     // estimation, (2 + templWidth)
00804     if (bestRes.maxpos.x > 1 + templWidth && bestRes.maxpos.x < 2*swidth+1-1-templWidth
00805         && bestRes.maxpos.y > 1+templWidth && bestRes.maxpos.y < 2*swidth+1-1-templWidth)
00806     {
00807         // subpixel estimation
00808         bestRes = subpixelMaxima(vigra::srcImageRange(bestCorrelation), bestRes.maxpos.toDiff2D());
00809         DEBUG_DEBUG("subpixel position: max:" << bestRes.maxi
00810                     << " at " << bestRes.maxpos << " under angle: " << bestAngle/M_PI*180);
00811     } else {
00812         // not enough values for subpixel estimation.
00813         DEBUG_ERROR("subpixel estimation not done, maxima to close to border");
00814     }
00815 
00816     bestRes.maxpos = bestRes.maxpos + searchUL;
00817     bestRes.maxAngle = bestAngle/M_PI*180;
00818     return bestRes;
00819 }
00820 #else
00821 // specialed version for FFT based correlation,
00822 // optimized correlation calculation to save reuseable results
00823 template <class SrcImage, class DestImage, class KernelImage>
00824 CorrelationResult correlateImageFastRotationFFT(SrcImage & src, DestImage & dest, KernelImage & unrotatedKernel,
00825     vigra::Diff2D templPos, int templWidth, vigra::Diff2D tmplSize,
00826     double startAngle, double stopAngle, int angleSteps)
00827 {
00828     const int sw = src.width();
00829     const int sh = src.height();
00830     std::vector<CorrelationResult> results(angleSteps);
00831     std::vector<DestImage> resultsImg(angleSteps, DestImage(sw, sh));
00832 
00833     vigra::FFTWComplexImage spatialSearch(sw, sh);
00834     vigra::FFTWComplexImage fourierSearch(sw, sh);
00835     //forward FFTW plan
00836     fftw_plan forwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)spatialSearch.begin(), (fftw_complex *)fourierSearch.begin(), FFTW_FORWARD, FFTW_ESTIMATE);
00837     //FFT of search image, we need it for all angles
00838     vigra::copyImage(srcImageRange(src), destImage(spatialSearch, vigra::FFTWWriteRealAccessor<vigra::FImage::value_type>()));
00839     fftw_execute(forwardPlan);
00840     // generate backwardPlan for inplace use
00841     fftw_plan backwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)fourierSearch.begin(), (fftw_complex *)fourierSearch.begin(), FFTW_BACKWARD, FFTW_ESTIMATE);
00842 
00843     // calculate look up sum tables
00844     // are used by all angles
00845     // use double instead of float!, otherwise there can be truncation errors
00846     vigra::DImage s(sw, sh);
00847     vigra::DImage s2(sw, sh);
00848     double val = src(0, 0);
00849     s(0, 0) = val;
00850     s2(0, 0) = val * val;
00851     // special treatment for first line
00852     for (int x = 1; x < sw; ++x)
00853     {
00854         val = src(x, 0);
00855         s(x, 0) = val + s(x - 1, 0);
00856         s2(x, 0) = val * val + s2(x - 1, 0);
00857     }
00858     // special treatment for first column
00859     for (int y = 1; y < sh; ++y)
00860     {
00861         val = src(0, y);
00862         s(0, y) = val + s(0, y - 1);
00863         s2(0, y) = val * val + s2(0, y - 1);
00864     }
00865     // final summation
00866     for (int y = 1; y < sh; ++y)
00867     {
00868         for (int x = 1; x < sw; ++x)
00869         {
00870             val = src(x, y);
00871             s(x, y) = val + s(x - 1, y) + s(x, y - 1) - s(x - 1, y - 1);
00872             s2(x, y) = val * val + s2(x - 1, y) + s2(x, y - 1) - s2(x - 1, y - 1);
00873         };
00874     };
00875 
00876     const double step = (stopAngle - startAngle) / (angleSteps - 1);
00877     const int kw = tmplSize.x;
00878     const int kh = tmplSize.y;
00879     const int yend = sh - 2 * templWidth;
00880     const int xend = sw - 2 * templWidth;
00881 #pragma omp parallel for
00882     for (int i = 0; i < angleSteps; ++i)
00883     {
00884         vigra::FImage kernel(tmplSize);
00885         vigra::FImage alpha(kernel.size());
00886         vigra_ext::PassThroughFunctor<float> nf;
00887         AppBase::DummyProgressDisplay dummy;
00888 
00889         RotateTransform t(startAngle + i * step, hugin_utils::FDiff2D(templWidth, templWidth), templPos);
00890         vigra_ext::transformImage(srcImageRange(unrotatedKernel, vigra::RGBToGrayAccessor<typename KernelImage::value_type>()), destImageRange(kernel),
00891             destImage(alpha), vigra::Diff2D(0, 0), t, nf, false, vigra_ext::INTERP_CUBIC, &dummy, true);
00892 
00893         // calculate mean and variance of kernel/template
00894         vigra::FindAverageAndVariance<vigra::FImage::PixelType> kMean;
00895         vigra::inspectImage(srcImageRange(kernel), kMean);
00896         // uniform patch, skip calculation
00897         if (kMean.variance(false) == 0)
00898         {
00899             continue;
00900         };
00901         // subtract mean from kernel/template
00902         vigra::transformImage(srcImageRange(kernel), destImage(kernel), vigra::functor::Arg1() - vigra::functor::Param(kMean.average()));
00903 
00904         vigra::FFTWComplexImage complexKernel(sw, sh);
00905         vigra::FFTWComplexImage fourierKernel(sw, sh);
00906         // copy kernel to complex data structure
00907         vigra::copyImage(srcImageRange(kernel), destImage(complexKernel, vigra::FFTWWriteRealAccessor<vigra::FImage::value_type>()));
00908         fftw_execute_dft(forwardPlan, (fftw_complex*)complexKernel.begin(), (fftw_complex*)fourierKernel.begin());
00909 
00910         // multiply SrcImage with conjugated kernel in frequency domain
00911         vigra::combineTwoImages(srcImageRange(fourierSearch), srcImage(fourierKernel), destImage(fourierKernel), &multiplyConjugate<vigra::FFTWComplex<vigra::FImage::value_type> >);
00912         // FFT back into spatial domain (inplace)
00913         fftw_execute_dft(backwardPlan, (fftw_complex*)fourierKernel.begin(), (fftw_complex*)fourierKernel.begin());
00914 
00915         // calculate constant part
00916         const double normFactor = 1.0 / (sw * sh * sqrt(kMean.variance(false)));
00917         for (int yr = 0; yr < yend; ++yr)
00918         {
00919             for (int xr = 0; xr < xend; ++xr)
00920             {
00921                 double value = fourierKernel(xr, yr).re() * normFactor;
00922                 // do final summation using the lookup tables
00923                 double sumF = s(xr + kw - 1, yr + kh - 1);
00924                 double sumF2 = s2(xr + kw - 1, yr + kh - 1);
00925                 if (xr > 0)
00926                 {
00927                     sumF -= s(xr - 1, yr + kh - 1);
00928                     sumF2 -= s2(xr - 1, yr + kh - 1);
00929                 };
00930                 if (yr > 0)
00931                 {
00932                     sumF -= s(xr + kw - 1, yr - 1);
00933                     sumF2 -= s2(xr + kw - 1, yr - 1);
00934                 };
00935                 if (xr > 0 && yr > 0)
00936                 {
00937                     sumF += s(xr - 1, yr - 1);
00938                     sumF2 += s2(xr - 1, yr - 1);
00939                 };
00940 
00941                 double den = sqrt((kw * kh * sumF2 - sumF * sumF));
00942                 // prevent division through zero
00943                 if (den != 0)
00944                 {
00945                     value /= den;
00946                     if (value > results[i].maxi)
00947                     {
00948                         results[i].maxi = value;
00949                         results[i].maxpos.x = xr + templWidth;
00950                         results[i].maxpos.y = yr + templWidth;
00951                     }
00952                     resultsImg[i](xr + templWidth, yr + templWidth) = vigra::NumericTraits<typename DestImage::value_type>::fromRealPromote(value);
00953                 };
00954             };
00955         };
00956     };
00957     fftw_destroy_plan(forwardPlan);
00958     fftw_destroy_plan(backwardPlan);
00959     int maxIndex = 0;
00960     double maxValue = 0;
00961     for (size_t i = 0; i < results.size(); ++i)
00962     {
00963         if (results[i].maxi > maxValue)
00964         {
00965             maxIndex = i;
00966             maxValue = results[i].maxi;
00967         };
00968     };
00969     CorrelationResult res(results[maxIndex]);
00970     res.maxAngle = startAngle + maxIndex * step;
00971     dest = resultsImg[maxIndex];
00972     return res;
00973 };
00974 
00975 template <class IMAGET, class IMAGES>
00976 CorrelationResult PointFineTuneRotSearch(const IMAGET & templImg,
00977     vigra::Diff2D templPos,
00978     int templSize,
00979     const IMAGES & searchImg,
00980     vigra::Diff2D searchPos,
00981     int sWidth,
00982     double startAngle,
00983     double stopAngle,
00984     int angleSteps)
00985 {
00986     DEBUG_TRACE("templPos: " << templPos << " searchPos: " << searchPos);
00987 
00988     // extract patch from template
00989     int templWidth = templSize / 2;
00990     vigra::Diff2D tmplUL(-templWidth, -templWidth);
00991     vigra::Diff2D tmplLR(templWidth, templWidth);
00992     // clip template
00993     if (tmplUL.x + templPos.x < 0) tmplUL.x = -templPos.x;
00994     if (tmplUL.y + templPos.y < 0) tmplUL.y = -templPos.y;
00995     if (tmplLR.x + templPos.x> templImg.width())
00996         tmplLR.x = templImg.width() - templPos.x;
00997     if (tmplLR.y + templPos.y > templImg.height())
00998         tmplLR.y = templImg.height() - templPos.y;
00999     vigra::Diff2D tmplSize = tmplLR - tmplUL + vigra::Diff2D(1, 1);
01000 
01001     // extract patch from search region
01002     // make search region bigger, so that interpolation can always be done
01003     int swidth = sWidth / 2 + (2 + templWidth);
01004     DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
01005     // clip search window
01006     if (searchPos.x < 0) searchPos.x = 0;
01007     if (searchPos.x >(int) searchImg.width()) searchPos.x = searchImg.width() - 1;
01008     if (searchPos.y < 0) searchPos.y = 0;
01009     if (searchPos.y >(int) searchImg.height()) searchPos.x = searchImg.height() - 1;
01010 
01011     vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
01012     vigra::Diff2D searchLR(searchPos.x + swidth, searchPos.y + swidth);
01013     // clip search window
01014     if (searchUL.x < 0) searchUL.x = 0;
01015     if (searchUL.x > searchImg.width()) searchUL.x = searchImg.width();
01016     if (searchUL.y < 0) searchUL.y = 0;
01017     if (searchUL.y > searchImg.height()) searchUL.y = searchImg.height();
01018     if (searchLR.x > searchImg.width()) searchLR.x = searchImg.width();
01019     if (searchLR.x < 0) searchLR.x = 0;
01020     if (searchLR.y > searchImg.height()) searchLR.y = searchImg.height();
01021     if (searchLR.y < 0) searchLR.y = 0;
01022     DEBUG_DEBUG("search borders: " << searchUL << "," << searchLR << " size: " << searchLR - searchUL);
01023 
01024     // correlation output
01025     vigra::Diff2D searchSize = searchLR - searchUL;
01026     vigra::FImage dest(searchSize);
01027     dest.init(-1);
01028 
01029     // source input
01030     vigra::FImage srcImage(searchLR - searchUL);
01031     vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
01032         searchImg.upperLeft() + searchLR,
01033         vigra::RGBToGrayAccessor<typename IMAGES::value_type>()),
01034         destImage(srcImage));
01035 
01036     CorrelationResult resCorrelate=correlateImageFastRotationFFT(srcImage, dest, templImg, templPos, templWidth, tmplSize,
01037         startAngle, stopAngle, angleSteps);
01038     CorrelationResult res;
01039 
01040     DEBUG_DEBUG("rotation search finished, max:" << resCorrelate.maxi
01041         << " at " << resCorrelate.maxpos << " angle:" << RAD_TO_DEG(resCorrelate.maxAngle) << "");
01042 
01043     // do a subpixel maxima estimation
01044     // check if the max is inside the pixel boundaries,
01045     // and there are enought correlation values for the subpixel
01046     // estimation, (2 + templWidth)
01047     if (resCorrelate.maxpos.x > 1 + templWidth && resCorrelate.maxpos.x < 2 * swidth + 1 - 1 - templWidth
01048         && resCorrelate.maxpos.y > 1 + templWidth && resCorrelate.maxpos.y < 2 * swidth + 1 - 1 - templWidth)
01049     {
01050         // subpixel estimation
01051         res = subpixelMaxima(vigra::srcImageRange(dest), resCorrelate.maxpos.toDiff2D());
01052 
01053         DEBUG_DEBUG("subpixel position: max:" << res.maxi
01054             << " at " << res.maxpos);
01055     }
01056     else {
01057         // not enough values for subpixel estimation.
01058         res = resCorrelate;
01059         DEBUG_ERROR("subpixel estimation not done, maxima to close to border");
01060     }
01061 
01062     res.maxpos = res.maxpos + searchUL;
01063     res.maxAngle = RAD_TO_DEG(resCorrelate.maxAngle);
01064     return res;
01065 }
01066 #endif
01067 
01077 template <class SrcIterator, class SrcAccessor,
01078           class DestIterator, class DestAccessor,
01079           class KernelIterator, class KernelAccessor>
01080 CorrelationResult correlateImage(SrcIterator sul, SrcIterator slr, SrcAccessor as,
01081                                  DestIterator dul, DestAccessor ad,
01082                                  KernelIterator ki, KernelAccessor ak,
01083                                  vigra::Diff2D kul, vigra::Diff2D klr,
01084                                  double threshold = 0.7 )
01085 {
01086     CorrelationResult res;
01087     vigra_precondition(kul.x <= 0 && kul.y <= 0,
01088                  "convolveImage(): coordinates of "
01089                  "kernel's upper left must be <= 0.");
01090     vigra_precondition(klr.x >= 0 && klr.y >= 0,
01091                  "convolveImage(): coordinates of "
01092                  "kernel's lower right must be >= 0.");
01093 
01094     // use traits to determine SumType as to prevent possible overflow
01095     typedef typename
01096         vigra::NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
01097     typedef typename
01098         vigra::NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
01099     typedef
01100         vigra::NumericTraits<typename DestAccessor::value_type> DestTraits;
01101 
01102     // calculate width and height of the image
01103     int w = slr.x - sul.x;
01104     int h = slr.y - sul.y;
01105     int wk = klr.x - kul.x +1;
01106     int hk = klr.y - kul.y +1;
01107 
01108 //    DEBUG_DEBUG("correlate Image srcSize " << (slr - sul).x << "," << (slr - sul).y
01109 //                << " tmpl size: " << wk << "," << hk);
01110 
01111     vigra_precondition(w >= wk && h >= hk,
01112                  "convolveImage(): kernel larger than image.");
01113 
01114     int x,y;
01115     int ystart = -kul.y;
01116     int yend   = h-klr.y;
01117     int xstart = -kul.x;
01118     int xend   = w-klr.x;
01119 
01120     // calculate template mean
01121     vigra::FindAverage<typename KernelAccessor::value_type> average;
01122     vigra::inspectImage(ki + kul, ki + klr + vigra::Diff2D(1,1),
01123                         ak,
01124                         average);
01125     KSumType kmean = average();
01126 
01127     // create y iterators, they iterate over the rows.
01128     DestIterator yd = dul + vigra::Diff2D(xstart, ystart);
01129     SrcIterator ys = sul + vigra::Diff2D(xstart, ystart);
01130 
01131 
01132 //    DEBUG_DEBUG("size: " << w << "," <<  h << " ystart: " << ystart <<", yend: " << yend);
01133     for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y)
01134     {
01135 
01136         // create x iterators, they iterate the coorelation over
01137         // the columns
01138         DestIterator xd(yd);
01139         SrcIterator xs(ys);
01140 
01141         for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x)
01142         {
01143             if (as(xd) < threshold) {
01144                 continue;
01145             }
01146 
01147             // init the sum
01148             SumType numerator = 0;
01149             SumType div1 = 0;
01150             SumType div2 = 0;
01151             SumType spixel = 0;
01152             KSumType kpixel = 0;
01153 
01154             // create inner y iterators
01155             // access to the source image
01156             SrcIterator yys = xs + kul;
01157             // access to the kernel image
01158             KernelIterator yk  = ki + kul;
01159 
01160             vigra::FindAverage<typename SrcAccessor::value_type> average;
01161             vigra::inspectImage(xs + kul, xs + klr + vigra::Diff2D(1,1), as, average);
01162             double mean = average();
01163 
01164             int xx, yy;
01165             for(yy=0; yy<hk; ++yy, ++yys.y, ++yk.y)
01166             {
01167                 // create inner x iterators
01168                 SrcIterator xxs = yys;
01169                 KernelIterator xk = yk;
01170 
01171                 for(xx=0; xx<wk; ++xx, ++xxs.x, ++xk.x)
01172                 {
01173                     spixel = as(xxs) - mean;
01174                     kpixel = ak(xk) - kmean;
01175                     numerator += kpixel * spixel;
01176                     div1 += kpixel * kpixel;
01177                     div2 += spixel * spixel;
01178                 }
01179             }
01180             numerator = (numerator/sqrt(div1 * div2));
01181             if (numerator > res.maxi) {
01182                 res.maxi = numerator;
01183                 res.maxpos.x = x;
01184                 res.maxpos.y = y;
01185             }
01186             numerator = numerator;
01187             // store correlation in destination pixel
01188             ad.set(DestTraits::fromRealPromote(numerator), xd);
01189         }
01190     }
01191     return res;
01192 }
01193 
01194 } // namespace
01195 
01196 #endif // VIGRA_EXT_CORRELATION_H

Generated on 22 May 2018 for Hugintrunk by  doxygen 1.4.7