Hugintrunk  0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Correlation.h
Go to the documentation of this file.
1 // -*- c-basic-offset: 4 -*-
24 #ifndef VIGRA_EXT_CORRELATION_H
25 #define VIGRA_EXT_CORRELATION_H
26 
27 #include <hugin_shared.h>
28 #include <vigra/stdimage.hxx>
29 #include <vigra/inspectimage.hxx>
30 #include <vigra/copyimage.hxx>
31 #include <vigra/resizeimage.hxx>
32 #include <vigra/transformimage.hxx>
33 
34 #include "hugin_utils/utils.h"
35 #include "hugin_math/hugin_math.h"
36 #include "vigra_ext/FitPolynom.h"
37 #include "vigra_ext/utils.h"
38 
39 // hmm.. not really great.. should be part of vigra_ext
41 #include "hugin_config.h"
42 #ifdef HAVE_FFTW
43 #include <vigra/fftw3.hxx>
44 #include <vigra/functorexpression.hxx>
46 #include <vector>
47 #else
48 #define VIGRA_EXT_USE_FAST_CORR
49 #endif
50 
51 //#define DEBUG_WRITE_FILES
52 
53 namespace vigra_ext{
54 
57 {
59  : maxi(-1), maxpos(0,0), curv(0,0), maxAngle(0)
60  { }
61  // value at correlation peak.
62  double maxi;
63  // position of maximum
65  // position of correlated point, used only when using projection aware routine
67  // curvature of the correlation peak
69  double maxAngle;
70 };
71 
72 #ifdef HAVE_FFTW
73 
75 static hugin_omp::Lock fftLock;
76 
77 // multiplication with conjugate number in Fourier space
78 template <class Value>
79 Value multiplyConjugate(Value const& v1, Value const& v2)
80 {
81  return v1 * vigra::conj(v2);
82 };
83 
92 template <class SrcImage, class DestImage, class KernelImage>
93 CorrelationResult correlateImageFastFFT(SrcImage & src, DestImage & dest, KernelImage & kernel, vigra::Diff2D kul, vigra::Diff2D klr)
94 {
95  vigra_precondition(kul.x <= 0 && kul.y <= 0,
96  "correlateImageFastFFT(): coordinates of kernel's upper left must be <= 0.");
97  vigra_precondition(klr.x >= 0 && klr.y >= 0,
98  "correlateImageFastFFT(): coordinates of kernel's lower right must be >= 0.");
99 
100  // calculate width and height of the image
101  const int kw = kernel.width();
102  const int kh = kernel.height();
103  const int sw = src.width();
104  const int sh = src.height();
105  vigra_precondition(sw >= kw && sh >= kh, "correlateImageFFT(): kernel larger than image.");
106 
107  CorrelationResult res;
108  if (sw < kw || sh < kh)
109  {
110  return res;
111  };
112 
113  // calculate mean and variance of kernel/template
114  vigra::FindAverageAndVariance<typename KernelImage::PixelType> kMean;
115  vigra::inspectImage(srcImageRange(kernel), kMean);
116  // uniform patch, skip calculation
117  if (kMean.variance(false) == 0)
118  {
119  return res;
120  };
121  // subtract mean from kernel/template
122  vigra::transformImage(srcImageRange(kernel), destImage(kernel), vigra::functor::Arg1() - vigra::functor::Param(kMean.average()));
123  // FFT: we are reusing the fftw_plan
124  vigra::FFTWComplexImage spatial(sw, sh);
125  vigra::FFTWComplexImage fourier(sw, sh);
126  vigra::FFTWComplexImage FKernel(sw, sh);
127  // copy kernel to complex data structure
128  vigra::copyImage(srcImageRange(kernel), destImage(spatial, vigra::FFTWWriteRealAccessor<typename KernelImage::PixelType>()));
129  // now do FFT of kernel
130  fftw_plan fftwplan;
131  {
132  // creation of plan is not thread safe, so use lock
133  hugin_omp::ScopedLock sl(fftLock);
134  fftwplan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(), FFTW_FORWARD, FFTW_ESTIMATE);
135  };
136  fftw_execute(fftwplan);
137  vigra::copyImage(srcImageRange(fourier), destImage(FKernel));
138  // now do FFT of search image, reuse fftw_plan
139  vigra::copyImage(srcImageRange(src), destImage(spatial, vigra::FFTWWriteRealAccessor<typename KernelImage::PixelType>()));
140  fftw_execute(fftwplan);
141  // give now not anymore needed memory free
142  spatial.resize(0, 0);
143 
144  // multiply SrcImage with conjugated kernel in frequency domain
145  vigra::combineTwoImages(srcImageRange(fourier), srcImage(FKernel), destImage(fourier), &multiplyConjugate<vigra::FFTWComplex<typename KernelImage::PixelType> >);
146  // FFT back into spatial domain (inplace)
147  fftw_plan backwardPlan;
148  {
149  // creation of plan is not thread safe, so use lock
150  hugin_omp::ScopedLock sl(fftLock);
151  backwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)fourier.begin(), (fftw_complex *)fourier.begin(), FFTW_BACKWARD, FFTW_ESTIMATE);
152  };
153  fftw_execute(backwardPlan);
154  {
155  // delayed destroy to use only one lock
156  hugin_omp::ScopedLock sl(fftLock);
157  fftw_destroy_plan(fftwplan);
158  fftw_destroy_plan(backwardPlan);
159  };
160 
161  // calculate look up sum tables
162  // use double instead of float!, otherwise there can be truncation errors
163  vigra::DImage s(sw, sh);
164  vigra::DImage s2(sw, sh);
165  double val = src(0, 0);
166  s(0, 0) = val;
167  s2(0, 0) = val * val;
168  // special treatment for first line
169  for (int x = 1; x < sw; ++x)
170  {
171  val = src(x, 0);
172  s(x, 0) = val + s(x - 1, 0);
173  s2(x, 0) = val * val + s2(x - 1, 0);
174  }
175  // special treatment for first column
176  for (int y = 1; y < sh; ++y)
177  {
178  val = src(0, y);
179  s(0, y) = val + s(0, y - 1);
180  s2(0, y) = val * val + s2(0, y - 1);
181  }
182  // final summation
183  for (int y = 1; y < sh; ++y)
184  {
185  for (int x = 1; x < sw; ++x)
186  {
187  val = src(x, y);
188  s(x, y) = val + s(x - 1, y) + s(x, y - 1) - s(x - 1, y - 1);
189  s2(x, y) = val * val + s2(x - 1, y) + s2(x, y - 1) - s2(x - 1, y - 1);
190  };
191  };
192  const int yend = sh - klr.y + kul.y;
193  const int xend = sw - klr.x + kul.x;
194  // calculate constant part
195  const double normFactor = 1.0 / (sw * sh * sqrt(kMean.variance(false)));
196  for (int yr = 0; yr < yend; ++yr)
197  {
198  for (int xr = 0; xr < xend; ++xr)
199  {
200  double value = fourier(xr, yr).re() * normFactor;
201  // do final summation using the lookup tables
202  double sumF = s(xr + kw - 1, yr + kh - 1);
203  double sumF2 = s2(xr + kw - 1, yr + kh - 1);
204  if (xr > 0)
205  {
206  sumF -= s(xr - 1, yr + kh - 1);
207  sumF2 -= s2(xr - 1, yr + kh - 1);
208  };
209  if (yr > 0)
210  {
211  sumF -= s(xr + kw - 1, yr - 1);
212  sumF2 -= s2(xr + kw - 1, yr - 1);
213  };
214  if (xr > 0 && yr > 0)
215  {
216  sumF += s(xr - 1, yr - 1);
217  sumF2 += s2(xr - 1, yr - 1);
218  };
219 
220  double den = sqrt((kw * kh * sumF2 - sumF * sumF));
221  // prevent division through zero
222  if (den != 0)
223  {
224  value /= den;
225  if (value > res.maxi)
226  {
227  res.maxi = value;
228  res.maxpos.x = xr - kul.x;
229  res.maxpos.y = yr - kul.y;
230  }
231  dest(xr - kul.x, yr - kul.y) = vigra::NumericTraits<typename DestImage::value_type>::fromRealPromote(value);
232  };
233  };
234  };
235  return res;
236 };
237 
238 #endif
239 
252 template <class SrcImage, class DestImage, class KernelImage>
254  DestImage & dest,
255  KernelImage & kernel,
256  vigra::Diff2D kul, vigra::Diff2D klr,
257  double threshold = 0.7 )
258 {
259  vigra_precondition(kul.x <= 0 && kul.y <= 0,
260  "convolveImage(): coordinates of "
261  "kernel's upper left must be <= 0.");
262  vigra_precondition(klr.x >= 0 && klr.y >= 0,
263  "convolveImage(): coordinates of "
264  "kernel's lower right must be >= 0.");
265 
266  // use traits to determine SumType as to prevent possible overflow
267  typedef typename
268  vigra::NumericTraits<typename SrcImage::value_type>::RealPromote SumType;
269  typedef typename
270  vigra::NumericTraits<typename KernelImage::value_type>::RealPromote KSumType;
271  typedef
272  vigra::NumericTraits<typename DestImage::value_type> DestTraits;
273 
274  // calculate width and height of the image
275  int w = src.width();
276  int h = src.height();
277  int wk = kernel.width();
278  int hk = kernel.height();
279 
280 // DEBUG_DEBUG("correlate Image srcSize " << (slr - sul).x << "," << (slr - sul).y
281 // << " tmpl size: " << wk << "," << hk);
282 
283  vigra_precondition(w >= wk && h >= hk,
284  "convolveImage(): kernel larger than image.");
285 
286  int ystart = -kul.y;
287  int yend = h-klr.y;
288  int xstart = -kul.x;
289  int xend = w-klr.x;
290 
291  // mean of kernel
292  KSumType kmean=0;
293  for(int y=0; y < hk; y++) {
294  for(int x=0; x < wk; x++) {
295  kmean += kernel(x,y);
296  }
297  }
298  kmean = kmean / (hk*wk);
299 
300  CorrelationResult res;
301 
302 // DEBUG_DEBUG("size: " << w << "," << h << " ystart: " << ystart <<", yend: " << yend);
303  int unifwarned=0;
304  for(int yr=ystart; yr < yend; ++yr)
305  {
306  for(int xr=xstart; xr < xend; ++xr)
307  {
308  if (dest(xr,yr) < threshold) {
309  continue;
310  }
311  // init the sum
312  SumType numerator = 0;
313  SumType div1 = 0;
314  SumType div2 = 0;
315  SumType spixel = 0;
316  KSumType kpixel = 0;
317 
318  // mean of image patch
319  KSumType mean=0;
320  int ym;
321  for(ym=yr+kul.y; ym <= yr+klr.y; ym++) {
322  for(int xm=xr+kul.x; xm <= xr+klr.x; xm++) {
323  mean += src(xm,ym);
324  }
325  }
326  mean = mean / (hk*wk);
327 
328  // perform correlation (inner loop)
329  ym=yr+kul.y;
330  int yk;
331  for(yk=0; yk < hk; yk++) {
332  int xm=xr+kul.x;
333  int xk;
334  for(xk=0; xk < wk; xk++) {
335  spixel = src(xm,ym) - mean;
336  kpixel = kernel(xk,yk) - kmean;
337  numerator += kpixel * spixel;
338  div1 += kpixel * kpixel;
339  div2 += spixel * spixel;
340  xm++;
341  }
342  ym++;
343  }
344  if (div1*div2 == 0) {
345  // This happens when one of the patches is perfectly uniform
346  numerator = 0; // Set correlation to zero since this is uninteresting
347  if (!unifwarned) {
348  DEBUG_DEBUG("Uniform patch(es) during correlation computation");
349  unifwarned=1;
350  }
351  } else
352  numerator = (numerator/sqrt(div1 * div2));
353 
354  if (numerator > res.maxi) {
355  res.maxi = numerator;
356  res.maxpos.x = xr;
357  res.maxpos.y = yr;
358  }
359  dest(xr,yr) = DestTraits::fromRealPromote(numerator);
360  }
361  }
362  return res;
363 }
364 
372 template <class Iterator, class Accessor>
373 CorrelationResult subpixelMaxima(vigra::triple<Iterator, Iterator, Accessor> img,
374  vigra::Diff2D max)
375 {
376  const int interpWidth = 1;
377  CorrelationResult res;
378  vigra_precondition(max.x > interpWidth && max.y > interpWidth,
379  "subpixelMaxima(): coordinates of "
380  "maxima must be > interpWidth, interpWidth.");
381  vigra::Diff2D sz = img.second - img.first;
382  vigra_precondition(sz.x - max.x >= interpWidth && sz.y - max.y >= interpWidth,
383  "subpixelMaxima(): coordinates of "
384  "maxima must interpWidth pixels from the border.");
385  typedef typename Accessor::value_type T;
386 
387  T x[2*interpWidth+1];
388  T zx[2*interpWidth+1];
389  T zy[2*interpWidth+1];
390 
391 #ifdef DEBUG_CORRELATION
392  exportImage(img,vigra::ImageExportInfo("test.tif"));
393 #endif
394 
395  Accessor acc = img.third;
396  Iterator begin=img.first;
397  for (int i=-interpWidth; i<=interpWidth; i++) {
398  // collect first row
399  x[i+interpWidth] = i;
400  zx[i+interpWidth] = acc(begin, max + vigra::Diff2D(i,0));
401  zy[i+interpWidth] = acc(begin, max + vigra::Diff2D(0,i));
402  }
403 
404  double a,b,c;
405  FitPolynom(x, x + 2*interpWidth+1, zx, a,b,c);
406  if (std::isnan(a) || std::isnan(b) || std::isnan(c)) {
407 #ifdef DEBUG_CORRELATION
408  exportImage(img,vigra::ImageExportInfo("test.tif"));
409 #endif
410  DEBUG_DEBUG("Bad polynomial fit results");
411  res.maxpos.x=max.x;
412  res.maxpos.y=max.y;
413  return res;
414  }
415 
416  // calculate extrema of x position by setting
417  // the 1st derivate to zero
418  // 2*c*x + b = 0
419  if (c==0)
420  res.maxpos.x=0;
421  else
422  res.maxpos.x = -b/(2*c);
423 
424  res.curv.x = 2*c;
425 
426  // calculate result at maxima
427  double maxx = c*res.maxpos.x*res.maxpos.x + b*res.maxpos.x + a;
428 
429  FitPolynom(x, x + 2*interpWidth+1, zy, a,b,c);
430  if (std::isnan(a) || std::isnan(b) || std::isnan(c))
431  {
432  DEBUG_DEBUG("Bad polynomial fit results");
433  res.maxpos.x = max.x;
434  res.maxpos.y = max.y;
435  return res;
436  }
437  // calculate extrema of y position
438  if (c==0)
439  res.maxpos.y=0;
440  else
441  res.maxpos.y = -b/(2*c);
442 
443  res.curv.y = 2*c;
444  double maxy = c*res.maxpos.y*res.maxpos.y + b*res.maxpos.y + a;
445 
446  // use mean of both maxima as new interpolation result
447  res.maxi = (maxx + maxy) / 2;
448  DEBUG_DEBUG("value at subpixel maxima (" << res.maxpos.x << " , "
449  <<res.maxpos.y << "): " << maxx << "," << maxy
450  << " mean:" << res.maxi << " coeff: a=" << a
451  << "; b=" << b << "; c=" << c);
452 
453  // test if the point has moved too much ( more than 1.5 pixel).
454  // actually, I should also test the confidence of the fit.
455  if (fabs(res.maxpos.x) > 1 || fabs(res.maxpos.y) > 1) {
456  DEBUG_NOTICE("subpixel Maxima has moved to much, ignoring");
457  res.maxpos.x = max.x;
458  res.maxpos.y = max.y;
459  } else {
460  res.maxpos.x = res.maxpos.x + max.x;
461  res.maxpos.y = res.maxpos.y + max.y;
462  }
463  return res;
464 }
465 
469 {
470 public:
472  : m_phi(phi), m_origin(origin), m_transl(transl)
473  { }
474 
475  bool transformImgCoord(double &destx, double &desty, double srcx, double srcy) const
476  {
477  srcx -= m_origin.x;
478  srcy -= m_origin.y;
479 
480  destx = srcx * cos(m_phi) + srcy * sin(m_phi);
481  desty = srcx * -sin(m_phi) + srcy * cos(m_phi);
482 
483  destx += m_transl.x;
484  desty += m_transl.y;
485  return true;
486  }
487 
488  double m_phi;
491 };
492 
503 template <class IMAGET, class ACCESSORT, class IMAGES, class ACCESSORS>
504 CorrelationResult PointFineTune(const IMAGET & templImg, ACCESSORT access_t,
505  vigra::Diff2D templPos,
506  int templSize,
507  const IMAGES & searchImg, ACCESSORS access_s,
508  vigra::Diff2D searchPos,
509  int sWidth)
510 {
511 // DEBUG_TRACE("templPos: " vigra::<< templPos << " searchPos: " vigra::<< searchPos);
512 
513  // extract patch from template
514 
515  int templWidth = templSize/2;
516  vigra::Diff2D tmplUL(templPos.x - templWidth, templPos.y-templWidth);
517  // lower right iterators "are past the end"
518  vigra::Diff2D tmplLR(templPos.x + templWidth + 1, templPos.y + templWidth + 1);
519  // clip corners to ensure the template is inside the image.
520  vigra::Diff2D tmplImgSize(templImg.size());
521  tmplUL = hugin_utils::simpleClipPoint(tmplUL, vigra::Diff2D(0,0), tmplImgSize);
522  tmplLR = hugin_utils::simpleClipPoint(tmplLR, vigra::Diff2D(0,0), tmplImgSize);
523  vigra::Diff2D tmplSize = tmplLR - tmplUL;
524  DEBUG_DEBUG("template position: " << templPos << " tmplUL: " << tmplUL
525  << " templLR:" << tmplLR << " size:" << tmplSize);
526 
527  // extract patch from search region
528  // make search region bigger, so that interpolation can always be done
529  int swidth = sWidth/2 +(2+templWidth);
530 // DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
531 // Diff2D subjPoint(searchPos);
532  // clip search window
533  if (searchPos.x < 0) searchPos.x = 0;
534  if (searchPos.x > (int) searchImg.width()) searchPos.x = searchImg.width()-1;
535  if (searchPos.y < 0) searchPos.y = 0;
536  if (searchPos.y > (int) searchImg.height()) searchPos.x = searchImg.height()-1;
537 
538  vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
539  // point past the end
540  vigra::Diff2D searchLR(searchPos.x + swidth+1, searchPos.y + swidth+1);
541  // clip search window
542  vigra::Diff2D srcImgSize(searchImg.size());
543  searchUL = hugin_utils::simpleClipPoint(searchUL, vigra::Diff2D(0,0), srcImgSize);
544  searchLR = hugin_utils::simpleClipPoint(searchLR, vigra::Diff2D(0,0), srcImgSize);
545 // DEBUG_DEBUG("search borders: " << searchLR.x << "," << searchLR.y);
546 
547  vigra::Diff2D searchSize = searchLR - searchUL;
548  // create output image
549 
550  // source input
551  vigra::FImage srcImage(searchLR-searchUL);
552  vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
553  searchImg.upperLeft() + searchLR,
554  access_s),
555  destImage(srcImage) );
556 
557  vigra::FImage templateImage(tmplSize);
558  vigra::copyImage(vigra::make_triple(templImg.upperLeft() + tmplUL,
559  templImg.upperLeft() + tmplLR,
560  access_t),
561  destImage(templateImage));
562 #ifdef DEBUG_WRITE_FILES
563  vigra::ImageExportInfo tmpli("hugin_templ.tif");
564  vigra::exportImage(vigra::srcImageRange(templateImage), tmpli);
565 
566  vigra::ImageExportInfo srci("hugin_searchregion.tif");
567  vigra::exportImage(vigra::srcImageRange(srcImage), srci);
568 #endif
569 
570  vigra::FImage dest(searchSize);
571  dest.init(-1);
572  // we could use the multiresolution version as well.
573  // but usually the region is quite small.
574  CorrelationResult res;
575 #ifdef HAVE_FFTW
576  DEBUG_DEBUG("+++++ starting fast correlation");
577  res = correlateImageFastFFT(srcImage, dest, templateImage,
578  tmplUL - templPos, tmplLR - templPos - vigra::Diff2D(1, 1));
579 #elif defined VIGRA_EXT_USE_FAST_CORR
580  DEBUG_DEBUG("+++++ starting fast correlation");
581  res = correlateImageFast(srcImage,
582  dest,
583  templateImage,
584  tmplUL - templPos, tmplLR - templPos - vigra::Diff2D(1, 1),
585  -1);
586 #else
587  DEBUG_DEBUG("+++++ starting normal correlation");
588  res = correlateImage(srcImage.upperLeft(),
589  srcImage.lowerRight(),
590  srcImage.accessor(),
591  dest.upperLeft(),
592  dest.accessor(),
593  templateImage.upperLeft() + templPos,
594  templateImage.accessor(),
595  tmplUL, tmplLR, -1);
596 #endif
597  DEBUG_DEBUG("normal search finished, max:" << res.maxi
598  << " at " << res.maxpos);
599  // do a subpixel maxima estimation
600  // check if the max is inside the pixel boundaries,
601  // and there are enought correlation values for the subpixel
602  // estimation, (2 + templWidth)
603  if (res.maxpos.x > 2 + templWidth && res.maxpos.x < 2*swidth+1-2-templWidth
604  && res.maxpos.y > 2+templWidth && res.maxpos.y < 2*swidth+1-2-templWidth)
605  {
606  // subpixel estimation
608  DEBUG_DEBUG("subpixel position: max:" << res.maxi
609  << " at " << res.maxpos);
610  } else {
611  // not enough values for subpixel estimation.
612  DEBUG_DEBUG("subpixel estimation not done, maxima too close to border");
613  }
614 
615  res.maxpos = res.maxpos + searchUL;
616  return res;
617 }
618 
629 #ifndef HAVE_FFTW
630 template <class IMAGET, class IMAGES>
632  vigra::Diff2D templPos,
633  int templSize,
634  const IMAGES & searchImg,
635  vigra::Diff2D searchPos,
636  int sWidth,
637  double startAngle,
638  double stopAngle,
639  int angleSteps)
640 {
641  DEBUG_TRACE("templPos: " << templPos << " searchPos: " << searchPos);
642 
643  // extract patch from template
644 
645  // make template size user configurable as well?
646  int templWidth = templSize/2;
647  vigra::Diff2D tmplUL(-templWidth, -templWidth);
648  vigra::Diff2D tmplLR(templWidth, templWidth);
649  // clip template
650  if (tmplUL.x + templPos.x < 0) tmplUL.x = -templPos.x;
651  if (tmplUL.y + templPos.y < 0) tmplUL.y = -templPos.y;
652  if (tmplLR.x + templPos.x> templImg.width())
653  tmplLR.x = templImg.width() - templPos.x;
654  if (tmplLR.y + templPos.y > templImg.height())
655  tmplLR.y = templImg.height() - templPos.y;
656  vigra::Diff2D tmplSize = tmplLR - tmplUL + vigra::Diff2D(1,1);
657 
658  // extract patch from search region
659  // make search region bigger, so that interpolation can always be done
660  int swidth = sWidth/2 +(2+templWidth);
661  DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
662 // Diff2D subjPoint(searchPos);
663  // clip search window
664  if (searchPos.x < 0) searchPos.x = 0;
665  if (searchPos.x > (int) searchImg.width()) searchPos.x = searchImg.width()-1;
666  if (searchPos.y < 0) searchPos.y = 0;
667  if (searchPos.y > (int) searchImg.height()) searchPos.x = searchImg.height()-1;
668 
669  vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
670  vigra::Diff2D searchLR(searchPos.x + swidth, searchPos.y + swidth);
671  // clip search window
672  if (searchUL.x < 0) searchUL.x = 0;
673  if (searchUL.x > searchImg.width()) searchUL.x = searchImg.width();
674  if (searchUL.y < 0) searchUL.y = 0;
675  if (searchUL.y > searchImg.height()) searchUL.y = searchImg.height();
676  if (searchLR.x > searchImg.width()) searchLR.x = searchImg.width();
677  if (searchLR.x < 0) searchLR.x = 0;
678  if (searchLR.y > searchImg.height()) searchLR.y = searchImg.height();
679  if (searchLR.y < 0) searchLR.y = 0;
680  DEBUG_DEBUG("search borders: " << searchUL << "," << searchLR << " size: " << searchLR -searchUL);
681 
682 
683 #ifdef DEBUG_WRITE_FILES
684  DEBUG_DEBUG("template area: " << templPos + tmplUL << " -> " << templPos + tmplLR);
685  vigra::exportImage(vigra::make_triple(templImg.upperLeft() + templPos + tmplUL,
686  templImg.upperLeft() + templPos + tmplLR + vigra::Diff2D(1,1),
687  templImg.accessor()),
688  vigra::ImageExportInfo("00_original_template.png"));
689 
690  vigra::exportImage(make_triple(searchImg.upperLeft() + searchUL,
691  searchImg.upperLeft() + searchLR,
692  searchImg.accessor()),
693  vigra::ImageExportInfo("00_searcharea.png"));
694 #endif
695 
696  // rotated template
697  vigra::FImage rotTemplate(tmplSize);
698  vigra::FImage alpha(tmplSize);
699 
700  // correlation output
701  vigra::Diff2D searchSize = searchLR - searchUL;
702  vigra::FImage dest(searchSize);
703  dest.init(1);
704  vigra::FImage bestCorrelation(searchSize);
705 
706  // source input
707  vigra::FImage srcImage(searchLR-searchUL);
708  vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
709  searchImg.upperLeft() + searchLR,
710  vigra::RGBToGrayAccessor<typename IMAGES::value_type>() ),
711  destImage(srcImage) );
712 
713  CorrelationResult bestRes;
714  bestRes.maxi = -1;
715  double bestAngle = 0;
716 
718  // test the image at rotation angles with 30 deg. steps.
719  double step = (stopAngle - startAngle)/(angleSteps-1);
720  double phi=startAngle;
722  for (double i=0; phi <= stopAngle; i++, phi += step) {
723  DEBUG_DEBUG("+++++ Rotating image, phi: " << RAD_TO_DEG(phi));
724  RotateTransform t(phi, hugin_utils::FDiff2D(templWidth, templWidth), templPos);
725  vigra_ext::transformImage(srcImageRange(templImg, vigra::RGBToGrayAccessor<typename IMAGET::value_type>()),
726  destImageRange(rotTemplate),
727  destImage(alpha),
728  vigra::Diff2D(0,0),
729  t,
730  nf,
731  false,
733  &dummy);
734  DEBUG_DEBUG("----- Image rotated");
735 
736  // force a search in at all points.
737  dest.init(-1);
738  DEBUG_DEBUG("+++++ starting correlation");
739 
740  CorrelationResult res;
741  // we could use the multiresolution version as well.
742  // but usually the region is quite small.
743 #ifdef HAVE_FFTW
744  res = correlateImageFastFFT(srcImage, dest, rotTemplate,
745  vigra::Diff2D(-templWidth, -templWidth), vigra::Diff2D(templWidth, templWidth));
746 #elif defined VIGRA_EXT_USE_FAST_CORR
747  res = correlateImageFast(srcImage,
748  dest,
749  rotTemplate,
750  vigra::Diff2D(-templWidth, -templWidth),
751  vigra::Diff2D(templWidth, templWidth),
752  -1);
753 #else
754  res = correlateImage(srcImage.upperLeft(),
755  srcImage.lowerRight(),
756  srcImage.accessor(),
757  dest.upperLeft(),
758  dest.accessor(),
759  rotTemplate.upperLeft() + vigra::Diff2D(templWidth, templWidth),
760  rotTemplate.accessor(),
761  vigra::Diff2D(-templWidth, -templWidth), vigra::Diff2D(templWidth, templWidth),
762  -1);
763 #endif
764  DEBUG_DEBUG("---- correlation finished");
765 
766 #ifdef DEBUG_WRITE_FILES
767  char fname[256];
768  vigra::BImage tmpImg(rotTemplate.size());
769  vigra::copyImage(srcImageRange(rotTemplate),destImage(tmpImg));
770  sprintf(fname, "%3.0f_rotated_template.png", phi/M_PI*180);
771  vigra::exportImage(srcImageRange(tmpImg), vigra::ImageExportInfo(fname));
772 
774  vigra::linearRangeMapping(
775  -1, 1, // src range
776  (unsigned char)0, (unsigned char)255) // dest range
777  );
778  tmpImg.resize(dest.size());
780  sprintf(fname, "%3.0f_corr_result.png", phi/M_PI*180);
781  vigra::exportImage(srcImageRange(tmpImg), vigra::ImageExportInfo(fname));
782 #endif
783 
784 
785  DEBUG_DEBUG("normal search finished, max:" << res.maxi
786  << " at " << res.maxpos << " angle:" << phi/M_PI*180 << "");
787 
788  if (res.maxi > bestRes.maxi) {
789  // remember best correlation.
790  bestCorrelation = dest;
791  bestRes = res;
792  bestAngle = phi;
793  }
794 
795  }
796 
797  DEBUG_DEBUG("rotation search finished, max:" << bestRes.maxi
798  << " at " << bestRes.maxpos << " angle:" << bestAngle/M_PI*180 << "");
799 
800  // do a subpixel maxima estimation
801  // check if the max is inside the pixel boundaries,
802  // and there are enought correlation values for the subpixel
803  // estimation, (2 + templWidth)
804  if (bestRes.maxpos.x > 1 + templWidth && bestRes.maxpos.x < 2*swidth+1-1-templWidth
805  && bestRes.maxpos.y > 1+templWidth && bestRes.maxpos.y < 2*swidth+1-1-templWidth)
806  {
807  // subpixel estimation
808  bestRes = subpixelMaxima(vigra::srcImageRange(bestCorrelation), bestRes.maxpos.toDiff2D());
809  DEBUG_DEBUG("subpixel position: max:" << bestRes.maxi
810  << " at " << bestRes.maxpos << " under angle: " << bestAngle/M_PI*180);
811  } else {
812  // not enough values for subpixel estimation.
813  DEBUG_ERROR("subpixel estimation not done, maxima to close to border");
814  }
815 
816  bestRes.maxpos = bestRes.maxpos + searchUL;
817  bestRes.maxAngle = bestAngle/M_PI*180;
818  return bestRes;
819 }
820 #else
821 // specialed version for FFT based correlation,
822 // optimized correlation calculation to save reuseable results
823 template <class SrcImage, class DestImage, class KernelImage>
824 CorrelationResult correlateImageFastRotationFFT(SrcImage & src, DestImage & dest, KernelImage & unrotatedKernel,
825  vigra::Diff2D templPos, int templWidth, vigra::Diff2D tmplSize,
826  double startAngle, double stopAngle, int angleSteps)
827 {
828  const int sw = src.width();
829  const int sh = src.height();
830  std::vector<CorrelationResult> results(angleSteps);
831  std::vector<DestImage> resultsImg(angleSteps, DestImage(sw, sh));
832 
833  vigra::FFTWComplexImage spatialSearch(sw, sh);
834  vigra::FFTWComplexImage fourierSearch(sw, sh);
835  //forward FFTW plan
836  fftw_plan forwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)spatialSearch.begin(), (fftw_complex *)fourierSearch.begin(), FFTW_FORWARD, FFTW_ESTIMATE);
837  //FFT of search image, we need it for all angles
838  vigra::copyImage(srcImageRange(src), destImage(spatialSearch, vigra::FFTWWriteRealAccessor<vigra::FImage::value_type>()));
839  fftw_execute(forwardPlan);
840  // generate backwardPlan for inplace use
841  fftw_plan backwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)fourierSearch.begin(), (fftw_complex *)fourierSearch.begin(), FFTW_BACKWARD, FFTW_ESTIMATE);
842 
843  // calculate look up sum tables
844  // are used by all angles
845  // use double instead of float!, otherwise there can be truncation errors
846  vigra::DImage s(sw, sh);
847  vigra::DImage s2(sw, sh);
848  double val = src(0, 0);
849  s(0, 0) = val;
850  s2(0, 0) = val * val;
851  // special treatment for first line
852  for (int x = 1; x < sw; ++x)
853  {
854  val = src(x, 0);
855  s(x, 0) = val + s(x - 1, 0);
856  s2(x, 0) = val * val + s2(x - 1, 0);
857  }
858  // special treatment for first column
859  for (int y = 1; y < sh; ++y)
860  {
861  val = src(0, y);
862  s(0, y) = val + s(0, y - 1);
863  s2(0, y) = val * val + s2(0, y - 1);
864  }
865  // final summation
866  for (int y = 1; y < sh; ++y)
867  {
868  for (int x = 1; x < sw; ++x)
869  {
870  val = src(x, y);
871  s(x, y) = val + s(x - 1, y) + s(x, y - 1) - s(x - 1, y - 1);
872  s2(x, y) = val * val + s2(x - 1, y) + s2(x, y - 1) - s2(x - 1, y - 1);
873  };
874  };
875 
876  const double step = (stopAngle - startAngle) / (angleSteps - 1);
877  const int kw = tmplSize.x;
878  const int kh = tmplSize.y;
879  const int yend = sh - 2 * templWidth;
880  const int xend = sw - 2 * templWidth;
881 #pragma omp parallel for
882  for (int i = 0; i < angleSteps; ++i)
883  {
884  vigra::FImage kernel(tmplSize);
885  vigra::FImage alpha(kernel.size());
888 
889  RotateTransform t(startAngle + i * step, hugin_utils::FDiff2D(templWidth, templWidth), templPos);
890  vigra_ext::transformImage(srcImageRange(unrotatedKernel, vigra::RGBToGrayAccessor<typename KernelImage::value_type>()), destImageRange(kernel),
891  destImage(alpha), vigra::Diff2D(0, 0), t, nf, false, vigra_ext::INTERP_CUBIC, &dummy, true);
892 
893  // calculate mean and variance of kernel/template
894  vigra::FindAverageAndVariance<vigra::FImage::PixelType> kMean;
895  vigra::inspectImage(srcImageRange(kernel), kMean);
896  // uniform patch, skip calculation
897  if (kMean.variance(false) == 0)
898  {
899  continue;
900  };
901  // subtract mean from kernel/template
902  vigra::transformImage(srcImageRange(kernel), destImage(kernel), vigra::functor::Arg1() - vigra::functor::Param(kMean.average()));
903 
904  vigra::FFTWComplexImage complexKernel(sw, sh);
905  vigra::FFTWComplexImage fourierKernel(sw, sh);
906  // copy kernel to complex data structure
907  vigra::copyImage(srcImageRange(kernel), destImage(complexKernel, vigra::FFTWWriteRealAccessor<vigra::FImage::value_type>()));
908  fftw_execute_dft(forwardPlan, (fftw_complex*)complexKernel.begin(), (fftw_complex*)fourierKernel.begin());
909 
910  // multiply SrcImage with conjugated kernel in frequency domain
911  vigra::combineTwoImages(srcImageRange(fourierSearch), srcImage(fourierKernel), destImage(fourierKernel), &multiplyConjugate<vigra::FFTWComplex<vigra::FImage::value_type> >);
912  // FFT back into spatial domain (inplace)
913  fftw_execute_dft(backwardPlan, (fftw_complex*)fourierKernel.begin(), (fftw_complex*)fourierKernel.begin());
914 
915  // calculate constant part
916  const double normFactor = 1.0 / (sw * sh * sqrt(kMean.variance(false)));
917  for (int yr = 0; yr < yend; ++yr)
918  {
919  for (int xr = 0; xr < xend; ++xr)
920  {
921  double value = fourierKernel(xr, yr).re() * normFactor;
922  // do final summation using the lookup tables
923  double sumF = s(xr + kw - 1, yr + kh - 1);
924  double sumF2 = s2(xr + kw - 1, yr + kh - 1);
925  if (xr > 0)
926  {
927  sumF -= s(xr - 1, yr + kh - 1);
928  sumF2 -= s2(xr - 1, yr + kh - 1);
929  };
930  if (yr > 0)
931  {
932  sumF -= s(xr + kw - 1, yr - 1);
933  sumF2 -= s2(xr + kw - 1, yr - 1);
934  };
935  if (xr > 0 && yr > 0)
936  {
937  sumF += s(xr - 1, yr - 1);
938  sumF2 += s2(xr - 1, yr - 1);
939  };
940 
941  double den = sqrt((kw * kh * sumF2 - sumF * sumF));
942  // prevent division through zero
943  if (den != 0)
944  {
945  value /= den;
946  if (value > results[i].maxi)
947  {
948  results[i].maxi = value;
949  results[i].maxpos.x = xr + templWidth;
950  results[i].maxpos.y = yr + templWidth;
951  }
952  resultsImg[i](xr + templWidth, yr + templWidth) = vigra::NumericTraits<typename DestImage::value_type>::fromRealPromote(value);
953  };
954  };
955  };
956  };
957  fftw_destroy_plan(forwardPlan);
958  fftw_destroy_plan(backwardPlan);
959  int maxIndex = 0;
960  double maxValue = 0;
961  for (size_t i = 0; i < results.size(); ++i)
962  {
963  if (results[i].maxi > maxValue)
964  {
965  maxIndex = i;
966  maxValue = results[i].maxi;
967  };
968  };
969  CorrelationResult res(results[maxIndex]);
970  res.maxAngle = startAngle + maxIndex * step;
971  dest = resultsImg[maxIndex];
972  return res;
973 };
974 
975 template <class IMAGET, class IMAGES>
976 CorrelationResult PointFineTuneRotSearch(const IMAGET & templImg,
977  vigra::Diff2D templPos,
978  int templSize,
979  const IMAGES & searchImg,
980  vigra::Diff2D searchPos,
981  int sWidth,
982  double startAngle,
983  double stopAngle,
984  int angleSteps)
985 {
986  DEBUG_TRACE("templPos: " << templPos << " searchPos: " << searchPos);
987 
988  // extract patch from template
989  int templWidth = templSize / 2;
990  vigra::Diff2D tmplUL(-templWidth, -templWidth);
991  vigra::Diff2D tmplLR(templWidth, templWidth);
992  // clip template
993  if (tmplUL.x + templPos.x < 0) tmplUL.x = -templPos.x;
994  if (tmplUL.y + templPos.y < 0) tmplUL.y = -templPos.y;
995  if (tmplLR.x + templPos.x> templImg.width())
996  tmplLR.x = templImg.width() - templPos.x;
997  if (tmplLR.y + templPos.y > templImg.height())
998  tmplLR.y = templImg.height() - templPos.y;
999  vigra::Diff2D tmplSize = tmplLR - tmplUL + vigra::Diff2D(1, 1);
1000 
1001  // extract patch from search region
1002  // make search region bigger, so that interpolation can always be done
1003  int swidth = sWidth / 2 + (2 + templWidth);
1004  DEBUG_DEBUG("search window half width/height: " << swidth << "x" << swidth);
1005  // clip search window
1006  if (searchPos.x < 0) searchPos.x = 0;
1007  if (searchPos.x >(int) searchImg.width()) searchPos.x = searchImg.width() - 1;
1008  if (searchPos.y < 0) searchPos.y = 0;
1009  if (searchPos.y >(int) searchImg.height()) searchPos.x = searchImg.height() - 1;
1010 
1011  vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
1012  vigra::Diff2D searchLR(searchPos.x + swidth, searchPos.y + swidth);
1013  // clip search window
1014  if (searchUL.x < 0) searchUL.x = 0;
1015  if (searchUL.x > searchImg.width()) searchUL.x = searchImg.width();
1016  if (searchUL.y < 0) searchUL.y = 0;
1017  if (searchUL.y > searchImg.height()) searchUL.y = searchImg.height();
1018  if (searchLR.x > searchImg.width()) searchLR.x = searchImg.width();
1019  if (searchLR.x < 0) searchLR.x = 0;
1020  if (searchLR.y > searchImg.height()) searchLR.y = searchImg.height();
1021  if (searchLR.y < 0) searchLR.y = 0;
1022  DEBUG_DEBUG("search borders: " << searchUL << "," << searchLR << " size: " << searchLR - searchUL);
1023 
1024  // correlation output
1025  vigra::Diff2D searchSize = searchLR - searchUL;
1026  vigra::FImage dest(searchSize);
1027  dest.init(-1);
1028 
1029  // source input
1030  vigra::FImage srcImage(searchLR - searchUL);
1031  vigra::copyImage(vigra::make_triple(searchImg.upperLeft() + searchUL,
1032  searchImg.upperLeft() + searchLR,
1033  vigra::RGBToGrayAccessor<typename IMAGES::value_type>()),
1034  destImage(srcImage));
1035 
1036  CorrelationResult resCorrelate=correlateImageFastRotationFFT(srcImage, dest, templImg, templPos, templWidth, tmplSize,
1037  startAngle, stopAngle, angleSteps);
1038  CorrelationResult res;
1039 
1040  DEBUG_DEBUG("rotation search finished, max:" << resCorrelate.maxi
1041  << " at " << resCorrelate.maxpos << " angle:" << RAD_TO_DEG(resCorrelate.maxAngle) << "");
1042 
1043  // do a subpixel maxima estimation
1044  // check if the max is inside the pixel boundaries,
1045  // and there are enought correlation values for the subpixel
1046  // estimation, (2 + templWidth)
1047  if (resCorrelate.maxpos.x > 1 + templWidth && resCorrelate.maxpos.x < 2 * swidth + 1 - 1 - templWidth
1048  && resCorrelate.maxpos.y > 1 + templWidth && resCorrelate.maxpos.y < 2 * swidth + 1 - 1 - templWidth)
1049  {
1050  // subpixel estimation
1051  res = subpixelMaxima(vigra::srcImageRange(dest), resCorrelate.maxpos.toDiff2D());
1052 
1053  DEBUG_DEBUG("subpixel position: max:" << res.maxi
1054  << " at " << res.maxpos);
1055  }
1056  else {
1057  // not enough values for subpixel estimation.
1058  res = resCorrelate;
1059  DEBUG_ERROR("subpixel estimation not done, maxima to close to border");
1060  }
1061 
1062  res.maxpos = res.maxpos + searchUL;
1063  res.maxAngle = RAD_TO_DEG(resCorrelate.maxAngle);
1064  return res;
1065 }
1066 #endif
1067 
1077 template <class SrcIterator, class SrcAccessor,
1078  class DestIterator, class DestAccessor,
1079  class KernelIterator, class KernelAccessor>
1080 CorrelationResult correlateImage(SrcIterator sul, SrcIterator slr, SrcAccessor as,
1081  DestIterator dul, DestAccessor ad,
1082  KernelIterator ki, KernelAccessor ak,
1083  vigra::Diff2D kul, vigra::Diff2D klr,
1084  double threshold = 0.7 )
1085 {
1086  CorrelationResult res;
1087  vigra_precondition(kul.x <= 0 && kul.y <= 0,
1088  "convolveImage(): coordinates of "
1089  "kernel's upper left must be <= 0.");
1090  vigra_precondition(klr.x >= 0 && klr.y >= 0,
1091  "convolveImage(): coordinates of "
1092  "kernel's lower right must be >= 0.");
1093 
1094  // use traits to determine SumType as to prevent possible overflow
1095  typedef typename
1096  vigra::NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
1097  typedef typename
1098  vigra::NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
1099  typedef
1100  vigra::NumericTraits<typename DestAccessor::value_type> DestTraits;
1101 
1102  // calculate width and height of the image
1103  int w = slr.x - sul.x;
1104  int h = slr.y - sul.y;
1105  int wk = klr.x - kul.x +1;
1106  int hk = klr.y - kul.y +1;
1107 
1108 // DEBUG_DEBUG("correlate Image srcSize " << (slr - sul).x << "," << (slr - sul).y
1109 // << " tmpl size: " << wk << "," << hk);
1110 
1111  vigra_precondition(w >= wk && h >= hk,
1112  "convolveImage(): kernel larger than image.");
1113 
1114  int x,y;
1115  int ystart = -kul.y;
1116  int yend = h-klr.y;
1117  int xstart = -kul.x;
1118  int xend = w-klr.x;
1119 
1120  // calculate template mean
1121  vigra::FindAverage<typename KernelAccessor::value_type> average;
1122  vigra::inspectImage(ki + kul, ki + klr + vigra::Diff2D(1,1),
1123  ak,
1124  average);
1125  KSumType kmean = average();
1126 
1127  // create y iterators, they iterate over the rows.
1128  DestIterator yd = dul + vigra::Diff2D(xstart, ystart);
1129  SrcIterator ys = sul + vigra::Diff2D(xstart, ystart);
1130 
1131 
1132 // DEBUG_DEBUG("size: " << w << "," << h << " ystart: " << ystart <<", yend: " << yend);
1133  for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y)
1134  {
1135 
1136  // create x iterators, they iterate the coorelation over
1137  // the columns
1138  DestIterator xd(yd);
1139  SrcIterator xs(ys);
1140 
1141  for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x)
1142  {
1143  if (as(xd) < threshold) {
1144  continue;
1145  }
1146 
1147  // init the sum
1148  SumType numerator = 0;
1149  SumType div1 = 0;
1150  SumType div2 = 0;
1151  SumType spixel = 0;
1152  KSumType kpixel = 0;
1153 
1154  // create inner y iterators
1155  // access to the source image
1156  SrcIterator yys = xs + kul;
1157  // access to the kernel image
1158  KernelIterator yk = ki + kul;
1159 
1160  vigra::FindAverage<typename SrcAccessor::value_type> average;
1161  vigra::inspectImage(xs + kul, xs + klr + vigra::Diff2D(1,1), as, average);
1162  double mean = average();
1163 
1164  int xx, yy;
1165  for(yy=0; yy<hk; ++yy, ++yys.y, ++yk.y)
1166  {
1167  // create inner x iterators
1168  SrcIterator xxs = yys;
1169  KernelIterator xk = yk;
1170 
1171  for(xx=0; xx<wk; ++xx, ++xxs.x, ++xk.x)
1172  {
1173  spixel = as(xxs) - mean;
1174  kpixel = ak(xk) - kmean;
1175  numerator += kpixel * spixel;
1176  div1 += kpixel * kpixel;
1177  div2 += spixel * spixel;
1178  }
1179  }
1180  numerator = (numerator/sqrt(div1 * div2));
1181  if (numerator > res.maxi) {
1182  res.maxi = numerator;
1183  res.maxpos.x = x;
1184  res.maxpos.y = y;
1185  }
1186  numerator = numerator;
1187  // store correlation in destination pixel
1188  ad.set(DestTraits::fromRealPromote(numerator), xd);
1189  }
1190  }
1191  return res;
1192 }
1193 
1194 } // namespace
1195 
1196 #endif // VIGRA_EXT_CORRELATION_H
RotateTransform(double phi, hugin_utils::FDiff2D origin, hugin_utils::FDiff2D transl)
Definition: Correlation.h:471
Dummy progress display, without output.
void transformImage(vigra::triple< SrcImageIterator, SrcImageIterator, SrcAccessor > src, vigra::triple< DestImageIterator, DestImageIterator, DestAccessor > dest, std::pair< AlphaImageIterator, AlphaAccessor > alpha, vigra::Diff2D destUL, TRANSFORM &transform, PixelTransform &pixelTransform, bool warparound, Interpolator interpol, AppBase::ProgressDisplay *progress, bool singleThreaded=false)
Transform an image into the panorama.
#define DEBUG_NOTICE(msg)
Definition: utils.h:70
misc math function &amp; classes used by other parts of the program
hugin_utils::FDiff2D corrPos
Definition: Correlation.h:66
CorrelationResult PointFineTune(const IMAGET &templImg, ACCESSORT access_t, vigra::Diff2D templPos, int templSize, const IMAGES &searchImg, ACCESSORS access_s, vigra::Diff2D searchPos, int sWidth)
fine tune a point with normalized cross correlation
Definition: Correlation.h:504
#define DEBUG_TRACE(msg)
Definition: utils.h:67
clockwise rotation around a origin point, and a translation afterwards.
Definition: Correlation.h:468
CorrelationResult subpixelMaxima(vigra::triple< Iterator, Iterator, Accessor > img, vigra::Diff2D max)
find the subpixel maxima by fitting 2nd order polynoms to x and y.
Definition: Correlation.h:373
void FitPolynom(T x, T xend, T y, double &a, double &b, double &c)
fit a second order polynom to a data set
Definition: FitPolynom.h:55
vigra::pair< typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::ImageConstAccessor > srcImage(const ROIImage< Image, Mask > &img)
Definition: ROIImage.h:300
functions to manage ROI&#39;s
CorrelationResult PointFineTuneRotSearch(const IMAGET &templImg, vigra::Diff2D templPos, int templSize, const IMAGES &searchImg, vigra::Diff2D searchPos, int sWidth, double startAngle, double stopAngle, int angleSteps)
fine tune a point with normalized cross correlation, searches x,y and phi (rotation around z) ...
Definition: Correlation.h:631
helper for OpenMP
CorrelationResult correlateImageFast(SrcImage &src, DestImage &dest, KernelImage &kernel, vigra::Diff2D kul, vigra::Diff2D klr, double threshold=0.7)
correlate a template with an image.
Definition: Correlation.h:253
void combineTwoImages(SrcImageIterator1 src1_upperleft, SrcImageIterator1 src1_lowerright, SrcAccessor1 src1_acc, SrcImageIterator2 src2_upperleft, SrcAccessor2 src2_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc, const Functor &func)
Definition: openmp_vigra.h:249
Maximum of correlation, position and value.
Definition: Correlation.h:56
hugin_utils::FDiff2D curv
Definition: Correlation.h:68
IMPEX double h[25][1024]
Definition: emor.cpp:169
vigra::pair< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImage(ROIImage< Image, Alpha > &img)
Definition: ROIImage.h:324
Contains functions to transform whole images.
#define DEBUG_ERROR(msg)
Definition: utils.h:76
vigra::triple< typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::ImageConstAccessor > srcImageRange(const ROIImage< Image, Mask > &img)
helper function for ROIImages
Definition: ROIImage.h:287
std::vector< deghosting::BImagePtr > threshold(const std::vector< deghosting::FImagePtr > &inputImages, const double threshold, const uint16_t flags)
Threshold function used for creating alpha masks for images.
Definition: threshold.h:41
hugin_utils::FDiff2D m_transl
Definition: Correlation.h:490
#define RAD_TO_DEG(x)
Definition: hugin_math.h:45
static T max(T x, T y)
Definition: svm.cpp:65
#define DEBUG_DEBUG(msg)
Definition: utils.h:68
hugin_utils::FDiff2D maxpos
Definition: Correlation.h:64
T simpleClipPoint(const T &point, const T &min, const T &max)
clip a point to fit int [min, max] does not do a mathematical clipping, just sets p...
Definition: hugin_math.h:158
vigra::triple< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImageRange(ROIImage< Image, Alpha > &img)
Definition: ROIImage.h:312
void copyImage(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor src_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc)
Definition: openmp_vigra.h:305
vigra::Diff2D toDiff2D() const
Definition: hugin_math.h:129
bool transformImgCoord(double &destx, double &desty, double srcx, double srcy) const
Definition: Correlation.h:475
#define M_PI
Definition: GaborFilter.cpp:34
CorrelationResult correlateImage(SrcIterator sul, SrcIterator slr, SrcAccessor as, DestIterator dul, DestAccessor ad, KernelIterator ki, KernelAccessor ak, vigra::Diff2D kul, vigra::Diff2D klr, double threshold=0.7)
correlate a template with an image.
Definition: Correlation.h:1080
hugin_utils::FDiff2D m_origin
Definition: Correlation.h:489