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 
00030 #include <vigra/stdimage.hxx>
00031 #include <vigra/transformimage.hxx>
00032 #include <vigra/inspectimage.hxx>
00033 #include <vigra/combineimages.hxx>
00034 #include <vigra/functorexpression.hxx>
00035 
00036 #include <vigra_ext/utils.h>
00037 
00038 #include <boost/random/mersenne_twister.hpp>
00039 
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     boost::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         boost::mt19937 &mt = const_cast<boost::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     typedef typename DestAccessor::value_type OPT;
00429 
00430     typedef PassThroughFunctor<RPT> RnF;
00431     typedef LinearTransformFunctor<RPT> LTF;
00432     LTF adjust(a,b);
00433 
00434     RFFT mean;
00435     if (division) {
00436         // calculate mean of flatfield image
00437         vigra::FindAverage<FFT> average;   // init functor
00438         vigra::inspectImage(ffImg.first, ffImg.first + (srcImg.second - srcImg.first), ffImg.second, average);
00439         mean = average();
00440     }
00441 
00442     if (gamma == 1.0) {
00443         RnF Rf;
00444         if (division) {
00445             combineTwoImagesDither(srcImg, ffImg, destImg,
00446                                           VigCorrFlatDivFunctor<PT, FFT, RnF, LTF>(mean, Rf, adjust),
00447                                           dither);
00448         } else {
00449             combineTwoImagesDither(srcImg, ffImg, destImg,
00450                                           VigCorrFlatAddFunctor<PT, FFT, RnF, LTF>(Rf, adjust),
00451                                           dither);
00452         }
00453     } else {
00454         GammaFunctor Rf(gamma, gammaMaxVal);
00455         if (division) {
00456             combineTwoImagesDither(srcImg, ffImg, destImg,
00457                                           VigCorrFlatDivFunctor<PT, FFT, GammaFunctor, LTF>(mean, Rf, adjust),
00458                                           dither);
00459         } else {
00460             combineTwoImagesDither(srcImg, ffImg, destImg,
00461                                           VigCorrFlatAddFunctor<PT, FFT, GammaFunctor, LTF>(Rf, adjust),
00462                                           dither);
00463         }
00464     }
00465 }
00466 
00467 
00468 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00469 void radialVigCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00470                          vigra::pair<DestIter, DestAccessor> destImg, double gamma, double gammaMaxVal,
00471                          const std::vector<double> & radCoeff, hugin_utils::FDiff2D center, bool division,
00472                          typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00473                          typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00474                          bool dither)
00475 {
00476     typedef typename ImgAccessor::value_type PT;
00477     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00478     typedef typename ImgAccessor::value_type OutPixelType;
00479 
00480     typedef PolySqDistFunctor<4> PolyF;
00481     typedef PassThroughFunctor<RPT> RnF;
00482     typedef LinearTransformFunctor<RPT> LTF;
00483 
00484     PolyF poly(radCoeff);
00485     LTF adjust(a,b);
00486 
00487     // adjust functor
00488 
00489     if (gamma == 1.0) {
00490         RnF Rf;
00491         if (division) {
00492             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00493                                            VigCorrDivFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00494                                            dither);
00495         } else {
00496             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00497                                            VigCorrAddFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00498                                            dither);
00499         }
00500     } else {
00501         GammaFunctor Rf(gamma, gammaMaxVal);
00502         if (division) {
00503             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00504                                            VigCorrDivFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00505                                            dither);
00506         } else {
00507             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00508                                            VigCorrAddFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00509                                            dither);
00510         }
00511     }
00512 }
00513 
00514 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00515 void applyBrightnessCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00516                                    vigra::pair<DestIter, DestAccessor> destImg,
00517                                    typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00518                                    typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b)
00519 {
00520     typedef typename ImgAccessor::value_type PT;
00521     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00522     
00523     vigra::transformImage(srcImg, destImg, LinearTransformFunctor<RPT>(a,b));
00524 }
00525 
00526 
00527 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00528 void applyGammaAndBrightCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00529                                    vigra::pair<DestIter, DestAccessor> destImg,
00530                                    double gamma, double maxGVal,
00531                                    typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00532                                    typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b)
00533 {
00534     typedef typename ImgAccessor::value_type PT;
00535     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00536     typedef LinearTransformFunctor<RPT> LTF;
00537 
00538     LTF ltf(a,b);
00539     GammaFunctor gf(gamma, maxGVal);
00540     vigra::transformImage(srcImg, destImg, NestFunctor<LTF, GammaFunctor>(ltf, gf));
00541 }
00542 
00543 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00544 void applyGammaCorrection(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00545                           vigra::pair<DestIter, DestAccessor> destImg,
00546                           double gamma, double maxGVal)
00547 {
00548     GammaFunctor gf(gamma, maxGVal);
00549     vigra::transformImage(srcImg, destImg, gf );
00550 }
00551 
00552 
00553 /*
00554 template <class ImgIter, class ImgAccessor, class DestIter, class DestAccessor>
00555 void correctRespVigExpInv(vigra::triple<ImgIter, ImgIter, ImgAccessor> srcImg,
00556                           vigra::pair<DestIter, DestAccessor> destImg, 
00557                           const std::vector<double> & EMoRCoeff, double maxGreyVal,
00558                           const std::vector<double> & radCoeff, hugin_utils::FDiff2D center, bool division,
00559                           typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote a,
00560                           typename vigra::NumericTraits<typename ImgAccessor::value_type>::RealPromote b,
00561                           bool dither)
00562 {
00563     typedef typename ImgAccessor::value_type PT;
00564     typedef typename vigra::NumericTraits<PT>::RealPromote RPT;
00565     typedef typename ImgAccessor::value_type OutPixelType;
00566 
00567     typedef PolySqDistFunctor<4> PolyF;
00568     typedef PassThroughFunctor<RPT> RnF;
00569     typedef LinearTransformFunctor<RPT> LTF;
00570 
00571     PolyF poly(radCoeff);
00572     LTF adjust(a,b);
00573 
00574     // adjust functor
00575 
00576     if (EMoRCoeff.size() == 0) {
00577         RnF Rf;
00578         if (division) {
00579             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y,
00580                                            VigCorrDivFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00581                                            dither);
00582         } else {
00583             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00584                                            VigCorrAddFunctor<PT, RnF, PolyF, LTF>(Rf, poly, adjust),
00585                                            dither);
00586         }
00587     } else {
00588         // prepare lookup table
00589         // create a camera response lut
00590         std::vector<float> lut10;
00591         vigra_ext::createEMoRLUT(params, lut10);
00592         // create inverse lut
00593         std::vector<float> invlut10(1<<10);
00594         vigra_ext::invertLUT(lut10, invlut10);
00595 
00596         GammaFunctor Rf(gamma, gammaMaxVal);
00597         if (division) {
00598             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00599                                            VigCorrDivFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00600                                            dither);
00601         } else {
00602             applyRadialVigCorrectionDither(srcImg, destImg, center.x, center.y, 
00603                                            VigCorrAddFunctor<PT, GammaFunctor, PolyF, LTF>(Rf, poly, adjust),
00604                                            dither);
00605         }
00606     }
00607 }
00608 
00609 */
00610 
00611 } // namespace
00612 
00613 #endif // VIGRA_EXT_VIGNETTING_CORRECTION_H

Generated on 31 Oct 2014 for Hugintrunk by  doxygen 1.4.7