00001
00024 #ifndef VIGRA_EXT_VIGNETTING_CORRECTION_H
00025 #define VIGRA_EXT_VIGNETTING_CORRECTION_H
00026
00027 #include <vector>
00028 #include <functional>
00029 #include "hugin_config.h"
00030
00031 #include <vigra/stdimage.hxx>
00032 #include <vigra/transformimage.hxx>
00033 #include <vigra/inspectimage.hxx>
00034 #include <vigra/combineimages.hxx>
00035 #include <vigra/functorexpression.hxx>
00036
00037 #include <vigra_ext/utils.h>
00038
00039 #include <random>
00040
00041
00042
00043 namespace vigra_ext{
00044
00045 template <class VT1, class VT2, class InvResp, class Adjust>
00046 class VigCorrFlatDivFunctor
00047 {
00048 public:
00051 typedef VT1 first_argument_type;
00052
00055 typedef VT2 second_argument_type;
00056
00059 typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00060
00061 typedef result_type RVT1;
00062 typedef typename vigra::NumericTraits<VT2>::RealPromote RVT2;
00063
00064 VigCorrFlatDivFunctor(RVT2 mean, const InvResp & fr, const Adjust & adj)
00065 : m_InvResp(fr), m_Adjust(adj), m_mean(mean)
00066 { }
00067
00068 InvResp m_InvResp;
00069 Adjust m_Adjust;
00070 RVT2 m_mean;
00071
00074 result_type operator()(first_argument_type const & v1, second_argument_type const & v2) const
00075 {
00076
00077
00078 RVT1 i(m_InvResp(v1));
00079 RVT2 corr((v2)/m_mean);
00080 i /= corr;
00081 return m_Adjust( i );
00082 }
00083 };
00084
00085 template <class VT1, class InvResp, class VigFunc, class Adjust>
00086 class VigCorrDivFunctor
00087 {
00088 public:
00091 typedef VT1 first_argument_type;
00092
00095 typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00096 typedef result_type RealVT1;
00097
00098 VigCorrDivFunctor(const InvResp & fr, const VigFunc & vf, const Adjust & adj)
00099 : m_InvResp(fr), m_VigFunc(vf), m_Adjust(adj)
00100 { }
00101
00102 InvResp m_InvResp;
00103 VigFunc m_VigFunc;
00104 Adjust m_Adjust;
00105
00108 result_type operator()(first_argument_type const & v1, float x, float y) const
00109 {
00110
00111
00112 return m_Adjust(m_InvResp(v1) / m_VigFunc(x,y));
00113 }
00114 };
00115
00116
00117
00118 template <class VT1, class VT2, class InvResp, class Adjust>
00119 class VigCorrFlatAddFunctor
00120 {
00121 public:
00124 typedef VT1 first_argument_type;
00125
00128 typedef VT2 second_argument_type;
00129
00132 typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00133
00134 typedef result_type RealVT1;
00135 typedef typename vigra::NumericTraits<VT2>::RealPromote RVT2;
00136
00137 VigCorrFlatAddFunctor(const InvResp & fr, const Adjust & adj)
00138 : m_InvResp(fr), m_Adjust(adj)
00139 { }
00140
00141 InvResp m_InvResp;
00142 Adjust m_Adjust;
00143
00146 result_type operator()(first_argument_type const & v1, second_argument_type const & v2) const
00147 {
00148
00149
00150 return m_Adjust(m_InvResp(v1) + v2);
00151 }
00152 };
00153
00154 template <class VT1, class InvResp, class VigFunc, class Adjust>
00155 class VigCorrAddFunctor
00156 {
00157 public:
00160 typedef VT1 first_argument_type;
00161
00164 typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00165
00166 typedef result_type RealVT1;
00167
00168 VigCorrAddFunctor(const InvResp & fr, const VigFunc & vf, const Adjust & adj)
00169 : m_InvResp(fr), m_VigFunc(vf), m_Adjust(adj)
00170 { }
00171
00172 InvResp m_InvResp;
00173 VigFunc m_VigFunc;
00174 Adjust m_Adjust;
00175
00178 result_type operator()(first_argument_type const & v1, float x, float y) const
00179 {
00180
00181
00182 return m_Adjust(m_InvResp(v1) + m_VigFunc(x,y));
00183 }
00184 };
00185
00186
00187 template <int NTERMS=4>
00188 struct PolySqDistFunctor
00189 {
00190 double m_coeff[NTERMS];
00191
00192 PolySqDistFunctor(const std::vector<double> & coeff)
00193 {
00194 for (unsigned int i=0; i<NTERMS; i++) m_coeff[i] = coeff[i];
00195 };
00196
00197 double operator()(double x, double y) const
00198 {
00199 double ret = m_coeff[0];
00200 double r2 = x*x + y*y;
00201 double r = r2;
00202 for (unsigned int i = 1; i < NTERMS; i++) {
00203 ret += m_coeff[i] * r;
00204 r *= r2;
00205 }
00206 return ret;
00207 }
00208 };
00209
00210 inline bool isTrueType(vigra::VigraFalseType) {
00211 return false;
00212 }
00213
00214 inline bool isTrueType(vigra::VigraTrueType) {
00215 return true;
00216 }
00217
00218 template<class T>
00219 bool ditheringNeeded(T const &)
00220 {
00221 typedef typename vigra::NumericTraits<T>::isIntegral is_integral;
00222 return isTrueType(is_integral());
00223 }
00224
00226
00227
00228
00229
00230
00231 template <class T>
00232 struct DitherFunctor
00233 {
00234 std::mt19937 Twister;
00235
00236 typedef T result_type;
00237
00238
00239
00240
00241
00242
00243 double dither(const double &v) const
00244 {
00245 std::mt19937 &mt = const_cast<std::mt19937 &>(Twister);
00246 double vFraction = v - floor(v);
00247
00248 if (vFraction > 0.25 && vFraction <= 0.75) {
00249
00250 double random = 0.5 * (double)mt() / UINT_MAX;
00251 if ((vFraction - 0.25) >= random) {
00252 return ceil(v);
00253 } else {
00254 return floor(v);
00255 }
00256 } else {
00257 return v;
00258 }
00259 }
00260
00261
00262 T dither(const T &v, vigra::VigraFalseType) const
00263 {
00264 T ret;
00265 for (size_t i=0; i < v.size(); i++) {
00266 ret[i] = dither(v[i]);
00267 }
00268 return ret;
00269 }
00270
00271
00272 T dither(const T & v, vigra::VigraTrueType) const
00273 {
00274 return dither(v);
00275 }
00276
00277 T operator()(const T &v) const
00278 {
00279 typedef typename vigra::NumericTraits<T>::isScalar is_scalar;
00280
00281 return dither(v, is_scalar());
00282 }
00283 };
00284
00285
00286
00287 struct GammaFunctor
00288 {
00289 double gamma;
00290 double maxval;
00291 GammaFunctor(double g, double m)
00292 : gamma(g), maxval(m)
00293 {}
00294
00295 template <class T>
00296 typename vigra::NumericTraits<T>::RealPromote operator()(T p) const
00297 {
00298 typedef typename vigra::NumericTraits<T>::RealPromote RT;
00299 RT pix(p);
00300 pix /= maxval;
00301 return vigra_ext::pow(pix, gamma)*maxval;
00302
00303
00304 }
00305 };
00306
00309 template <class PT>
00310 struct LinearTransformFunctor
00311 {
00312 typedef PT result_type;
00313
00314 PT m_a;
00315 PT m_b;
00316 LinearTransformFunctor(PT a, PT b)
00317 : m_a(a), m_b(b)
00318 { };
00319
00320 PT operator()(PT p) const
00321 {
00322 return p * m_a + m_b;
00323 }
00324 };
00325
00326
00327 template <class SrcImageIterator, class SrcAccessor,
00328 class DestImageIterator, class DestAccessor, class Functor>
00329 void
00330 applyRadialVigCorrection(vigra::triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00331 vigra::pair<DestImageIterator, DestAccessor> dest,
00332 double cx, double cy,
00333 Functor const & f)
00334 {
00335 applyRadialVigCorrection(src.first, src.second, src.third, dest.first, dest.second,
00336 cx, cy, f);
00337 }
00338
00339 template <class SrcImageIterator, class SrcAccessor,
00340 class DestImageIterator, class DestAccessor, class Functor>
00341 void
00342 applyRadialVigCorrection(SrcImageIterator src_upperleft,
00343 SrcImageIterator src_lowerright, SrcAccessor sa,
00344 DestImageIterator dest_upperleft, DestAccessor da,
00345 double cx, double cy,
00346 Functor const & f)
00347 {
00348 vigra::Diff2D destSize = src_lowerright - src_upperleft;
00349
00350 double sf = 1.0/sqrt(destSize.x/2.0*destSize.x/2.0 + destSize.y/2.0*destSize.y/2.0);
00351 double ynorm = -cy * sf;
00352
00353 for(; src_upperleft.y < src_lowerright.y; ++src_upperleft.y, ++dest_upperleft.y)
00354 {
00355 typename SrcImageIterator::row_iterator s(src_upperleft.rowIterator());
00356 typename SrcImageIterator::row_iterator send(s+ destSize.x);
00357 typename DestImageIterator::row_iterator d(dest_upperleft.rowIterator());
00358 double xnorm = -cx*sf;
00359 for(; s != send; ++s, ++d) {
00360 da.set(f(sa(s), xnorm, ynorm), d);
00361 xnorm +=sf;
00362 }
00363 ynorm += sf;
00364 }
00365 }
00366
00371 template <class SrcImageIterator, class SrcAccessor,
00372 class DestImageIterator, class DestAccessor, class Functor>
00373 void
00374 applyRadialVigCorrectionDither(vigra::triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00375 vigra::pair<DestImageIterator, DestAccessor> dest,
00376 double cx, double cy,
00377 Functor const & f, bool dither)
00378 {
00379 typedef DitherFunctor<typename vigra::NumericTraits<typename SrcAccessor::value_type>::RealPromote> DF;
00380 if (dither) {
00381 DF df;
00382 NestFunctor<DF, Functor> nf(df, f);
00383 applyRadialVigCorrection(src, dest, cx, cy, nf );
00384 } else {
00385 applyRadialVigCorrection(src, dest, cx, cy, f);
00386 }
00387 }
00388
00389
00394 template <class ImgIter, class ImgAccessor,
00395 class FFIter, class FFAccessor,
00396 class DestIter, class DestAccessor,
00397 class Functor>
00398 void
00399 combineTwoImagesDither(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00400 vigra::pair<FFIter, FFAccessor> ffImg,
00401 vigra::pair<DestIter, DestAccessor> destImg,
00402 Functor const & f, bool dither)
00403 {
00404 typedef DitherFunctor<typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote> DF;
00405 if (dither) {
00406 DF df;
00407 vigra::combineTwoImages(srcImg, ffImg, destImg,
00408 NestFunctor<DF, Functor>(df, f) );
00409 } else {
00410 vigra::combineTwoImages(srcImg, ffImg, destImg, f);
00411 }
00412 }
00413
00414
00415 template <class ImgIter, class ImgAccessor, class FFIter, class FFAccessor, class DestIter, class DestAccessor>
00416 void flatfieldVigCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00417 vigra::pair<FFIter, FFAccessor> ffImg,
00418 vigra::pair<DestIter, DestAccessor> destImg,
00419 double gamma, double gammaMaxVal, bool division,
00420 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00421 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00422 bool dither)
00423 {
00424 typedef typename ImgAccessor::value_type PT;
00425 typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00426 typedef typename FFAccessor::value_type FFT;
00427 typedef typename vigra::NumericTraits<FFT>::RealPromote RFFT;
00428
00429 typedef PassThroughFunctor<RPT> RnF;
00430 typedef LinearTransformFunctor<RPT> LTF;
00431 LTF adjust(a,b);
00432
00433 RFFT mean;
00434 if (division) {
00435
00436 vigra::FindAverage<FFT> average;
00437 vigra::inspectImage(ffImg.first, ffImg.first + (srcImg.second - srcImg.first), ffImg.second, average);
00438 mean = average();
00439 }
00440
00441 if (gamma == 1.0) {
00442 RnF Rf;
00443 if (division) {
00444 combineTwoImagesDither(srcImg, ffImg, destImg,
00445 VigCorrFlatDivFunctor<PT, FFT, RnF, LTF>(mean, Rf, adjust),
00446 dither);
00447 } else {
00448 combineTwoImagesDither(srcImg, ffImg, destImg,
00449 VigCorrFlatAddFunctor<PT, FFT, RnF, LTF>(Rf, adjust),
00450 dither);
00451 }
00452 } else {
00453 GammaFunctor Rf(gamma, gammaMaxVal);
00454 if (division) {
00455 combineTwoImagesDither(srcImg, ffImg, destImg,
00456 VigCorrFlatDivFunctor<PT, FFT, GammaFunctor, LTF>(mean, Rf, adjust),
00457 dither);
00458 } else {
00459 combineTwoImagesDither(srcImg, ffImg, destImg,
00460 VigCorrFlatAddFunctor<PT, FFT, GammaFunctor, LTF>(Rf, adjust),
00461 dither);
00462 }
00463 }
00464 }
00465
00466
00467 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00468 void radialVigCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00469 vigra::pair<DestIter, DestAccessor> destImg, double gamma, double gammaMaxVal,
00470 const std::vector<double> & radCoeff, hugin_utils::FDiff2D center, bool division,
00471 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00472 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00473 bool dither)
00474 {
00475 typedef typename ImgAccessor::value_type PT;
00476 typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00477
00478 typedef PolySqDistFunctor<4> PolyF;
00479 typedef PassThroughFunctor<RPT> RnF;
00480 typedef LinearTransformFunctor<RPT> LTF;
00481
00482 PolyF poly(radCoeff);
00483 LTF adjust(a,b);
00484
00485
00486
00487 if (gamma == 1.0) {
00488 RnF Rf;
00489 if (division) {
00490 applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00491 VigCorrDivFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00492 dither);
00493 } else {
00494 applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00495 VigCorrAddFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00496 dither);
00497 }
00498 } else {
00499 GammaFunctor Rf(gamma, gammaMaxVal);
00500 if (division) {
00501 applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00502 VigCorrDivFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00503 dither);
00504 } else {
00505 applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00506 VigCorrAddFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00507 dither);
00508 }
00509 }
00510 }
00511
00512 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00513 void applyBrightnessCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00514 vigra::pair<DestIter, DestAccessor> destImg,
00515 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00516 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b)
00517 {
00518 typedef typename ImgAccessor::value_type PT;
00519 typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00520
00521 vigra::transformImage(srcImg, destImg, LinearTransformFunctor<RPT>(a,b));
00522 }
00523
00524
00525 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00526 void applyGammaAndBrightCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00527 vigra::pair<DestIter, DestAccessor> destImg,
00528 double gamma, double maxGVal,
00529 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00530 typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b)
00531 {
00532 typedef typename ImgAccessor::value_type PT;
00533 typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00534 typedef LinearTransformFunctor<RPT> LTF;
00535
00536 LTF ltf(a,b);
00537 GammaFunctor gf(gamma, maxGVal);
00538 vigra::transformImage(srcImg, destImg, NestFunctor<LTF, GammaFunctor>(ltf, gf));
00539 }
00540
00541 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00542 void applyGammaCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00543 vigra::pair<DestIter, DestAccessor> destImg,
00544 double gamma, double maxGVal)
00545 {
00546 GammaFunctor gf(gamma, maxGVal);
00547 vigra::transformImage(srcImg, destImg, gf );
00548 }
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 }
00610
00611 #endif // VIGRA_EXT_VIGNETTING_CORRECTION_H