VignettingCorrection.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
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 #ifdef HAVE_CXX11
00040 #include <random>
00041 #else
00042 #include <boost/random/mersenne_twister.hpp>
00043 #endif
00044 
00045 //#define DEBUG_WRITE_FILES
00046 
00047 namespace vigra_ext{
00048 
00049 template <class VT1, class VT2, class InvResp, class Adjust>
00050 class VigCorrFlatDivFunctor
00051 {
00052 public:
00055     typedef VT1 first_argument_type;
00056     
00059     typedef VT2 second_argument_type;
00060      
00063     typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00064 
00065     typedef result_type RVT1;
00066     typedef typename vigra::NumericTraits<VT2>::RealPromote RVT2;
00067 
00068     VigCorrFlatDivFunctor(RVT2 mean, const InvResp & fr, const Adjust & adj)
00069         : m_InvResp(fr), m_Adjust(adj), m_mean(mean)
00070         { }
00071 
00072     InvResp m_InvResp;
00073     Adjust m_Adjust;
00074     RVT2 m_mean;
00075     
00078     result_type operator()(first_argument_type const & v1, second_argument_type const & v2) const
00079     {
00080         // apply inverse response/gamma correction, vignetting correction by
00081         // division and possible brightness adjust
00082         RVT1 i(m_InvResp(v1));
00083         RVT2 corr((v2)/m_mean);
00084         i /= corr;
00085         return m_Adjust( i );
00086     }
00087 };
00088 
00089 template <class VT1, class InvResp, class VigFunc, class Adjust>
00090 class VigCorrDivFunctor
00091 {
00092 public:
00095     typedef VT1 first_argument_type;
00096     
00099     typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00100     typedef result_type RealVT1;
00101 
00102     VigCorrDivFunctor(const InvResp & fr, const VigFunc & vf, const Adjust & adj)
00103         : m_InvResp(fr), m_VigFunc(vf), m_Adjust(adj)
00104         { }
00105 
00106     InvResp m_InvResp;
00107     VigFunc m_VigFunc;
00108     Adjust m_Adjust;
00109     
00112     result_type operator()(first_argument_type const & v1, float x, float y) const
00113     {
00114         // apply inverse response/gamma correction, vignetting correction by
00115         // division and possible brightness adjust
00116         return m_Adjust(m_InvResp(v1) / m_VigFunc(x,y));
00117     }
00118 };
00119 
00120 
00121 
00122 template <class VT1, class VT2, class InvResp, class Adjust>
00123 class VigCorrFlatAddFunctor
00124 {
00125 public:
00128     typedef VT1 first_argument_type;
00129     
00132     typedef VT2 second_argument_type;
00133      
00136     typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00137 
00138     typedef result_type RealVT1;
00139     typedef typename vigra::NumericTraits<VT2>::RealPromote RVT2;
00140 
00141     VigCorrFlatAddFunctor(const InvResp & fr, const Adjust & adj)
00142         : m_InvResp(fr), m_Adjust(adj)
00143         { }
00144 
00145     InvResp m_InvResp;
00146     Adjust m_Adjust;
00147     
00150     result_type operator()(first_argument_type const & v1, second_argument_type const & v2) const
00151     {
00152         // apply inverse response/gamma correction, vignetting correction by
00153         // division and possible brightness adjust
00154         return m_Adjust(m_InvResp(v1) + v2);
00155     }
00156 };
00157 
00158 template <class VT1, class InvResp, class VigFunc, class Adjust>
00159 class VigCorrAddFunctor
00160 {
00161 public:
00164     typedef VT1 first_argument_type;
00165 
00168     typedef typename vigra::NumericTraits<VT1>::RealPromote result_type;
00169 
00170     typedef result_type RealVT1;
00171 
00172     VigCorrAddFunctor(const InvResp & fr, const VigFunc & vf, const Adjust & adj)
00173         : m_InvResp(fr), m_VigFunc(vf), m_Adjust(adj)
00174         { }
00175 
00176     InvResp m_InvResp;
00177     VigFunc m_VigFunc;
00178     Adjust m_Adjust;
00179     
00182     result_type operator()(first_argument_type const & v1, float x, float y) const
00183     {
00184         // apply inverse response/gamma correction, vignetting correction by
00185         // division and possible brightness adjust
00186         return m_Adjust(m_InvResp(v1) + m_VigFunc(x,y));
00187     }
00188 };
00189 
00190 
00191 template <int NTERMS=4>
00192 struct PolySqDistFunctor
00193 {
00194     double m_coeff[NTERMS];
00195 
00196     PolySqDistFunctor(const std::vector<double> & coeff)
00197     { 
00198         for (unsigned int i=0; i<NTERMS; i++) m_coeff[i] = coeff[i];
00199     };
00200 
00201     double operator()(double x, double y) const
00202     {
00203         double ret = m_coeff[0];
00204         double r2 = x*x + y*y;
00205         double r = r2;
00206         for (unsigned int i = 1; i < NTERMS; i++) {
00207             ret += m_coeff[i] * r;
00208             r *= r2;
00209         }
00210         return ret;
00211     }
00212 };
00213 
00214 inline bool isTrueType(vigra::VigraFalseType) {
00215     return false;
00216 }
00217 
00218 inline bool isTrueType(vigra::VigraTrueType) {
00219     return true;
00220 }
00221 
00222 template<class T>
00223 bool ditheringNeeded(T const &)
00224 {
00225     typedef typename vigra::NumericTraits<T>::isIntegral is_integral;
00226     return isTrueType(is_integral());
00227 }
00228 
00230     // Dithering is used to fool the eye into seeing gradients that are finer
00231     // than the precision of the pixel type.
00232     // This prevents the occurence of cleanly-bordered regions in the output where
00233     // the pixel values suddenly change from N to N+1.
00234     // Such regions are especially objectionable in the green channel of 8-bit images.
00235 template <class T>
00236 struct DitherFunctor
00237 {
00238 #ifdef HAVE_CXX11
00239     std::mt19937 Twister;
00240 #else
00241     boost::mt19937 Twister;
00242 #endif
00243 
00244     typedef T result_type;
00245 
00246     // Dithering is used to fool the eye into seeing gradients that are finer
00247     // than the precision of the pixel type.
00248     // This prevents the occurence of cleanly-bordered regions in the output where
00249     // the pixel values suddenly change from N to N+1.
00250     // Such regions are especially objectionable in the green channel of 8-bit images.
00251     double dither(const double &v) const
00252     {
00253 #ifdef HAVE_CXX11
00254         std::mt19937 &mt = const_cast<std::mt19937 &>(Twister);
00255 #else
00256         boost::mt19937 &mt = const_cast<boost::mt19937 &>(Twister);
00257 #endif
00258         double vFraction = v - floor(v);
00259         // Only dither values within a certain range of the rounding cutoff point.
00260         if (vFraction > 0.25 && vFraction <= 0.75) {
00261             // Generate a random number between 0 and 0.5.
00262             double random = 0.5 * (double)mt() / UINT_MAX;
00263             if ((vFraction - 0.25) >= random) {
00264                 return ceil(v);
00265             } else {
00266                 return floor(v);
00267             }
00268         } else {
00269             return v;
00270         }
00271     }
00272 
00273     // dither vector
00274     T dither(const T &v, vigra::VigraFalseType) const
00275     {
00276         T ret;
00277         for (size_t i=0; i < v.size(); i++) {
00278             ret[i] = dither(v[i]);
00279         }
00280         return ret;
00281     }
00282 
00283     // dither scalar type
00284     T dither(const T & v, vigra::VigraTrueType) const
00285     {
00286         return dither(v);
00287     }
00288 
00289     T operator()(const T &v) const
00290     {
00291         typedef typename vigra::NumericTraits<T>::isScalar is_scalar;
00292 
00293         return dither(v, is_scalar());
00294     }
00295 };
00296 
00297 
00298 
00299 struct GammaFunctor
00300 {
00301     double gamma;
00302     double maxval;
00303     GammaFunctor(double g, double m)
00304     : gamma(g), maxval(m)
00305     {}
00306 
00307     template <class T>
00308     typename vigra::NumericTraits<T>::RealPromote operator()(T p) const
00309     {
00310         typedef typename vigra::NumericTraits<T>::RealPromote RT;
00311         RT pix(p);
00312         pix /= maxval;
00313         return vigra_ext::pow(pix, gamma)*maxval;
00314 //        return vigra::NumericTraits<T>::fromRealPromote(vigra_ext::pow(pix, gamma)*maxval);
00315 //        return pow(pix, gamma)*maxval;
00316     }
00317 };
00318 
00321 template <class PT>
00322 struct LinearTransformFunctor
00323 {
00324     typedef PT result_type;
00325 
00326     PT m_a;
00327     PT m_b;
00328     LinearTransformFunctor(PT a, PT b)
00329     : m_a(a), m_b(b)
00330     {  };
00331 
00332     PT operator()(PT p) const
00333     {
00334         return p * m_a + m_b;
00335     }
00336 };
00337 
00338 
00339 template <class SrcImageIterator, class SrcAccessor,
00340           class DestImageIterator, class DestAccessor, class Functor>
00341 void
00342 applyRadialVigCorrection(vigra::triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00343                          vigra::pair<DestImageIterator, DestAccessor> dest,
00344                          double cx, double cy,
00345                          Functor const & f) 
00346 {
00347     applyRadialVigCorrection(src.first, src.second, src.third, dest.first, dest.second,
00348                              cx, cy, f);
00349 }
00350 
00351 template <class SrcImageIterator, class SrcAccessor,
00352           class DestImageIterator, class DestAccessor, class Functor>
00353 void
00354 applyRadialVigCorrection(SrcImageIterator src_upperleft,
00355                          SrcImageIterator src_lowerright, SrcAccessor sa,
00356                          DestImageIterator dest_upperleft, DestAccessor da,
00357                          double cx, double cy,
00358                          Functor const & f) 
00359 {
00360     vigra::Diff2D destSize = src_lowerright - src_upperleft;
00361 
00362     double sf = 1.0/sqrt(destSize.x/2.0*destSize.x/2.0 + destSize.y/2.0*destSize.y/2.0);
00363     double ynorm = -cy * sf;
00364 
00365     for(; src_upperleft.y < src_lowerright.y; ++src_upperleft.y, ++dest_upperleft.y)
00366     {
00367         typename SrcImageIterator::row_iterator s(src_upperleft.rowIterator());
00368         typename SrcImageIterator::row_iterator send(s+ destSize.x);
00369         typename DestImageIterator::row_iterator d(dest_upperleft.rowIterator());
00370         double xnorm = -cx*sf;
00371         for(; s != send; ++s, ++d) {
00372             da.set(f(sa(s), xnorm, ynorm), d);
00373             xnorm +=sf;
00374         }
00375         ynorm += sf;
00376     } 
00377 }
00378 
00383 template <class SrcImageIterator, class SrcAccessor,
00384           class DestImageIterator, class DestAccessor, class Functor>
00385 void
00386 applyRadialVigCorrectionDither(vigra::triple<SrcImageIterator, SrcImageIterator,  SrcAccessor> src,
00387                                vigra::pair<DestImageIterator, DestAccessor> dest,
00388                                double cx, double cy,
00389                                Functor const & f, bool dither) 
00390 {
00391     typedef DitherFunctor<typename vigra::NumericTraits<typename SrcAccessor::value_type>::RealPromote> DF;
00392     if (dither) {
00393         DF df;
00394         NestFunctor<DF, Functor> nf(df, f);
00395         applyRadialVigCorrection(src, dest, cx, cy, nf );
00396     } else {
00397         applyRadialVigCorrection(src, dest, cx, cy, f);
00398     }
00399 }
00400 
00401 
00406 template <class ImgIter, class ImgAccessor, 
00407           class FFIter, class FFAccessor,
00408           class DestIter, class DestAccessor,
00409           class Functor>
00410 void 
00411 combineTwoImagesDither(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00412                        vigra::pair<FFIter, FFAccessor> ffImg,
00413                        vigra::pair<DestIter, DestAccessor> destImg, 
00414                        Functor const & f, bool dither)
00415 {
00416     typedef DitherFunctor<typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote> DF;
00417     if (dither) {
00418         DF df;
00419         vigra::combineTwoImages(srcImg, ffImg, destImg, 
00420                                 NestFunctor<DF, Functor>(df, f) );
00421     } else {
00422         vigra::combineTwoImages(srcImg, ffImg, destImg, f);
00423     }
00424 }
00425 
00426 
00427 template <class ImgIter, class ImgAccessor, class FFIter, class FFAccessor, class DestIter, class DestAccessor>
00428 void flatfieldVigCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00429                             vigra::pair<FFIter, FFAccessor> ffImg,
00430                             vigra::pair<DestIter, DestAccessor> destImg,
00431                             double gamma, double gammaMaxVal, bool division,
00432                             typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00433                             typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00434                             bool dither)
00435 {
00436     typedef typename ImgAccessor::value_type PT;
00437     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00438     typedef typename FFAccessor::value_type FFT;
00439     typedef typename vigra::NumericTraits<FFT>::RealPromote RFFT;
00440     typedef typename DestAccessor::value_type OPT;
00441 
00442     typedef PassThroughFunctor<RPT> RnF;
00443     typedef LinearTransformFunctor<RPT> LTF;
00444     LTF adjust(a,b);
00445 
00446     RFFT mean;
00447     if (division) {
00448         // calculate mean of flatfield image
00449         vigra::FindAverage<FFT> average;   // init functor
00450         vigra::inspectImage(ffImg.first, ffImg.first + (srcImg.second - srcImg.first), ffImg.second, average);
00451         mean = average();
00452     }
00453 
00454     if (gamma == 1.0) {
00455         RnF Rf;
00456         if (division) {
00457             combineTwoImagesDither(srcImg, ffImg, destImg,
00458                                           VigCorrFlatDivFunctor<PT, FFT, RnF, LTF>(mean, Rf, adjust),
00459                                           dither);
00460         } else {
00461             combineTwoImagesDither(srcImg, ffImg, destImg,
00462                                           VigCorrFlatAddFunctor<PT, FFT, RnF, LTF>(Rf, adjust),
00463                                           dither);
00464         }
00465     } else {
00466         GammaFunctor Rf(gamma, gammaMaxVal);
00467         if (division) {
00468             combineTwoImagesDither(srcImg, ffImg, destImg,
00469                                           VigCorrFlatDivFunctor<PT, FFT, GammaFunctor, LTF>(mean, Rf, adjust),
00470                                           dither);
00471         } else {
00472             combineTwoImagesDither(srcImg, ffImg, destImg,
00473                                           VigCorrFlatAddFunctor<PT, FFT, GammaFunctor, LTF>(Rf, adjust),
00474                                           dither);
00475         }
00476     }
00477 }
00478 
00479 
00480 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00481 void radialVigCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00482                          vigra::pair<DestIter, DestAccessor> destImg, double gamma, double gammaMaxVal,
00483                          const std::vector<double> & radCoeff, hugin_utils::FDiff2D center, bool division,
00484                          typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00485                          typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00486                          bool dither)
00487 {
00488     typedef typename ImgAccessor::value_type PT;
00489     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00490     typedef typename ImgAccessor::value_type OutPixelType;
00491 
00492     typedef PolySqDistFunctor<4> PolyF;
00493     typedef PassThroughFunctor<RPT> RnF;
00494     typedef LinearTransformFunctor<RPT> LTF;
00495 
00496     PolyF poly(radCoeff);
00497     LTF adjust(a,b);
00498 
00499     // adjust functor
00500 
00501     if (gamma == 1.0) {
00502         RnF Rf;
00503         if (division) {
00504             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00505                                            VigCorrDivFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00506                                            dither);
00507         } else {
00508             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00509                                            VigCorrAddFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00510                                            dither);
00511         }
00512     } else {
00513         GammaFunctor Rf(gamma, gammaMaxVal);
00514         if (division) {
00515             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00516                                            VigCorrDivFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00517                                            dither);
00518         } else {
00519             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00520                                            VigCorrAddFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00521                                            dither);
00522         }
00523     }
00524 }
00525 
00526 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00527 void applyBrightnessCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00528                                    vigra::pair<DestIter, DestAccessor> destImg,
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     
00535     vigra::transformImage(srcImg, destImg, LinearTransformFunctor<RPT>(a,b));
00536 }
00537 
00538 
00539 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00540 void applyGammaAndBrightCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00541                                    vigra::pair<DestIter, DestAccessor> destImg,
00542                                    double gamma, double maxGVal,
00543                                    typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00544                                    typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b)
00545 {
00546     typedef typename ImgAccessor::value_type PT;
00547     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00548     typedef LinearTransformFunctor<RPT> LTF;
00549 
00550     LTF ltf(a,b);
00551     GammaFunctor gf(gamma, maxGVal);
00552     vigra::transformImage(srcImg, destImg, NestFunctor<LTF, GammaFunctor>(ltf, gf));
00553 }
00554 
00555 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00556 void applyGammaCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00557                           vigra::pair<DestIter, DestAccessor> destImg,
00558                           double gamma, double maxGVal)
00559 {
00560     GammaFunctor gf(gamma, maxGVal);
00561     vigra::transformImage(srcImg, destImg, gf );
00562 }
00563 
00564 
00565 /*
00566 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00567 void correctRespVigExpInv(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00568                           vigra::pair<DestIter, DestAccessor> destImg, 
00569                           const std::vector<double> & EMoRCoeff, double maxGreyVal,
00570                           const std::vector<double> & radCoeff, hugin_utils::FDiff2D center, bool division,
00571                           typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00572                           typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00573                           bool dither)
00574 {
00575     typedef typename ImgAccessor::value_type PT;
00576     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00577     typedef typename ImgAccessor::value_type OutPixelType;
00578 
00579     typedef PolySqDistFunctor<4> PolyF;
00580     typedef PassThroughFunctor<RPT> RnF;
00581     typedef LinearTransformFunctor<RPT> LTF;
00582 
00583     PolyF poly(radCoeff);
00584     LTF adjust(a,b);
00585 
00586     // adjust functor
00587 
00588     if (EMoRCoeff.size() == 0) {
00589         RnF Rf;
00590         if (division) {
00591             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00592                                            VigCorrDivFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00593                                            dither);
00594         } else {
00595             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00596                                            VigCorrAddFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00597                                            dither);
00598         }
00599     } else {
00600         // prepare lookup table
00601         // create a camera response lut
00602         std::vector<float> lut10;
00603         vigra_ext::createEMoRLUT(params, lut10);
00604         // create inverse lut
00605         std::vector<float> invlut10(1<<10);
00606         vigra_ext::invertLUT(lut10, invlut10);
00607 
00608         GammaFunctor Rf(gamma, gammaMaxVal);
00609         if (division) {
00610             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00611                                            VigCorrDivFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00612                                            dither);
00613         } else {
00614             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00615                                            VigCorrAddFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00616                                            dither);
00617         }
00618     }
00619 }
00620 
00621 */
00622 
00623 } // namespace
00624 
00625 #endif // VIGRA_EXT_VIGNETTING_CORRECTION_H

Generated on 31 Aug 2015 for Hugintrunk by  doxygen 1.4.7