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 #include <random>
00040 
00041 //#define DEBUG_WRITE_FILES
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         // apply inverse response/gamma correction, vignetting correction by
00077         // division and possible brightness adjust
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         // apply inverse response/gamma correction, vignetting correction by
00111         // division and possible brightness adjust
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         // apply inverse response/gamma correction, vignetting correction by
00149         // division and possible brightness adjust
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         // apply inverse response/gamma correction, vignetting correction by
00181         // division and possible brightness adjust
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     // Dithering is used to fool the eye into seeing gradients that are finer
00227     // than the precision of the pixel type.
00228     // This prevents the occurence of cleanly-bordered regions in the output where
00229     // the pixel values suddenly change from N to N+1.
00230     // Such regions are especially objectionable in the green channel of 8-bit images.
00231 template <class T>
00232 struct DitherFunctor
00233 {
00234     std::mt19937 Twister;
00235 
00236     typedef T result_type;
00237 
00238     // Dithering is used to fool the eye into seeing gradients that are finer
00239     // than the precision of the pixel type.
00240     // This prevents the occurence of cleanly-bordered regions in the output where
00241     // the pixel values suddenly change from N to N+1.
00242     // Such regions are especially objectionable in the green channel of 8-bit images.
00243     double dither(const double &v) const
00244     {
00245         std::mt19937 &mt = const_cast<std::mt19937 &>(Twister);
00246         double vFraction = v - floor(v);
00247         // Only dither values within a certain range of the rounding cutoff point.
00248         if (vFraction > 0.25 && vFraction <= 0.75) {
00249             // Generate a random number between 0 and 0.5.
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     // dither vector
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     // dither scalar type
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 //        return vigra::NumericTraits<T>::fromRealPromote(vigra_ext::pow(pix, gamma)*maxval);
00303 //        return pow(pix, gamma)*maxval;
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         // calculate mean of flatfield image
00436         vigra::FindAverage<FFT> average;   // init functor
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     // adjust functor
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 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00553 void correctRespVigExpInv(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00554                           vigra::pair<DestIter, DestAccessor> destImg, 
00555                           const std::vector<double> & EMoRCoeff, double maxGreyVal,
00556                           const std::vector<double> & radCoeff, hugin_utils::FDiff2D center, bool division,
00557                           typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00558                           typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00559                           bool dither)
00560 {
00561     typedef typename ImgAccessor::value_type PT;
00562     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00563     typedef typename ImgAccessor::value_type OutPixelType;
00564 
00565     typedef PolySqDistFunctor<4> PolyF;
00566     typedef PassThroughFunctor<RPT> RnF;
00567     typedef LinearTransformFunctor<RPT> LTF;
00568 
00569     PolyF poly(radCoeff);
00570     LTF adjust(a,b);
00571 
00572     // adjust functor
00573 
00574     if (EMoRCoeff.size() == 0) {
00575         RnF Rf;
00576         if (division) {
00577             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00578                                            VigCorrDivFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00579                                            dither);
00580         } else {
00581             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00582                                            VigCorrAddFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00583                                            dither);
00584         }
00585     } else {
00586         // prepare lookup table
00587         // create a camera response lut
00588         std::vector<float> lut10;
00589         vigra_ext::createEMoRLUT(params, lut10);
00590         // create inverse lut
00591         std::vector<float> invlut10(1<<10);
00592         vigra_ext::invertLUT(lut10, invlut10);
00593 
00594         GammaFunctor Rf(gamma, gammaMaxVal);
00595         if (division) {
00596             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00597                                            VigCorrDivFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00598                                            dither);
00599         } else {
00600             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00601                                            VigCorrAddFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00602                                            dither);
00603         }
00604     }
00605 }
00606 
00607 */
00608 
00609 } // namespace
00610 
00611 #endif // VIGRA_EXT_VIGNETTING_CORRECTION_H

Generated on 7 Dec 2016 for Hugintrunk by  doxygen 1.4.7