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 <vigra/impex.hxx>
00035 
00036 #include "hugin_utils/utils.h"
00037 #include "hugin_math/hugin_math.h"
00038 #include "vigra_ext/Pyramid.h"
00039 #include "vigra_ext/FitPolynom.h"
00040 #include "vigra_ext/utils.h"
00041 
00042 // hmm.. not really great.. should be part of vigra_ext
00043 #include "vigra_ext/ImageTransforms.h"
00044 
00045 #define VIGRA_EXT_USE_FAST_CORR
00046 
00047 //#define DEBUG_WRITE_FILES
00048 
00049 namespace vigra_ext{
00050 
00052 struct CorrelationResult
00053 {
00054     CorrelationResult()
00055         : maxi(-1), maxpos(0,0), curv(0,0), maxAngle(0)
00056         { }
00057     // value at correlation peak.
00058     double maxi;
00059     hugin_utils::FDiff2D maxpos;
00060     // curvature of the correlation peak
00061     hugin_utils::FDiff2D curv;
00062     double maxAngle;
00063 };
00064 
00077 template <class SrcImage, class DestImage, class KernelImage>
00078 CorrelationResult correlateImageFast(SrcImage & src,
00079                                      DestImage & dest,
00080                                      KernelImage & kernel,
00081                                      vigra::Diff2D kul, vigra::Diff2D klr,
00082                                      double threshold = 0.7 )
00083 {
00084     vigra_precondition(kul.x <= 0 && kul.y <= 0,
00085                  "convolveImage(): coordinates of "
00086                  "kernel's upper left must be <= 0.");
00087     vigra_precondition(klr.x >= 0 && klr.y >= 0,
00088                  "convolveImage(): coordinates of "
00089                  "kernel's lower right must be >= 0.");
00090 
00091     // use traits to determine SumType as to prevent possible overflow
00092     typedef typename
00093         vigra::NumericTraits<typename SrcImage::value_type>::RealPromote SumType;
00094     typedef typename
00095         vigra::NumericTraits<typename KernelImage::value_type>::RealPromote KSumType;
00096     typedef
00097         vigra::NumericTraits<typename DestImage::value_type> DestTraits;
00098 
00099     // calculate width and height of the image
00100     int w = src.width();
00101     int h = src.height();
00102     int wk = kernel.width();
00103     int hk = kernel.height();
00104 
00105 //    DEBUG_DEBUG("correlate Image srcSize " << (slr - sul).x << "," << (slr - sul).y
00106 //                << " tmpl size: " << wk << "," << hk);
00107 
00108     vigra_precondition(w >= wk && h >= hk,
00109                  "convolveImage(): kernel larger than image.");
00110 
00111 //    int x,y;
00112     int ystart = -kul.y;
00113     int yend   = h-klr.y;
00114     int xstart = -kul.x;
00115     int xend   = w-klr.x;
00116 
00117     // mean of kernel
00118     KSumType kmean=0;
00119     for(int y=0; y < hk; y++) {
00120         for(int x=0; x < wk; x++) {
00121             kmean += kernel(x,y);
00122         }
00123     }
00124     kmean = kmean / (hk*wk);
00125 
00126     CorrelationResult res;
00127 
00128     // create y iterators, they iterate over the rows.
00129 //    DestIterator yd = dul + vigra::Diff2D(xstart, ystart);
00130 //    SrcIterator ys = sul + vigra::Diff2D(xstart, ystart);
00131 
00132 
00133 //    DEBUG_DEBUG("size: " << w << "," <<  h << " ystart: " << ystart <<", yend: " << yend);
00134     int unifwarned=0;
00135     for(int yr=ystart; yr < yend; ++yr)
00136     {
00137         // create x iterators, they iterate the coorelation over
00138         // the columns
00139 //        DestIterator xd(yd);
00140 //        SrcIterator xs(ys);
00141 
00142         for(int xr=xstart; xr < xend; ++xr)
00143         {
00144             if (dest(xr,yr) < threshold) {
00145                 continue;
00146             }
00147 //            int x0, y0, x1, y1;
00148 
00149 //            y0 = kul.y;
00150 //            y1 = klr.y;
00151 //            x0 = kul.x;
00152 //            x1 = klr.x;;
00153 
00154             // init the sum
00155             SumType numerator = 0;
00156             SumType div1 = 0;
00157             SumType div2 = 0;
00158             SumType spixel = 0;
00159             KSumType kpixel = 0;
00160 
00161             // create inner y iterators
00162             // access to the source image
00163 //            SrcIterator yys = xs + kul;
00164             // access to the kernel image
00165 //            KernelIterator yk  = ki + kul;
00166 
00167             // mean of image patch
00168             KSumType mean=0;
00169                         int ym;
00170             for(ym=yr+kul.y; ym <= yr+klr.y; ym++) {
00171                 for(int xm=xr+kul.x; xm <= xr+klr.x; xm++) {
00172                     mean += src(xm,ym);
00173                 }
00174             }
00175             mean = mean / (hk*wk);
00176 
00177             // perform correlation (inner loop)
00178             ym=yr+kul.y;
00179             int yk;
00180             for(yk=0; yk < hk; yk++) {
00181                 int xm=xr+kul.x;
00182                 int xk;
00183                 for(xk=0; xk < wk; xk++) {
00184                     spixel = src(xm,ym) - mean;
00185                     kpixel = kernel(xk,yk) - kmean;
00186                     numerator += kpixel * spixel;
00187                     div1 += kpixel * kpixel;
00188                     div2 += spixel * spixel;
00189                     xm++;
00190                 }
00191                 ym++;
00192             }
00193             if (div1*div2 == 0) {
00194                 // This happens when one of the patches is perfectly uniform
00195                 numerator = 0;  // Set correlation to zero since this is uninteresting
00196                 if (!unifwarned) {
00197                     DEBUG_DEBUG("Uniform patch(es) during correlation computation");
00198                     unifwarned=1;
00199                 }
00200             } else
00201                 numerator = (numerator/sqrt(div1 * div2));
00202             
00203             if (numerator > res.maxi) {
00204                 res.maxi = numerator;
00205                 res.maxpos.x = xr;
00206                 res.maxpos.y = yr;
00207             }
00208             dest(xr,yr) = DestTraits::fromRealPromote(numerator);
00209         }
00210     }
00211     return res;
00212 }
00213 
00221 template <class Iterator, class Accessor>
00222 CorrelationResult subpixelMaxima(vigra::triple<Iterator, Iterator, Accessor> img,
00223                                  vigra::Diff2D max)
00224 {
00225     const int interpWidth = 1;
00226     CorrelationResult res;
00227     vigra_precondition(max.x > interpWidth && max.y > interpWidth,
00228                  "subpixelMaxima(): coordinates of "
00229                  "maxima must be > interpWidth, interpWidth.");
00230     vigra::Diff2D sz = img.second - img.first;
00231     vigra_precondition(sz.x - max.x >= interpWidth && sz.y - max.y >= interpWidth,
00232                  "subpixelMaxima(): coordinates of "
00233                  "maxima must interpWidth pixels from the border.");
00234     typedef typename Accessor::value_type T;
00235 
00236     T x[2*interpWidth+1];
00237     T zx[2*interpWidth+1];
00238     T zy[2*interpWidth+1];
00239 
00240 #ifdef DEBUG_CORRELATION
00241     exportImage(img,vigra::ImageExportInfo("test.tif"));
00242 #endif
00243     
00244     Accessor acc = img.third;
00245     Iterator begin=img.first;
00246     for (int i=-interpWidth; i<=interpWidth; i++) {
00247         // collect first row
00248         x[i+interpWidth] = i;
00249         zx[i+interpWidth] = acc(begin, max + vigra::Diff2D(i,0));
00250         zy[i+interpWidth] = acc(begin, max + vigra::Diff2D(0,i));
00251     }
00252 
00253     double a,b,c;
00254     FitPolynom(x, x + 2*interpWidth+1, zx, a,b,c);
00255     if (hugin_utils::isnan(a) || hugin_utils::isnan(b) || hugin_utils::isnan(c)) {
00256         exportImage(img,vigra::ImageExportInfo("test.tif"));
00257         DEBUG_ERROR("Bad polynomial fit results");
00258         res.maxpos.x=max.x;
00259         res.maxpos.y=max.y;
00260         return res;
00261     }
00262 
00263     // calculate extrema of x position by setting
00264     // the 1st derivate to zero
00265     // 2*c*x + b = 0
00266     if (c==0)
00267         res.maxpos.x=0;
00268     else
00269         res.maxpos.x = -b/(2*c);
00270 
00271     res.curv.x = 2*c;
00272 
00273     // calculate result at maxima
00274     double maxx = c*res.maxpos.x*res.maxpos.x + b*res.maxpos.x + a;
00275 
00276     FitPolynom(x, x + 2*interpWidth+1, zy, a,b,c);
00277     // calculate extrema of y position
00278     if (c==0)
00279         res.maxpos.y=0;
00280     else
00281         res.maxpos.y = -b/(2*c);
00282 
00283     res.curv.y = 2*c;
00284     double maxy = c*res.maxpos.y*res.maxpos.y + b*res.maxpos.y + a;
00285 
00286     // use mean of both maxima as new interpolation result
00287     res.maxi = (maxx + maxy) / 2;
00288     DEBUG_DEBUG("value at subpixel maxima (" << res.maxpos.x << " , "
00289                 <<res.maxpos.y << "): " << maxx << "," << maxy
00290                 << " mean:" << res.maxi << " coeff: a=" << a
00291                 << "; b=" << b << "; c=" << c);
00292 
00293     // test if the point has moved too much ( more than 1.5 pixel).
00294     // actually, I should also test the confidence of the fit.
00295     if (fabs(res.maxpos.x) > 1 || fabs(res.maxpos.y) > 1) {
00296         DEBUG_NOTICE("subpixel Maxima has moved to much, ignoring");
00297         res.maxpos.x = max.x;
00298         res.maxpos.y = max.y;
00299     } else {
00300         res.maxpos.x = res.maxpos.x + max.x;
00301         res.maxpos.y = res.maxpos.y + max.y;
00302     }
00303     return res;
00304 }
00305 
00308 class RotateTransform
00309 {
00310 public:
00311     RotateTransform(double phi, hugin_utils::FDiff2D origin, hugin_utils::FDiff2D transl)
00312         : m_phi(phi), m_origin(origin), m_transl(transl)
00313         { }
00314 
00315     bool transformImgCoord(double &destx, double &desty, double srcx, double srcy) const
00316     {
00317         srcx -= m_origin.x;
00318         srcy -= m_origin.y;
00319 
00320         destx = srcx * cos(m_phi) + srcy * sin(m_phi);
00321         desty = srcx * -sin(m_phi) + srcy * cos(m_phi);
00322 
00323         destx += m_transl.x;
00324         desty += m_transl.y;
00325         return true;
00326     }
00327 
00328     double m_phi;
00329     hugin_utils::FDiff2D m_origin;
00330     hugin_utils::FDiff2D m_transl;
00331 };
00332 
00343 template <class IMAGET, class IMAGES>
00344 CorrelationResult PointFineTune(const IMAGET & templImg,
00345                                 vigra::Diff2D templPos,
00346                                 int templSize,
00347                                 const IMAGES & searchImg,
00348                                 vigra::Diff2D searchPos,
00349                                 int sWidth)
00350 {
00351 //    DEBUG_TRACE("templPos: " vigra::<< templPos << " searchPos: " vigra::<< searchPos);
00352 
00353     // extract patch from template
00354 
00355     int templWidth = templSize/2;
00356     vigra::Diff2D tmplUL(templPos.x - templWidth, templPos.y-templWidth);
00357     // lower right iterators "are past the end"
00358     vigra::Diff2D tmplLR(templPos.x + templWidth + 1, templPos.y + templWidth + 1);
00359     // clip corners to ensure the template is inside the image.
00360     vigra::Diff2D tmplImgSize(templImg.size());
00361     tmplUL = hugin_utils::simpleClipPoint(tmplUL, vigra::Diff2D(0,0), tmplImgSize);
00362     tmplLR = hugin_utils::simpleClipPoint(tmplLR, vigra::Diff2D(0,0), tmplImgSize);
00363     vigra::Diff2D tmplSize = tmplLR - tmplUL;
00364     DEBUG_DEBUG("template position: " << templPos << "  tmplUL: " << tmplUL
00365                     << "  templLR:" << tmplLR << "  size:" << tmplSize);
00366 
00367     // extract patch from search region
00368     // make search region bigger, so that interpolation can always be done
00369     int swidth = sWidth/2 +(2+templWidth);
00370 //    DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
00371 //    Diff2D subjPoint(searchPos);
00372     // clip search window
00373     if (searchPos.x < 0) searchPos.x = 0;
00374     if (searchPos.x > (int) searchImg.width()) searchPos.x = searchImg.width()-1;
00375     if (searchPos.y < 0) searchPos.y = 0;
00376     if (searchPos.y > (int) searchImg.height()) searchPos.x = searchImg.height()-1;
00377 
00378     vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
00379     // point past the end
00380     vigra::Diff2D searchLR(searchPos.x + swidth+1, searchPos.y + swidth+1);
00381     // clip search window
00382     vigra::Diff2D srcImgSize(searchImg.size());
00383     searchUL = hugin_utils::simpleClipPoint(searchUL, vigra::Diff2D(0,0), srcImgSize);
00384     searchLR = hugin_utils::simpleClipPoint(searchLR, vigra::Diff2D(0,0), srcImgSize);
00385 //    DEBUG_DEBUG("search borders: " << searchLR.x << "," << searchLR.y);
00386 
00387     vigra::Diff2D searchSize = searchLR - searchUL;
00388     // create output image
00389 
00390 //#ifdef VIGRA_EXT_USE_FAST_CORR
00391     // source input
00392     vigra::FImage srcImage(searchLR-searchUL);
00393     vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
00394                                         searchImg.upperLeft() + searchLR,
00395                                         vigra::RGBToGrayAccessor<typename IMAGES::value_type>() ),
00396                      destImage(srcImage) );
00397 
00398     vigra::FImage templateImage(tmplSize);
00399     vigra::copyImage(vigra::make_triple(templImg.upperLeft() + tmplUL,
00400                                         templImg.upperLeft() + tmplLR,
00401                                         vigra::RGBToGrayAccessor<typename IMAGET::value_type>()),
00402                      destImage(templateImage));
00403 #ifdef DEBUG_WRITE_FILES
00404     vigra::ImageExportInfo tmpli("hugin_templ.tif");
00405     vigra::exportImage(vigra::srcImageRange(templateImage), tmpli);
00406 
00407     vigra::ImageExportInfo srci("hugin_searchregion.tif");
00408     vigra::exportImage(vigra::srcImageRange(srcImage), srci);
00409 #endif
00410 
00411 //#endif
00412 
00413     vigra::FImage dest(searchSize);
00414     dest.init(-1);
00415     // we could use the multiresolution version as well.
00416     // but usually the region is quite small.
00417     CorrelationResult res;
00418 #ifdef VIGRA_EXT_USE_FAST_CORR
00419     DEBUG_DEBUG("+++++ starting fast correlation");
00420     res = correlateImageFast(srcImage,
00421                              dest,
00422                              templateImage,
00423                              tmplUL-templPos, tmplLR-templPos - vigra::Diff2D(1,1),
00424                              -1);
00425 #else
00426     DEBUG_DEBUG("+++++ starting normal correlation");
00427     res = correlateImage(srcImage.upperLeft(),
00428                          srcImage.lowerRight(),
00429                          srcImage.accessor(),
00430                          dest.upperLeft(),
00431                          dest.accessor(),
00432                          templateImage.upperLeft() + templPos,
00433                          templateImage.accessor(),
00434                          tmplUL, tmplLR, -1);
00435 
00436 //     res = correlateImage(searchImg.upperLeft() + searchUL,
00437 //                          searchImg.upperLeft() + searchLR,
00438 //                          searchImg.accessor(),
00439 //                          dest.upperLeft(),
00440 //                          dest.accessor(),
00441 //                          templImg.upperLeft() + templPos,
00442 //                          templImg.accessor(),
00443 //                          tmplUL, tmplLR, -1);
00444 #endif
00445     DEBUG_DEBUG("normal search finished, max:" << res.maxi
00446                 << " at " << res.maxpos);
00447     // do a subpixel maxima estimation
00448     // check if the max is inside the pixel boundaries,
00449     // and there are enought correlation values for the subpixel
00450     // estimation, (2 + templWidth)
00451     if (res.maxpos.x > 2 + templWidth && res.maxpos.x < 2*swidth+1-2-templWidth
00452         && res.maxpos.y > 2+templWidth && res.maxpos.y < 2*swidth+1-2-templWidth)
00453     {
00454         // subpixel estimation
00455         res = subpixelMaxima(vigra::srcImageRange(dest), res.maxpos.toDiff2D());
00456         DEBUG_DEBUG("subpixel position: max:" << res.maxi
00457                     << " at " << res.maxpos);
00458     } else {
00459         // not enough values for subpixel estimation.
00460         DEBUG_DEBUG("subpixel estimation not done, maxima too close to border");
00461     }
00462 
00463     res.maxpos = res.maxpos + searchUL;
00464     return res;
00465 }
00466 
00477 template <class IMAGET, class IMAGES>
00478 CorrelationResult PointFineTuneRotSearch(const IMAGET & templImg,
00479                                          vigra::Diff2D templPos,
00480                                          int templSize,
00481                                          const IMAGES & searchImg,
00482                                          vigra::Diff2D searchPos,
00483                                          int sWidth,
00484                                          double startAngle,
00485                                          double stopAngle,
00486                                          int angleSteps)
00487 {
00488     DEBUG_TRACE("templPos: " << templPos << " searchPos: " << searchPos);
00489 
00490     // extract patch from template
00491 
00492     // make template size user configurable as well?
00493     int templWidth = templSize/2;
00494     vigra::Diff2D tmplUL(-templWidth, -templWidth);
00495     vigra::Diff2D tmplLR(templWidth, templWidth);
00496     // clip template
00497     if (tmplUL.x + templPos.x < 0) tmplUL.x = -templPos.x;
00498     if (tmplUL.y + templPos.y < 0) tmplUL.y = -templPos.y;
00499     if (tmplLR.x + templPos.x> templImg.width())
00500         tmplLR.x = templImg.width() - templPos.x;
00501     if (tmplLR.y + templPos.y > templImg.height())
00502         tmplLR.y = templImg.height() - templPos.y;
00503     vigra::Diff2D tmplSize = tmplLR - tmplUL + vigra::Diff2D(1,1);
00504 
00505     // extract patch from search region
00506     // make search region bigger, so that interpolation can always be done
00507     int swidth = sWidth/2 +(2+templWidth);
00508     DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
00509 //    Diff2D subjPoint(searchPos);
00510     // clip search window
00511     if (searchPos.x < 0) searchPos.x = 0;
00512     if (searchPos.x > (int) searchImg.width()) searchPos.x = searchImg.width()-1;
00513     if (searchPos.y < 0) searchPos.y = 0;
00514     if (searchPos.y > (int) searchImg.height()) searchPos.x = searchImg.height()-1;
00515 
00516     vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
00517     vigra::Diff2D searchLR(searchPos.x + swidth, searchPos.y + swidth);
00518     // clip search window
00519     if (searchUL.x < 0) searchUL.x = 0;
00520     if (searchUL.x > searchImg.width()) searchUL.x = searchImg.width();
00521     if (searchUL.y < 0) searchUL.y = 0;
00522     if (searchUL.y > searchImg.height()) searchUL.y = searchImg.height();
00523     if (searchLR.x > searchImg.width()) searchLR.x = searchImg.width();
00524     if (searchLR.x < 0) searchLR.x = 0;
00525     if (searchLR.y > searchImg.height()) searchLR.y = searchImg.height();
00526     if (searchLR.y < 0) searchLR.y = 0;
00527     DEBUG_DEBUG("search borders: " << searchUL << "," << searchLR << " size: " << searchLR -searchUL);
00528 
00529 
00530 #ifdef DEBUG_WRITE_FILES
00531     DEBUG_DEBUG("template area: " << templPos + tmplUL << " -> " << templPos + tmplLR);
00532     vigra::exportImage(vigra::make_triple(templImg.upperLeft() + templPos + tmplUL,
00533                                           templImg.upperLeft() + templPos + tmplLR + vigra::Diff2D(1,1),
00534                                           templImg.accessor()),
00535                        vigra::ImageExportInfo("00_original_template.png"));
00536 
00537     vigra::exportImage(make_triple(searchImg.upperLeft() + searchUL,
00538                                    searchImg.upperLeft() + searchLR,
00539                                    searchImg.accessor()),
00540                        vigra::ImageExportInfo("00_searcharea.png"));
00541 #endif
00542 
00543     // rotated template
00544     vigra::FImage rotTemplate(tmplSize);
00545     vigra::FImage alpha(tmplSize);
00546 
00547     // correlation output
00548     vigra::Diff2D searchSize = searchLR - searchUL;
00549     vigra::FImage dest(searchSize);
00550     dest.init(1);
00551     vigra::FImage bestCorrelation(searchSize);
00552 
00553 //#ifdef VIGRA_EXT_USE_FAST_CORR
00554     // source input
00555     vigra::FImage srcImage(searchLR-searchUL);
00556     vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
00557                                         searchImg.upperLeft() + searchLR,
00558                                         vigra::RGBToGrayAccessor<typename IMAGES::value_type>() ),
00559                      destImage(srcImage) );
00560 //#endif
00561 
00562     CorrelationResult bestRes;
00563     bestRes.maxi = -1;
00564     double bestAngle = 0;
00565 
00566     AppBase::MultiProgressDisplay dummy;
00567     // test the image at rotation angles with 30 deg. steps.
00568     double step = (stopAngle - startAngle)/(angleSteps-1);
00569     double phi=startAngle;
00570     vigra_ext::PassThroughFunctor<float> nf;
00571     for (double i=0; phi <= stopAngle; i++, phi += step) {
00572         DEBUG_DEBUG("+++++ Rotating image, phi: " << RAD_TO_DEG(phi));
00573         RotateTransform t(phi, hugin_utils::FDiff2D(templWidth, templWidth), templPos);
00574         vigra_ext::transformImage(srcImageRange(templImg, vigra::RGBToGrayAccessor<typename IMAGET::value_type>()),
00575                            destImageRange(rotTemplate),
00576                            destImage(alpha),
00577                            vigra::Diff2D(0,0),
00578                            t,
00579                            nf,
00580                            false,
00581                            vigra_ext::INTERP_CUBIC,
00582                            dummy);
00583         DEBUG_DEBUG("----- Image rotated");
00584 
00585         // force a search in at all points.
00586         dest.init(-1);
00587         DEBUG_DEBUG("+++++ starting correlation");
00588 
00589         CorrelationResult res;
00590         // we could use the multiresolution version as well.
00591         // but usually the region is quite small.
00592 #ifdef VIGRA_EXT_USE_FAST_CORR
00593         res = correlateImageFast(srcImage,
00594                                  dest,
00595                                  rotTemplate,
00596                                  vigra::Diff2D(-templWidth, -templWidth),
00597                                  vigra::Diff2D(templWidth, templWidth),
00598                                  -1);
00599 #else
00600         res = correlateImage(srcImage.upperLeft(),
00601                              srcImage.lowerRight(),
00602                              srcImage.accessor(),
00603                              dest.upperLeft(),
00604                              dest.accessor(),
00605                              rotTemplate.upperLeft() + vigra::Diff2D(templWidth, templWidth),
00606                              rotTemplate.accessor(),
00607                              vigra::Diff2D(-templWidth, -templWidth), vigra::Diff2D(templWidth, templWidth),
00608                              -1);
00609 #endif
00610         DEBUG_DEBUG("---- correlation finished");
00611 
00612 #ifdef DEBUG_WRITE_FILES
00613         char fname[256];
00614         vigra::BImage tmpImg(rotTemplate.size());
00615         vigra::copyImage(srcImageRange(rotTemplate),destImage(tmpImg));
00616         sprintf(fname, "%3.0f_rotated_template.png", phi/M_PI*180);
00617         vigra::exportImage(srcImageRange(tmpImg), vigra::ImageExportInfo(fname));
00618 
00619         vigra::transformImage(vigra::srcImageRange(dest), vigra::destImage(dest),
00620                               vigra::linearRangeMapping(
00621                                   -1, 1,               // src range
00622                                   (unsigned char)0, (unsigned char)255) // dest range
00623             );
00624         tmpImg.resize(dest.size());
00625         vigra::copyImage(srcImageRange(dest),destImage(tmpImg));
00626         sprintf(fname, "%3.0f_corr_result.png", phi/M_PI*180);
00627         vigra::exportImage(srcImageRange(tmpImg), vigra::ImageExportInfo(fname));
00628 #endif
00629 
00630 
00631         DEBUG_DEBUG("normal search finished, max:" << res.maxi
00632                     << " at " << res.maxpos << " angle:" << phi/M_PI*180 << "");
00633 
00634         if (res.maxi > bestRes.maxi) {
00635             // remember best correlation.
00636             bestCorrelation = dest;
00637             bestRes = res;
00638             bestAngle = phi;
00639         }
00640 
00641     }
00642 
00643     DEBUG_DEBUG("rotation search finished, max:" << bestRes.maxi
00644                 << " at " << bestRes.maxpos << " angle:" << bestAngle/M_PI*180 << "");
00645 
00646     // do a subpixel maxima estimation
00647     // check if the max is inside the pixel boundaries,
00648     // and there are enought correlation values for the subpixel
00649     // estimation, (2 + templWidth)
00650     if (bestRes.maxpos.x > 1 + templWidth && bestRes.maxpos.x < 2*swidth+1-1-templWidth
00651         && bestRes.maxpos.y > 1+templWidth && bestRes.maxpos.y < 2*swidth+1-1-templWidth)
00652     {
00653         // subpixel estimation
00654         bestRes = subpixelMaxima(vigra::srcImageRange(bestCorrelation), bestRes.maxpos.toDiff2D());
00655         DEBUG_DEBUG("subpixel position: max:" << bestRes.maxi
00656                     << " at " << bestRes.maxpos << " under angle: " << bestAngle/M_PI*180);
00657     } else {
00658         // not enough values for subpixel estimation.
00659         DEBUG_ERROR("subpixel estimation not done, maxima to close to border");
00660     }
00661 
00662     bestRes.maxpos = bestRes.maxpos + searchUL;
00663     bestRes.maxAngle = bestAngle/M_PI*180;
00664     return bestRes;
00665 }
00666 
00676 template <class SrcIterator, class SrcAccessor,
00677           class DestIterator, class DestAccessor,
00678           class KernelIterator, class KernelAccessor>
00679 CorrelationResult correlateImage(SrcIterator sul, SrcIterator slr, SrcAccessor as,
00680                                  DestIterator dul, DestAccessor ad,
00681                                  KernelIterator ki, KernelAccessor ak,
00682                                  vigra::Diff2D kul, vigra::Diff2D klr,
00683                                  double threshold = 0.7 )
00684 {
00685     CorrelationResult res;
00686     vigra_precondition(kul.x <= 0 && kul.y <= 0,
00687                  "convolveImage(): coordinates of "
00688                  "kernel's upper left must be <= 0.");
00689     vigra_precondition(klr.x >= 0 && klr.y >= 0,
00690                  "convolveImage(): coordinates of "
00691                  "kernel's lower right must be >= 0.");
00692 
00693     // use traits to determine SumType as to prevent possible overflow
00694     typedef typename
00695         vigra::NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
00696     typedef typename
00697         vigra::NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
00698     typedef
00699         vigra::NumericTraits<typename DestAccessor::value_type> DestTraits;
00700 
00701     // calculate width and height of the image
00702     int w = slr.x - sul.x;
00703     int h = slr.y - sul.y;
00704     int wk = klr.x - kul.x +1;
00705     int hk = klr.y - kul.y +1;
00706 
00707 //    DEBUG_DEBUG("correlate Image srcSize " << (slr - sul).x << "," << (slr - sul).y
00708 //                << " tmpl size: " << wk << "," << hk);
00709 
00710     vigra_precondition(w >= wk && h >= hk,
00711                  "convolveImage(): kernel larger than image.");
00712 
00713     int x,y;
00714     int ystart = -kul.y;
00715     int yend   = h-klr.y;
00716     int xstart = -kul.x;
00717     int xend   = w-klr.x;
00718 
00719     // calculate template mean
00720     vigra::FindAverage<typename KernelAccessor::value_type> average;
00721     vigra::inspectImage(ki + kul, ki + klr + vigra::Diff2D(1,1),
00722                         ak,
00723                         average);
00724     KSumType kmean = average();
00725 
00726     // create y iterators, they iterate over the rows.
00727     DestIterator yd = dul + vigra::Diff2D(xstart, ystart);
00728     SrcIterator ys = sul + vigra::Diff2D(xstart, ystart);
00729 
00730 
00731 //    DEBUG_DEBUG("size: " << w << "," <<  h << " ystart: " << ystart <<", yend: " << yend);
00732     for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y)
00733     {
00734 
00735         // create x iterators, they iterate the coorelation over
00736         // the columns
00737         DestIterator xd(yd);
00738         SrcIterator xs(ys);
00739 
00740         for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x)
00741         {
00742             if (as(xd) < threshold) {
00743                 continue;
00744             }
00745 //            int x0, y0, x1, y1;
00746 
00747 //            y0 = kul.y;
00748 //            y1 = klr.y;
00749 //            x0 = kul.x;
00750 //            x1 = klr.x;;
00751 
00752             // init the sum
00753             SumType numerator = 0;
00754             SumType div1 = 0;
00755             SumType div2 = 0;
00756             SumType spixel = 0;
00757             KSumType kpixel = 0;
00758 
00759             // create inner y iterators
00760             // access to the source image
00761             SrcIterator yys = xs + kul;
00762             // access to the kernel image
00763             KernelIterator yk  = ki + kul;
00764 
00765             vigra::FindAverage<typename SrcAccessor::value_type> average;
00766             vigra::inspectImage(xs + kul, xs + klr + vigra::Diff2D(1,1), as, average);
00767             double mean = average();
00768 
00769             int xx, yy;
00770             for(yy=0; yy<hk; ++yy, ++yys.y, ++yk.y)
00771             {
00772                 // create inner x iterators
00773                 SrcIterator xxs = yys;
00774                 KernelIterator xk = yk;
00775 
00776                 for(xx=0; xx<wk; ++xx, ++xxs.x, ++xk.x)
00777                 {
00778                     spixel = as(xxs) - mean;
00779                     kpixel = ak(xk) - kmean;
00780                     numerator += kpixel * spixel;
00781                     div1 += kpixel * kpixel;
00782                     div2 += spixel * spixel;
00783                 }
00784             }
00785             numerator = (numerator/sqrt(div1 * div2));
00786             if (numerator > res.maxi) {
00787                 res.maxi = numerator;
00788                 res.maxpos.x = x;
00789                 res.maxpos.y = y;
00790             }
00791             numerator = numerator;
00792             // store correlation in destination pixel
00793             ad.set(DestTraits::fromRealPromote(numerator), xd);
00794         }
00795     }
00796     return res;
00797 }
00798 
00799 
00800 
00801 #if 0
00802 
00805 template <class Image>
00806 bool findTemplate(const Image & templ, const std::string & filename,
00807                   CorrelationResult & res)
00808 {
00809 //    DEBUG_TRACE("");
00810     // calculate pyramid level based on template size
00811     // the template should have be at least 1 pixel wide or high
00812     // (in the lowest resolution).
00813     int tw = templ.width();
00814     int th = templ.height();
00815     int ts = th<tw ? th : tw;
00816 
00817     // calculate the lowest level that will make sense
00818     int maxLevel = (int)floor(sqrt((double)ts/4.0));
00819     int levels = maxLevel+1;
00820     DEBUG_DEBUG("starting search on pyramid level " << maxLevel);
00821 
00822         vigra::BImage *templs = new vigra::BImage[levels];
00823     templs[0].resize(templ.width(), templ.height());
00824     vigra::copyImage(srcImageRange(templ), destImage(templs[0]));
00825 
00826     // create other levels
00827     for(int i=1; i<levels; i++) {
00828         vigra_ext::reduceToNextLevel(templs[i-1], templs[i]);
00829     }
00830 
00831     // this image will store the correlation coefficients
00832     // it will also be used by correlateImage as a mask image.
00833     // only the correlation for pixels with resImage(x,y) > theshold
00834     // will be calcuated. this avoids searching in completely uninteresting
00835     // image parts. (saves a lot of time)
00836 
00837     vigra::FImage * prevMask = new vigra::FImage(2,2,0.9);
00838     vigra::FImage * currentMask = new vigra::FImage();
00839     // start matching with highest level (lowest resolution images)
00840     float threshold=0;
00841     for (int i=maxLevel; i>=0; i--) {
00842 //    for (int i=0; i>=0; i--) {
00843         DEBUG_DEBUG("correlation level " << i);
00844         std::stringstream t;
00845         t << filename << "_" << i << "_";
00846         std::string basename = t.str();
00847 //        const vigra::BImage & s = ImageCache::getInstance().getPyramidImage(filename, i);
00848         DEBUG_DEBUG("subject size: " << s.width() << "," << s.height()
00849                     << "template size: " << templs[i].width() << "," << templs[i].height());
00850 //        saveImage(templs[i], basename + "template.pnm");
00851 //        saveImage(s, basename + "subject.pnm");
00852         currentMask->resize(s.width(),s.height());
00853         // scale up dest/mask Image
00854         // probably this should also use a threshold and extend the
00855         // resulting mask to the neighborhood pixels, just to catch unfortunate
00856         // correlations.
00857         vigra::resizeImageNoInterpolation(vigra::srcImageRange(*prevMask),vigra::destImageRange(*currentMask));
00858 
00859 //        saveScaledImage(*currentMask, -1, 1, basename + "mask.pnm");
00860 
00861         res = correlateImage(s.upperLeft(),
00862                              s.lowerRight(),
00863                              vigra::StandardValueAccessor<unsigned char>(),
00864                              currentMask->upperLeft(),
00865                              vigra::StandardValueAccessor<float>(),
00866                              templs[i].upperLeft(),
00867                              vigra::StandardValueAccessor<unsigned char>(),
00868                              vigra::Diff2D(0,0),
00869                              templs[i].size() - vigra::Diff2D(1,1),
00870                              threshold
00871             );
00872 //        saveScaledImage(*currentMask, -1, 1, basename + "result.pnm");
00873         DEBUG_DEBUG("Correlation result at level " << i << ":  max:" << res.maxi << " at: " << res.maxpos);
00874         // simple adaptive threshold.
00875         // FIXME, make this configurable? or select it according to some
00876         // statistical value
00877         threshold = res.maxi - 0.07;
00878 
00879         vigra::FImage *tmp = prevMask;
00880         prevMask = currentMask;
00881         currentMask = tmp;
00882     }
00883     delete prevMask;
00884     delete currentMask;
00885     return true;
00886 }
00887 #endif
00888 
00889 
00890 } // namespace
00891 
00892 #endif // VIGRA_EXT_CORRELATION_H

Generated on Thu Sep 18 01:25:37 2014 for Hugintrunk by  doxygen 1.3.9.1