24 #ifndef VIGRA_EXT_CORRELATION_H
25 #define VIGRA_EXT_CORRELATION_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>
41 #include "hugin_config.h"
43 #include <vigra/fftw3.hxx>
44 #include <vigra/functorexpression.hxx>
48 #define VIGRA_EXT_USE_FAST_CORR
78 template <
class Value>
79 Value multiplyConjugate(Value
const& v1, Value
const& v2)
81 return v1 * vigra::conj(v2);
92 template <
class SrcImage,
class DestImage,
class KernelImage>
93 CorrelationResult correlateImageFastFFT(SrcImage & src, DestImage & dest, KernelImage & kernel, vigra::Diff2D kul, vigra::Diff2D klr)
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.");
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.");
107 CorrelationResult res;
108 if (sw < kw || sh < kh)
114 vigra::FindAverageAndVariance<typename KernelImage::PixelType> kMean;
117 if (kMean.variance(
false) == 0)
124 vigra::FFTWComplexImage spatial(sw, sh);
125 vigra::FFTWComplexImage fourier(sw, sh);
126 vigra::FFTWComplexImage FKernel(sw, sh);
134 fftwplan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(), FFTW_FORWARD, FFTW_ESTIMATE);
136 fftw_execute(fftwplan);
140 fftw_execute(fftwplan);
142 spatial.resize(0, 0);
147 fftw_plan backwardPlan;
151 backwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)fourier.begin(), (fftw_complex *)fourier.begin(), FFTW_BACKWARD, FFTW_ESTIMATE);
153 fftw_execute(backwardPlan);
157 fftw_destroy_plan(fftwplan);
158 fftw_destroy_plan(backwardPlan);
163 vigra::DImage s(sw, sh);
164 vigra::DImage s2(sw, sh);
165 double val = src(0, 0);
167 s2(0, 0) = val * val;
169 for (
int x = 1; x < sw; ++x)
172 s(x, 0) = val + s(x - 1, 0);
173 s2(x, 0) = val * val + s2(x - 1, 0);
176 for (
int y = 1; y < sh; ++y)
179 s(0, y) = val + s(0, y - 1);
180 s2(0, y) = val * val + s2(0, y - 1);
183 for (
int y = 1; y < sh; ++y)
185 for (
int x = 1; x < sw; ++x)
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);
192 const int yend = sh - klr.y + kul.y;
193 const int xend = sw - klr.x + kul.x;
195 const double normFactor = 1.0 / (sw * sh * sqrt(kMean.variance(
false)));
196 for (
int yr = 0; yr < yend; ++yr)
198 for (
int xr = 0; xr < xend; ++xr)
200 double value = fourier(xr, yr).re() * normFactor;
202 double sumF = s(xr + kw - 1, yr + kh - 1);
203 double sumF2 = s2(xr + kw - 1, yr + kh - 1);
206 sumF -= s(xr - 1, yr + kh - 1);
207 sumF2 -= s2(xr - 1, yr + kh - 1);
211 sumF -= s(xr + kw - 1, yr - 1);
212 sumF2 -= s2(xr + kw - 1, yr - 1);
214 if (xr > 0 && yr > 0)
216 sumF += s(xr - 1, yr - 1);
217 sumF2 += s2(xr - 1, yr - 1);
220 double den = sqrt((kw * kh * sumF2 - sumF * sumF));
225 if (value > res.maxi)
228 res.maxpos.x = xr - kul.x;
229 res.maxpos.y = yr - kul.y;
231 dest(xr - kul.x, yr - kul.y) = vigra::NumericTraits<typename DestImage::value_type>::fromRealPromote(value);
252 template <
class SrcImage,
class DestImage,
class KernelImage>
255 KernelImage & kernel,
256 vigra::Diff2D kul, vigra::Diff2D klr,
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.");
268 vigra::NumericTraits<typename SrcImage::value_type>::RealPromote SumType;
270 vigra::NumericTraits<typename KernelImage::value_type>::RealPromote KSumType;
272 vigra::NumericTraits<typename DestImage::value_type> DestTraits;
276 int h = src.height();
277 int wk = kernel.width();
278 int hk = kernel.height();
283 vigra_precondition(w >= wk && h >= hk,
284 "convolveImage(): kernel larger than image.");
293 for(
int y=0; y < hk; y++) {
294 for(
int x=0; x < wk; x++) {
295 kmean += kernel(x,y);
298 kmean = kmean / (hk*wk);
304 for(
int yr=ystart; yr < yend; ++yr)
306 for(
int xr=xstart; xr < xend; ++xr)
312 SumType numerator = 0;
321 for(ym=yr+kul.y; ym <= yr+klr.y; ym++) {
322 for(
int xm=xr+kul.x; xm <= xr+klr.x; xm++) {
326 mean = mean / (hk*wk);
331 for(yk=0; yk < hk; yk++) {
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;
344 if (div1*div2 == 0) {
348 DEBUG_DEBUG(
"Uniform patch(es) during correlation computation");
352 numerator = (numerator/sqrt(div1 * div2));
354 if (numerator > res.
maxi) {
355 res.
maxi = numerator;
359 dest(xr,yr) = DestTraits::fromRealPromote(numerator);
372 template <
class Iterator,
class Accessor>
376 const int interpWidth = 1;
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;
387 T x[2*interpWidth+1];
388 T zx[2*interpWidth+1];
389 T zy[2*interpWidth+1];
391 #ifdef DEBUG_CORRELATION
392 exportImage(img,vigra::ImageExportInfo(
"test.tif"));
395 Accessor acc = img.third;
396 Iterator begin=img.first;
397 for (
int i=-interpWidth; i<=interpWidth; i++) {
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));
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"));
429 FitPolynom(x, x + 2*interpWidth+1, zy, a,b,c);
430 if (std::isnan(a) || std::isnan(b) || std::isnan(c))
447 res.
maxi = (maxx + maxy) / 2;
449 <<res.
maxpos.
y <<
"): " << maxx <<
"," << maxy
450 <<
" mean:" << res.
maxi <<
" coeff: a=" << a
451 <<
"; b=" << b <<
"; c=" << c);
456 DEBUG_NOTICE(
"subpixel Maxima has moved to much, ignoring");
480 destx = srcx * cos(
m_phi) + srcy * sin(
m_phi);
481 desty = srcx * -sin(
m_phi) + srcy * cos(
m_phi);
503 template <
class IMAGET,
class ACCESSORT,
class IMAGES,
class ACCESSORS>
505 vigra::Diff2D templPos,
507 const IMAGES & searchImg, ACCESSORS access_s,
508 vigra::Diff2D searchPos,
515 int templWidth = templSize/2;
516 vigra::Diff2D tmplUL(templPos.x - templWidth, templPos.y-templWidth);
518 vigra::Diff2D tmplLR(templPos.x + templWidth + 1, templPos.y + templWidth + 1);
520 vigra::Diff2D tmplImgSize(templImg.size());
523 vigra::Diff2D tmplSize = tmplLR - tmplUL;
524 DEBUG_DEBUG(
"template position: " << templPos <<
" tmplUL: " << tmplUL
525 <<
" templLR:" << tmplLR <<
" size:" << tmplSize);
529 int swidth = sWidth/2 +(2+templWidth);
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;
538 vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
540 vigra::Diff2D searchLR(searchPos.x + swidth+1, searchPos.y + swidth+1);
542 vigra::Diff2D srcImgSize(searchImg.size());
547 vigra::Diff2D searchSize = searchLR - searchUL;
551 vigra::FImage
srcImage(searchLR-searchUL);
553 searchImg.upperLeft() + searchLR,
557 vigra::FImage templateImage(tmplSize);
559 templImg.upperLeft() + tmplLR,
562 #ifdef DEBUG_WRITE_FILES
563 vigra::ImageExportInfo tmpli(
"hugin_templ.tif");
566 vigra::ImageExportInfo srci(
"hugin_searchregion.tif");
570 vigra::FImage dest(searchSize);
577 res = correlateImageFastFFT(srcImage, dest, templateImage,
578 tmplUL - templPos, tmplLR - templPos - vigra::Diff2D(1, 1));
579 #elif defined VIGRA_EXT_USE_FAST_CORR
584 tmplUL - templPos, tmplLR - templPos - vigra::Diff2D(1, 1),
589 srcImage.lowerRight(),
593 templateImage.upperLeft() + templPos,
594 templateImage.accessor(),
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)
612 DEBUG_DEBUG(
"subpixel estimation not done, maxima too close to border");
630 template <
class IMAGET,
class IMAGES>
632 vigra::Diff2D templPos,
634 const IMAGES & searchImg,
635 vigra::Diff2D searchPos,
641 DEBUG_TRACE(
"templPos: " << templPos <<
" searchPos: " << searchPos);
646 int templWidth = templSize/2;
647 vigra::Diff2D tmplUL(-templWidth, -templWidth);
648 vigra::Diff2D tmplLR(templWidth, templWidth);
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);
660 int swidth = sWidth/2 +(2+templWidth);
661 DEBUG_DEBUG(
"search window half width/height: " << swidth <<
"x" << swidth);
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;
669 vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
670 vigra::Diff2D searchLR(searchPos.x + swidth, searchPos.y + swidth);
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);
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"));
690 vigra::exportImage(make_triple(searchImg.upperLeft() + searchUL,
691 searchImg.upperLeft() + searchLR,
692 searchImg.accessor()),
693 vigra::ImageExportInfo(
"00_searcharea.png"));
697 vigra::FImage rotTemplate(tmplSize);
698 vigra::FImage alpha(tmplSize);
701 vigra::Diff2D searchSize = searchLR - searchUL;
702 vigra::FImage dest(searchSize);
704 vigra::FImage bestCorrelation(searchSize);
707 vigra::FImage
srcImage(searchLR-searchUL);
709 searchImg.upperLeft() + searchLR,
710 vigra::RGBToGrayAccessor<typename IMAGES::value_type>() ),
715 double bestAngle = 0;
719 double step = (stopAngle - startAngle)/(angleSteps-1);
720 double phi=startAngle;
722 for (
double i=0; phi <= stopAngle; i++, phi += step) {
744 res = correlateImageFastFFT(srcImage, dest, rotTemplate,
745 vigra::Diff2D(-templWidth, -templWidth), vigra::Diff2D(templWidth, templWidth));
746 #elif defined VIGRA_EXT_USE_FAST_CORR
750 vigra::Diff2D(-templWidth, -templWidth),
751 vigra::Diff2D(templWidth, templWidth),
755 srcImage.lowerRight(),
759 rotTemplate.upperLeft() + vigra::Diff2D(templWidth, templWidth),
760 rotTemplate.accessor(),
761 vigra::Diff2D(-templWidth, -templWidth), vigra::Diff2D(templWidth, templWidth),
766 #ifdef DEBUG_WRITE_FILES
768 vigra::BImage tmpImg(rotTemplate.size());
770 sprintf(fname,
"%3.0f_rotated_template.png", phi/
M_PI*180);
771 vigra::exportImage(
srcImageRange(tmpImg), vigra::ImageExportInfo(fname));
774 vigra::linearRangeMapping(
776 (
unsigned char)0, (
unsigned char)255)
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));
786 <<
" at " << res.
maxpos <<
" angle:" << phi/
M_PI*180 <<
"");
788 if (res.
maxi > bestRes.maxi) {
790 bestCorrelation = dest;
797 DEBUG_DEBUG(
"rotation search finished, max:" << bestRes.maxi
798 <<
" at " << bestRes.maxpos <<
" angle:" << bestAngle/
M_PI*180 <<
"");
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)
809 DEBUG_DEBUG(
"subpixel position: max:" << bestRes.maxi
810 <<
" at " << bestRes.maxpos <<
" under angle: " << bestAngle/
M_PI*180);
813 DEBUG_ERROR(
"subpixel estimation not done, maxima to close to border");
816 bestRes.maxpos = bestRes.maxpos + searchUL;
817 bestRes.maxAngle = bestAngle/
M_PI*180;
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)
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));
833 vigra::FFTWComplexImage spatialSearch(sw, sh);
834 vigra::FFTWComplexImage fourierSearch(sw, sh);
836 fftw_plan forwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)spatialSearch.begin(), (fftw_complex *)fourierSearch.begin(), FFTW_FORWARD, FFTW_ESTIMATE);
839 fftw_execute(forwardPlan);
841 fftw_plan backwardPlan = fftw_plan_dft_2d(sh, sw, (fftw_complex *)fourierSearch.begin(), (fftw_complex *)fourierSearch.begin(), FFTW_BACKWARD, FFTW_ESTIMATE);
846 vigra::DImage s(sw, sh);
847 vigra::DImage s2(sw, sh);
848 double val = src(0, 0);
850 s2(0, 0) = val * val;
852 for (
int x = 1; x < sw; ++x)
855 s(x, 0) = val + s(x - 1, 0);
856 s2(x, 0) = val * val + s2(x - 1, 0);
859 for (
int y = 1; y < sh; ++y)
862 s(0, y) = val + s(0, y - 1);
863 s2(0, y) = val * val + s2(0, y - 1);
866 for (
int y = 1; y < sh; ++y)
868 for (
int x = 1; x < sw; ++x)
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);
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)
884 vigra::FImage kernel(tmplSize);
885 vigra::FImage alpha(kernel.size());
889 RotateTransform t(startAngle + i * step,
hugin_utils::FDiff2D(templWidth, templWidth), templPos);
894 vigra::FindAverageAndVariance<vigra::FImage::PixelType> kMean;
897 if (kMean.variance(
false) == 0)
904 vigra::FFTWComplexImage complexKernel(sw, sh);
905 vigra::FFTWComplexImage fourierKernel(sw, sh);
908 fftw_execute_dft(forwardPlan, (fftw_complex*)complexKernel.begin(), (fftw_complex*)fourierKernel.begin());
913 fftw_execute_dft(backwardPlan, (fftw_complex*)fourierKernel.begin(), (fftw_complex*)fourierKernel.begin());
916 const double normFactor = 1.0 / (sw * sh * sqrt(kMean.variance(
false)));
917 for (
int yr = 0; yr < yend; ++yr)
919 for (
int xr = 0; xr < xend; ++xr)
921 double value = fourierKernel(xr, yr).re() * normFactor;
923 double sumF = s(xr + kw - 1, yr + kh - 1);
924 double sumF2 = s2(xr + kw - 1, yr + kh - 1);
927 sumF -= s(xr - 1, yr + kh - 1);
928 sumF2 -= s2(xr - 1, yr + kh - 1);
932 sumF -= s(xr + kw - 1, yr - 1);
933 sumF2 -= s2(xr + kw - 1, yr - 1);
935 if (xr > 0 && yr > 0)
937 sumF += s(xr - 1, yr - 1);
938 sumF2 += s2(xr - 1, yr - 1);
941 double den = sqrt((kw * kh * sumF2 - sumF * sumF));
946 if (value > results[i].maxi)
948 results[i].maxi = value;
949 results[i].maxpos.x = xr + templWidth;
950 results[i].maxpos.y = yr + templWidth;
952 resultsImg[i](xr + templWidth, yr + templWidth) = vigra::NumericTraits<typename DestImage::value_type>::fromRealPromote(value);
957 fftw_destroy_plan(forwardPlan);
958 fftw_destroy_plan(backwardPlan);
961 for (
size_t i = 0; i < results.size(); ++i)
963 if (results[i].maxi > maxValue)
966 maxValue = results[i].maxi;
969 CorrelationResult res(results[maxIndex]);
970 res.maxAngle = startAngle + maxIndex * step;
971 dest = resultsImg[maxIndex];
975 template <
class IMAGET,
class IMAGES>
977 vigra::Diff2D templPos,
979 const IMAGES & searchImg,
980 vigra::Diff2D searchPos,
986 DEBUG_TRACE(
"templPos: " << templPos <<
" searchPos: " << searchPos);
989 int templWidth = templSize / 2;
990 vigra::Diff2D tmplUL(-templWidth, -templWidth);
991 vigra::Diff2D tmplLR(templWidth, templWidth);
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);
1003 int swidth = sWidth / 2 + (2 + templWidth);
1004 DEBUG_DEBUG(
"search window half width/height: " << swidth <<
"x" << swidth);
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;
1011 vigra::Diff2D searchUL(searchPos.x - swidth, searchPos.y - swidth);
1012 vigra::Diff2D searchLR(searchPos.x + swidth, searchPos.y + swidth);
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);
1025 vigra::Diff2D searchSize = searchLR - searchUL;
1026 vigra::FImage dest(searchSize);
1030 vigra::FImage
srcImage(searchLR - searchUL);
1032 searchImg.upperLeft() + searchLR,
1033 vigra::RGBToGrayAccessor<typename IMAGES::value_type>()),
1036 CorrelationResult resCorrelate=correlateImageFastRotationFFT(
srcImage, dest, templImg, templPos, templWidth, tmplSize,
1037 startAngle, stopAngle, angleSteps);
1038 CorrelationResult res;
1040 DEBUG_DEBUG(
"rotation search finished, max:" << resCorrelate.maxi
1041 <<
" at " << resCorrelate.maxpos <<
" angle:" <<
RAD_TO_DEG(resCorrelate.maxAngle) <<
"");
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)
1054 <<
" at " << res.maxpos);
1059 DEBUG_ERROR(
"subpixel estimation not done, maxima to close to border");
1062 res.maxpos = res.maxpos + searchUL;
1063 res.maxAngle =
RAD_TO_DEG(resCorrelate.maxAngle);
1077 template <
class SrcIterator,
class SrcAccessor,
1078 class DestIterator,
class DestAccessor,
1079 class KernelIterator,
class KernelAccessor>
1081 DestIterator dul, DestAccessor ad,
1082 KernelIterator ki, KernelAccessor ak,
1083 vigra::Diff2D kul, vigra::Diff2D klr,
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.");
1096 vigra::NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
1098 vigra::NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
1100 vigra::NumericTraits<typename DestAccessor::value_type> DestTraits;
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;
1111 vigra_precondition(w >= wk && h >= hk,
1112 "convolveImage(): kernel larger than image.");
1115 int ystart = -kul.y;
1117 int xstart = -kul.x;
1121 vigra::FindAverage<typename KernelAccessor::value_type> average;
1122 vigra::inspectImage(ki + kul, ki + klr + vigra::Diff2D(1,1),
1125 KSumType kmean = average();
1128 DestIterator yd = dul + vigra::Diff2D(xstart, ystart);
1129 SrcIterator ys = sul + vigra::Diff2D(xstart, ystart);
1133 for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y)
1138 DestIterator xd(yd);
1141 for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x)
1148 SumType numerator = 0;
1152 KSumType kpixel = 0;
1156 SrcIterator yys = xs + kul;
1158 KernelIterator yk = ki + kul;
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();
1165 for(yy=0; yy<hk; ++yy, ++yys.y, ++yk.y)
1168 SrcIterator xxs = yys;
1169 KernelIterator xk = yk;
1171 for(xx=0; xx<wk; ++xx, ++xxs.x, ++xk.x)
1173 spixel = as(xxs) - mean;
1174 kpixel = ak(xk) - kmean;
1175 numerator += kpixel * spixel;
1176 div1 += kpixel * kpixel;
1177 div2 += spixel * spixel;
1180 numerator = (numerator/sqrt(div1 * div2));
1181 if (numerator > res.
maxi) {
1182 res.
maxi = numerator;
1186 numerator = numerator;
1188 ad.set(DestTraits::fromRealPromote(numerator), xd);
1196 #endif // VIGRA_EXT_CORRELATION_H
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)
misc math function & classes used by other parts of the program
hugin_utils::FDiff2D corrPos
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
CorrelationResult subpixelMaxima(vigra::triple< Iterator, Iterator, Accessor > img, vigra::Diff2D max)
find the subpixel maxima by fitting 2nd order polynoms to x and y.
void FitPolynom(T x, T xend, T y, double &a, double &b, double &c)
fit a second order polynom to a data set
vigra::pair< typename ROIImage< Image, Mask >::image_const_traverser, typename ROIImage< Image, Mask >::ImageConstAccessor > srcImage(const ROIImage< Image, Mask > &img)
functions to manage ROI'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) ...
CorrelationResult correlateImageFast(SrcImage &src, DestImage &dest, KernelImage &kernel, vigra::Diff2D kul, vigra::Diff2D klr, double threshold=0.7)
correlate a template with an image.
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)
Maximum of correlation, position and value.
hugin_utils::FDiff2D curv
vigra::pair< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImage(ROIImage< Image, Alpha > &img)
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
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.
hugin_utils::FDiff2D maxpos
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...
vigra::triple< typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::image_traverser, typename ROIImage< Image, Alpha >::ImageAccessor > destImageRange(ROIImage< Image, Alpha > &img)
void copyImage(SrcImageIterator src_upperleft, SrcImageIterator src_lowerright, SrcAccessor src_acc, DestImageIterator dest_upperleft, DestAccessor dest_acc)
vigra::Diff2D toDiff2D() const
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.