"Fossies" - the Fresh Open Source Software Archive

Member "pfstools-2.2.0/src/tmo/durand02/tmo_durand02.cpp" (12 Aug 2021, 5517 Bytes) of package /linux/privat/pfstools-2.2.0.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 /**
    2  * @file tmo_bilateral.cpp
    3  * @brief Local tone mapping operator based on bilateral filtering.
    4  * Durand et al. 2002
    5  *
    6  * Fast Bilateral Filtering for the Display of High-Dynamic-Range Images.
    7  * F. Durand and J. Dorsey.
    8  * In ACM Transactions on Graphics, 2002.
    9  *
   10  * 
   11  * This file is a part of PFSTMO package.
   12  * ---------------------------------------------------------------------- 
   13  * Copyright (C) 2003,2004 Grzegorz Krawczyk
   14  * 
   15  *  This program is free software; you can redistribute it and/or modify
   16  *  it under the terms of the GNU General Public License as published by
   17  *  the Free Software Foundation; either version 2 of the License, or
   18  *  (at your option) any later version.
   19  *
   20  *  This program is distributed in the hope that it will be useful,
   21  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   22  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   23  *  GNU General Public License for more details.
   24  *
   25  *  You should have received a copy of the GNU General Public License
   26  *  along with this program; if not, write to the Free Software
   27  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
   28  * ---------------------------------------------------------------------- 
   29  * 
   30  * @author Grzegorz Krawczyk, <krawczyk@mpi-sb.mpg.de>
   31  *
   32  * $Id: tmo_durand02.cpp,v 1.6 2009/02/23 19:09:41 rafm Exp $
   33  */
   34 
   35 #include <config.h>
   36 
   37 #include <iostream>
   38 #include <vector>
   39 #include <algorithm>
   40 #include <math.h>
   41 
   42 #include "pfstmo.h"
   43 
   44 //#undef HAVE_FFTW3F
   45 
   46 #ifdef HAVE_FFTW3F
   47 #include "fastbilateral.h"
   48 #else
   49 #include "bilateral.h"
   50 #endif
   51 
   52 
   53 static void findMaxMinPercentile(pfstmo::Array2D* I, float minPrct, float& minLum, 
   54   float maxPrct, float& maxLum);
   55 
   56 
   57 /*
   58 
   59 From Durand's webpage:
   60 <http://graphics.lcs.mit.edu/~fredo/PUBLI/Siggraph2002/>
   61 
   62 Here is the high-level set of operation that you need to do in order
   63 to perform contrast reduction
   64 
   65 input intensity= 1/61*(R*20+G*40+B)
   66 r=R/(input intensity), g=G/input intensity, B=B/input intensity
   67 log(base)=Bilateral(log(input intensity))
   68 log(detail)=log(input intensity)-log(base)
   69 log (output intensity)=log(base)*compressionfactor+log(detail)
   70 R output = r*exp(log(output intensity)), etc.
   71 
   72 */
   73 
   74 void tmo_durand02(unsigned int width, unsigned int height,
   75   float *nR, float *nG, float *nB,
   76   float sigma_s, float sigma_r, float baseContrast, int downsample,
   77   const bool color_correction,
   78   pfstmo_progress_callback progress_cb ) 
   79 {
   80   pfstmo::Array2D* R = new pfstmo::Array2D(width, height, nR);
   81   pfstmo::Array2D* G = new pfstmo::Array2D(width, height, nG);
   82   pfstmo::Array2D* B = new pfstmo::Array2D(width, height, nB);
   83 
   84   int i;
   85   int w = R->getCols();
   86   int h = R->getRows();
   87   int size = w*h;
   88   pfstmo::Array2D* I = new pfstmo::Array2D(w,h); // intensities
   89   pfstmo::Array2D* BASE = new pfstmo::Array2D(w,h); // base layer
   90   pfstmo::Array2D* DETAIL = new pfstmo::Array2D(w,h); // detail layer
   91 
   92   float min_pos = 1e10f; // minimum positive value (to avoid log(0))
   93   for( i=0 ; i<size ; i++ )
   94   {
   95     (*I)(i) = 1.0f/61.0f * ( 20.0f*(*R)(i) + 40.0f*(*G)(i) + (*B)(i) );
   96     if( unlikely((*I)(i) < min_pos && (*I)(i) > 0) )
   97       min_pos = (*I)(i);
   98   }
   99   
  100   for( i=0 ; i<size ; i++ )
  101   {
  102     float L = (*I)(i);
  103     if( unlikely( L <= 0 ) )
  104       L = min_pos;
  105     
  106     (*R)(i) /= L;
  107     (*G)(i) /= L;
  108     (*B)(i) /= L;
  109 
  110     (*I)(i) = logf( L );
  111   }
  112 
  113 #ifdef HAVE_FFTW3F
  114   fastBilateralFilter( I, BASE, sigma_s, sigma_r, downsample, progress_cb );
  115 #else
  116   bilateralFilter( I, BASE, sigma_s, sigma_r, progress_cb );
  117 #endif
  118 
  119   //!! FIX: find minimum and maximum luminance, but skip 1% of outliers
  120   float maxB,minB;
  121   findMaxMinPercentile(BASE, 0.01f, minB, 0.995f, maxB);
  122 
  123   DEBUG_STR << "Base contrast: " << "maxB=" << maxB << " minB=" << minB
  124             << " c=" << maxB-minB << std::endl;
  125 
  126   float compressionfactor = baseContrast / (maxB-minB);
  127 
  128   DEBUG_STR << "Base contrast (compressed): " << "maxB=" << maxB*compressionfactor
  129             << " minB=" << minB*compressionfactor
  130             << " c=" << (maxB-minB)*compressionfactor << std::endl;
  131 
  132   // Color correction factor
  133   const float k1 = 1.48;
  134   const float k2 = 0.82;
  135   const float s = ( (1 + k1)*pow(compressionfactor,k2) )/( 1 + k1*pow(compressionfactor,k2) );
  136   
  137   for( i=0 ; i<size ; i++ )
  138   {
  139     (*DETAIL)(i) = (*I)(i) - (*BASE)(i);
  140     // maxB will be mapped to log(0)=1
  141     (*I)(i) = ((*BASE)(i)-maxB) * compressionfactor + (*DETAIL)(i);
  142 
  143     if( likely( color_correction ) ) {
  144       (*R)(i) =  powf( (*R)(i), s ) *  expf( (*I)(i) );
  145       (*G)(i) =  powf( (*G)(i), s ) *  expf( (*I)(i) );
  146       (*B)(i) =  powf( (*B)(i), s ) *  expf( (*I)(i) );
  147     } else {
  148       (*R)(i) *= expf( (*I)(i) );
  149       (*G)(i) *= expf( (*I)(i) );
  150       (*B)(i) *= expf( (*I)(i) );
  151     }
  152   }
  153 
  154   delete I;
  155   delete BASE;
  156   delete DETAIL;
  157 
  158   delete B;
  159   delete G;
  160   delete R;
  161 
  162   progress_cb( 100 );
  163 
  164 }
  165 
  166 
  167 
  168 /**
  169  * @brief Find minimum and maximum value skipping the extreems
  170  *
  171  */
  172 static void findMaxMinPercentile(pfstmo::Array2D* I, float minPrct, float& minLum, 
  173   float maxPrct, float& maxLum)
  174 {
  175   int size = I->getRows() * I->getCols();
  176   std::vector<float> vI;
  177 
  178   for( int i=0 ; i<size ; i++ )
  179     if( (*I)(i)!=0.0f )
  180       vI.push_back((*I)(i));
  181       
  182   std::sort(vI.begin(), vI.end());
  183 
  184   minLum = vI.at( int(minPrct*vI.size()) );
  185   maxLum = vI.at( int(maxPrct*vI.size()) );
  186 }