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

Generated on 27 Jul 2015 for Hugintrunk by  doxygen 1.4.7