|
TuttleOFX
1
|
00001 #include "NLMDenoiserDefinitions.hpp" 00002 #include "NLMDenoiserPlugin.hpp" 00003 #include "imageUtils/noiseAnalysis.hpp" 00004 00005 #include <tuttle/common/utils/global.hpp> 00006 #include <tuttle/common/math/rectOp.hpp> 00007 #include <tuttle/plugin/ImageGilProcessor.hpp> 00008 #include <tuttle/plugin/IProgress.hpp> 00009 #include <tuttle/plugin/exceptions.hpp> 00010 #include <terry/globals.hpp> 00011 #include <terry/basic_colors.hpp> 00012 #include <terry/channel.hpp> 00013 00014 #include <ofxsImageEffect.h> 00015 #include <ofxsMultiThread.h> 00016 00017 #include <boost/gil/gil_all.hpp> 00018 #include <boost/scoped_ptr.hpp> 00019 00020 #include <cassert> 00021 #include <cmath> 00022 #include <vector> 00023 #include <sstream> 00024 00025 00026 namespace tuttle { 00027 namespace plugin { 00028 namespace nlmDenoiser { 00029 00030 namespace bgil = boost::gil; 00031 00032 template<class View> 00033 NLMDenoiserProcess<View>::NLMDenoiserProcess( NLMDenoiserPlugin & instance ) 00034 : ImageGilProcessor<View>( instance, eImageOrientationIndependant ) 00035 , _plugin( instance ) 00036 { 00037 _paramRedStrength = instance.fetchDoubleParam( kParamRedStrength ); 00038 _paramGreenStrength = instance.fetchDoubleParam( kParamGreenStrength ); 00039 _paramBlueStrength = instance.fetchDoubleParam( kParamBlueStrength ); 00040 00041 _paramRedGrainSize = instance.fetchDoubleParam( kParamRedGrainSize ); 00042 _paramGreenGrainSize = instance.fetchDoubleParam( kParamGreenGrainSize ); 00043 _paramBlueGrainSize = instance.fetchDoubleParam( kParamBlueGrainSize ); 00044 00045 _paramPatchRadius = instance.fetchIntParam( kParamPatchRadius ); 00046 _paramRegionRadius = instance.fetchIntParam( kParamRegionRadius ); 00047 _paramDepth = instance.fetchIntParam( kParamDepth ); 00048 _paramPreBlurring = instance.fetchDoubleParam( kParamPreBlurring ); 00049 00050 _paramOptimized = instance.fetchBooleanParam( kParamOptimization ); 00051 } 00052 00053 template<class View> 00054 NLMDenoiserProcess<View>::~NLMDenoiserProcess() 00055 { 00056 } 00057 00058 template<class View> 00059 void NLMDenoiserProcess<View>::addFrame( const OfxRectI & dBounds, const int dstBitDepth, const int dstComponents, const double time, const int z ) 00060 { 00061 // Fetch main input image 00062 TUTTLE_TLOG( TUTTLE_INFO, "NLMDenoiserProcess<View>::addFrame time:" << time ); 00063 OFX::Image *img = _plugin._clipSrc->fetchImage( time ); 00064 if( !img ) 00065 BOOST_THROW_EXCEPTION( exception::ImageNotReady() ); 00066 _srcImgs.push_back( img ); 00067 OfxRectI bounds = img->getBounds(); 00068 const int w = bounds.x2 - bounds.x1; 00069 const int h = bounds.y2 - bounds.y1; 00070 00071 memcpy( &_upScaledBounds, &bounds, sizeof (OfxRectI ) ); 00072 00073 // Build views 00074 _srcViews.push_back( 00075 bgil::interleaved_view( w, h, (Pixel *) img->getPixelData(), img->getRowDistanceBytes() ) ); 00076 00077 // Make sure bit depths are same 00078 OFX::EBitDepth srcBitDepth = img->getPixelDepth(); 00079 OFX::EPixelComponent srcComponents = img->getPixelComponents(); 00080 00081 // See if they have the same bit depths 00082 if( srcBitDepth != dstBitDepth || srcComponents != dstComponents ) 00083 { 00084 BOOST_THROW_EXCEPTION( exception::BitDepthMismatch() ); 00085 } 00086 } 00087 00088 template<class View> 00089 void NLMDenoiserProcess<View>::setup( const OFX::RenderArguments &args ) 00090 { 00091 const int depth = _paramDepth->getValue(); 00092 00093 _srcViews.clear(); 00094 _srcImgs.clear(); 00095 00096 this->_dst.reset( _plugin._clipDst->fetchImage( args.time ) ); 00097 // Fetch output image 00098 if( !this->_dst.get() ) 00099 BOOST_THROW_EXCEPTION( exception::ImageNotReady() ); 00100 const OfxRectI dBounds = this->_dst->getBounds(); 00101 00102 const OFX::EBitDepth dstBitDepth = this->_dst->getPixelDepth(); 00103 const OFX::EPixelComponent dstComponents = this->_dst->getPixelComponents(); 00104 00105 this->_dstView = bgil::interleaved_view( dBounds.x2 - dBounds.x1, dBounds.y2 - dBounds.y1, 00106 ( Pixel * ) this->_dst->getPixelData(), 00107 this->_dst->getRowDistanceBytes() ); 00108 00109 // Get render frame range 00110 const OfxRangeD clipFullRange = _plugin._clipSrc->getFrameRange(); 00111 OfxRangeD requestedRange; 00112 requestedRange.min = args.time - depth; 00113 requestedRange.max = args.time + depth; 00114 OfxRangeD realRange; 00115 realRange.min = std::max( requestedRange.min, clipFullRange.min ); 00116 realRange.max = std::min( requestedRange.max, clipFullRange.max ); 00117 00118 int i = 0; 00119 addFrame( dBounds, dstBitDepth, dstComponents, (int) ( args.time ), i++ ); 00120 00121 for( OfxTime t = args.time - 1; t >= realRange.min; --t ) 00122 { 00123 TUTTLE_TLOG_VAR2( TUTTLE_INFO, args.time, t ); 00124 addFrame( dBounds, dstBitDepth, dstComponents, t, i++ ); 00125 } 00126 00127 for( OfxTime t = args.time + 1; t <= realRange.max; ++t ) 00128 { 00129 TUTTLE_TLOG_VAR2( TUTTLE_INFO, args.time, t ); 00130 addFrame( dBounds, dstBitDepth, dstComponents, t, i++ ); 00131 } 00132 } 00133 00134 template<class View> 00135 void NLMDenoiserProcess<View>::preProcess() 00136 { 00137 // Initialize progress bar 00138 if( _paramOptimized->getValue() ) 00139 { 00140 const int min_xpi = std::min( (int) _paramRegionRadius->getValue() * 2, (int) _srcViews[0].width() ), 00141 min_ypi = std::min( (int) _paramRegionRadius->getValue() * 2, (int) _srcViews[0].height() ); 00142 std::stringstream msg; 00143 msg << "NL-Means algorithm in progress (automatic bandwidth = " << computeBandwidth() << ")."; 00144 this->progressBegin( _paramDepth->getValue() * min_xpi * min_ypi, msg.str() ); 00145 } 00146 else 00147 { 00148 this->progressBegin( (int) ( ( _srcViews[0].width() * this->_renderArgs.renderScale.x ) * 00149 ( _srcViews[0].height() * this->_renderArgs.renderScale.y ) / 00150 _plugin._clipSrc->getPixelAspectRatio() ), "NL-Means algorithm in progress" ); 00151 } 00152 } 00153 00154 /** 00155 * @brief Function called by rendering thread each time a process must be done. 00156 * @param[in] procWindowRoW Processing window in RoW 00157 */ 00158 template<class View> 00159 void NLMDenoiserProcess<View>::multiThreadProcessImages( const OfxRectI& procWindowRoW ) 00160 { 00161 NlmParams params; 00162 // Change reference point 00163 params.mix[0] = (float) _paramRedStrength->getValue(); 00164 params.mix[1] = (float) _paramGreenStrength->getValue(); 00165 params.mix[2] = (float) _paramBlueStrength->getValue(); 00166 params.mix[2] = 1.0; 00167 00168 params.bws[0] = (float) _paramRedGrainSize->getValue(); 00169 params.bws[1] = (float) _paramGreenGrainSize->getValue(); 00170 params.bws[2] = (float) _paramBlueGrainSize->getValue(); 00171 params.bws[3] = 1.0; 00172 00173 params.patchRadius = _paramPatchRadius->getValue(); 00174 params.regionRadius = _paramRegionRadius->getValue(); 00175 params.preBlurring = (float) _paramPreBlurring->getValue(); 00176 00177 // Destination subview cropped by the procwindow 00178 View subDst = bgil::subimage_view( this->_dstView, 00179 procWindowRoW.x1 - this->_renderArgs.renderWindow.x1, 00180 procWindowRoW.y1 - this->_renderArgs.renderWindow.y1, 00181 procWindowRoW.x2 - procWindowRoW.x1, 00182 procWindowRoW.y2 - procWindowRoW.y1 ); 00183 nlMeans( subDst, procWindowRoW, params ); 00184 } 00185 00186 template<class View> 00187 double NLMDenoiserProcess<View>::computeBandwidth() 00188 { 00189 double bandwidth = 0.0; 00190 // Assume that we work in RGB 00191 const int nc = 3; 00192 double np = ( 2 * _paramPatchRadius->getValue() ) * ( 2 * _paramPatchRadius->getValue() ) * nc; 00193 if( np < 100 ) 00194 bandwidth = ( ( ( ( ( ( 1.1785e-12 * np - 5.1827e-10 ) * np + 9.5946e-08 ) * np - 9.7798e-06 ) * np + 6.0756e-04 ) * np - 0.0248 ) * np + 1.9203 ) * np + 7.9599; 00195 else 00196 bandwidth = ( -7.2611e-04 * np + 1.3213 ) * np + 15.2726; 00197 00198 return bandwidth / _paramDepth->getValue(); 00199 } 00200 00201 /** 00202 * @brief Function called to apply nl-means denoising 00203 * 00204 * @param[out] dst Destination image view 00205 * @param[in] procWindow 00206 * @param[in] params 00207 * 00208 */ 00209 template<class View> 00210 void NLMDenoiserProcess<View>::nlMeans( View &dst, const OfxRectI & procWindow, const NlmParams & params ) 00211 { 00212 using namespace boost::gil; 00213 using namespace terry; 00214 00215 typedef typename rgba32f_view_t::x_iterator WeightIt; 00216 typedef typename View::x_iterator x_iterator; 00217 typedef typename View::locator Locator; 00218 00219 const int w = dst.width(); 00220 const int h = dst.height(); 00221 00222 static const int nc = num_channels<Pixel>::value; 00223 boost::array<float, nc> mix; 00224 00225 for( std::size_t i = 0; i < mix.size(); ++i ) 00226 { 00227 mix[i] = params.mix[i] * channel_traits< Channel >::max_value(); 00228 } 00229 00230 std::vector<View> subSrcViews; 00231 typename std::vector<View>::iterator it; 00232 const int margin = params.regionRadius + params.patchRadius; 00233 00234 // Upscale process window 00235 OfxRectI upscaledProcWindow, tUpscaledProcWindow; 00236 upscaledProcWindow.x1 = procWindow.x1 - margin - 1; 00237 upscaledProcWindow.y1 = procWindow.y1 - margin - 1; 00238 upscaledProcWindow.x2 = procWindow.x2 + margin + 1; 00239 upscaledProcWindow.y2 = procWindow.y2 + margin + 1; 00240 upscaledProcWindow = rectanglesIntersection( upscaledProcWindow, _upScaledBounds ); 00241 00242 // Translate 00243 tUpscaledProcWindow.x1 = upscaledProcWindow.x1 - _upScaledBounds.x1; 00244 tUpscaledProcWindow.y1 = upscaledProcWindow.y1 - _upScaledBounds.y1; 00245 tUpscaledProcWindow.x2 = upscaledProcWindow.x2 - _upScaledBounds.x1; 00246 tUpscaledProcWindow.y2 = upscaledProcWindow.y2 - _upScaledBounds.y1; 00247 00248 // Create up-scaled-proc-windowed subviews sequence 00249 for( it = _srcViews.begin(); it != _srcViews.end(); ++it ) 00250 { 00251 subSrcViews.push_back( subimage_view( *it, tUpscaledProcWindow.x1, tUpscaledProcWindow.y1, 00252 tUpscaledProcWindow.x2 - tUpscaledProcWindow.x1, 00253 tUpscaledProcWindow.y2 - tUpscaledProcWindow.y1 ) ); 00254 } 00255 00256 // Allocate average buffers (assimilate this two buffers as 2D float buffers, not images !) 00257 rgba32f_image_t weight_cumul( w, h ); /// @todo: use memory allocated by host using memorySuite 00258 rgba32f_image_t weight_norm( w, h ); 00259 rgba32f_view_t view_wc( view( weight_cumul ) ); 00260 rgba32f_view_t view_norm( view( weight_norm ) ); 00261 00262 // Initialize averaging buffer to 0 00263 fill_black( view_wc ); 00264 fill_black( view_norm ); 00265 00266 // Bugs bunny is here 00267 OfxRectI nProcWindow; 00268 nProcWindow.x1 = ( procWindow.x1 - upscaledProcWindow.x1 ); 00269 nProcWindow.y1 = ( procWindow.y1 - upscaledProcWindow.y1 ); 00270 nProcWindow.x2 = nProcWindow.x1 + w; 00271 nProcWindow.y2 = nProcWindow.y1 + h; 00272 00273 computeWeights( subSrcViews, nProcWindow, view_wc, view_norm, params ); 00274 00275 if( !_plugin.abort() ) 00276 { 00277 View procView = subimage_view( subSrcViews[0], nProcWindow.x1, nProcWindow.y1, w, h ); 00278 for( int yj = 0; yj < h; ++yj ) 00279 { 00280 WeightIt wcIter = view_wc.row_begin( yj ); 00281 WeightIt wnIter = view_norm.row_begin( yj ); 00282 x_iterator src_it = procView.row_begin( yj ); 00283 x_iterator dst_it = dst.row_begin( yj ); 00284 for( int xj = 0; xj < w; ++xj ) 00285 { 00286 // Final estimate 00287 for( int v = 0; v < nc; ++v ) 00288 { 00289 ( *dst_it )[v] = ( ( ( *wcIter )[v] * mix[v] + 1.0f * ( *src_it )[v] ) / ( ( *wnIter )[v] * mix[v] + 1.0f ) ); 00290 } 00291 00292 // Fill alpha 00293 assign_channel_if_exists_t<Pixel, alpha_t>()( *src_it, *dst_it ); 00294 00295 ++src_it; 00296 ++dst_it; 00297 ++wcIter; 00298 ++wnIter; 00299 } 00300 this->progressForward( w ); 00301 } 00302 } 00303 } 00304 00305 template<class View> 00306 void NLMDenoiserProcess<View>::computeWeights( const std::vector<View> & srcViews, 00307 const OfxRectI & procWindow, 00308 bgil::rgba32f_view_t & view_wc, 00309 bgil::rgba32f_view_t & view_norm, 00310 const NlmParams & params ) 00311 { 00312 typedef typename View::x_iterator sIterator; 00313 typedef typename bgil::rgba32f_view_t::x_iterator wIterator; 00314 typedef typename bgil::channel_type<View>::type dpix_t; 00315 typedef typename View::locator Loc; 00316 typedef typename bgil::rgba32f_view_t::locator WLoc; 00317 00318 const int patchRadius = params.patchRadius; 00319 const int depth = srcViews.size(); 00320 00321 const int wi = srcViews[0].width(); 00322 const int hi = srcViews[0].height(); 00323 00324 // Noise variance estimation 00325 const double nv = imageUtils::noise_variance( srcViews[0] ); 00326 const double sigma = std::sqrt( nv < 0 ? 0 : nv ); 00327 Loc loc1, loc2; 00328 WLoc wcLoc, wnLoc; 00329 00330 // Optimisation based on: AN IMPROVED NON-LOCAL DENOISING ALGORITHM, LNLA 2008 00331 // Define the size of the neighborhood 00332 const int min_xpi = std::min( params.regionRadius, wi / 2 ); 00333 const int min_ypi = std::min( params.regionRadius, hi / 2 ); 00334 00335 static const int nc = boost::mpl::min< boost::mpl::int_<3>, typename bgil::num_channels<Pixel>::type >::type::value; 00336 boost::array<float,nc> bws; 00337 for( std::size_t i = 0; i < bws.size(); ++i ) 00338 { 00339 bws[i] = params.bws[i]; 00340 } 00341 int lbound, hbound; 00342 // [Kervrann] notations 00343 std::vector<double> h1( nc ); 00344 std::vector<double> h2( nc ); 00345 00346 for( int i = 0; i < nc; ++i ) 00347 { 00348 if( bws[i] < 0 ) 00349 bws[i] = (float) computeBandwidth(); 00350 h1[i] = bws[i] * sigma; 00351 h2[i] = 1.0 / ( h1[i] * h1[i] ); 00352 } 00353 00354 double abs_e, eucl_dist, weigth, e; 00355 00356 // For zi (displacment) 00357 for( int zi = 0; zi < depth; ++zi ) 00358 { 00359 // For yi (displacment) 00360 for( int yi = -min_ypi; yi <= min_ypi; ++yi ) 00361 { 00362 // For xi (displacment) 00363 for( int xi = -min_xpi; xi <= min_xpi; ++xi ) 00364 { 00365 // If not 0 displacment 00366 if( xi != 0 || yi != 0 ) 00367 { 00368 const int xl = xi < 0 ? std::abs( xi ) : 0; 00369 const int xh = wi + xi > wi ? wi - xi : wi; 00370 const int yl = yi < 0 ? std::abs( yi ) : 0; 00371 const int yh = hi + yi > hi ? hi - yi : hi; 00372 00373 // For yj 00374 for( int yj = yl; yj < yh; ++yj ) 00375 { 00376 // Initialize patch euclidian distance 00377 eucl_dist = 0.0f; 00378 int j = yj + yi; 00379 // Horizontal & vertical averaging bounds 00380 if( yi >= 0 ) 00381 { 00382 lbound = -patchRadius + yj < 0 ? std::max( -yj, -patchRadius ) : -patchRadius; 00383 hbound = lbound + ( patchRadius * 2 + 1 ); 00384 if( hbound + j > hi ) 00385 hbound = std::min( hi - j, patchRadius ); 00386 } 00387 else 00388 { 00389 hbound = patchRadius + yj > hi ? std::min( patchRadius, ( hi - yj ) ) : patchRadius; 00390 lbound = hbound - ( patchRadius * 2 + 1 ); 00391 if( lbound + j < 0 ) 00392 lbound = 0; 00393 } 00394 00395 // Warmup: initial average accumulation 00396 int xl_bound = std::max( xl - patchRadius, 0 ); 00397 int xr_bound = std::min( wi, xl_bound + patchRadius * 2 ); 00398 loc1 = srcViews[zi].xy_at( xi, j ); 00399 loc2 = srcViews[ 0].xy_at( 0, yj ); 00400 for( int xj = xl_bound; xj < xr_bound; ++xj ) 00401 { 00402 // following "if" is bad but simplify the code 00403 if( ( xj + xi ) >= 0 ) 00404 { 00405 for( int k = lbound; k < hbound; ++k ) 00406 { 00407 for( int v = 0; v < nc; ++v ) 00408 { 00409 e = loc1( xj, k )[v] - loc2( xj, k )[v]; 00410 eucl_dist += e * e; 00411 } 00412 } 00413 } 00414 } 00415 00416 // For xj 00417 for( int xj = xl; xj < xh; ++xj ) 00418 { 00419 int i = xj + xi; 00420 00421 loc1 = srcViews[zi].xy_at( i, j ); 00422 loc2 = srcViews[0].xy_at( xj, yj ); 00423 00424 // 2D centered patch sliding average 00425 if( patchRadius + i < wi && patchRadius + xj < wi ) 00426 { 00427 for( int k = lbound; k < hbound; ++k ) 00428 { 00429 for( int v = 0; v < nc; ++v ) 00430 { 00431 e = loc1( patchRadius, k )[v] - loc2( patchRadius, k )[v]; 00432 eucl_dist += e * e; 00433 } 00434 } 00435 } 00436 00437 if( xj - patchRadius - 1 >= 0 && i - patchRadius - 1 >= 0 ) 00438 { 00439 for( int k = lbound; k < hbound; ++k ) 00440 { 00441 for( int v = 0; v < nc; ++v ) 00442 { 00443 e = loc1( -patchRadius - 1, k )[v] - loc2( -patchRadius - 1, k )[v]; 00444 eucl_dist -= e * e; 00445 } 00446 } 00447 } 00448 00449 // Symetric weigthening will be computed 00450 bool w1Pass = ( zi == 0 && i >= procWindow.x1 && i < procWindow.x2 && 00451 j >= procWindow.y1 && j < procWindow.y2 ); 00452 00453 // Weigthening will be computed 00454 bool w2Pass = ( xj >= procWindow.x1 && xj < procWindow.x2 && 00455 yj >= procWindow.y1 && yj < procWindow.y2 ); 00456 00457 // Weight computation (Modified Bisquare weightening function) 00458 if( w1Pass || w2Pass ) 00459 { 00460 wcLoc = view_wc.xy_at( xj - procWindow.x1, yj - procWindow.y1 ); 00461 wnLoc = view_norm.xy_at( xj - procWindow.x1, yj - procWindow.y1 ); 00462 00463 abs_e = std::abs( eucl_dist ); 00464 for( int v = 0; v < nc; ++v ) 00465 { 00466 if( abs_e <= h1[v] ) 00467 { 00468 weigth = 1.0 - ( abs_e * abs_e * h2[v] ); 00469 // Powerize to 8 00470 weigth *= weigth; 00471 weigth *= weigth; 00472 weigth *= weigth; 00473 00474 // Weight accumulation 00475 if( w1Pass ) 00476 { 00477 wcLoc( xi, yi )[v] += weigth * loc2( 0, 0 )[v]; 00478 wnLoc( xi, yi )[v] += weigth; 00479 } 00480 if( w2Pass ) 00481 { 00482 // Symmetry 00483 wcLoc( 0, 0 )[v] += weigth * loc1( 0, 0 )[v]; 00484 wnLoc( 0, 0 )[v] += weigth; 00485 } 00486 } 00487 } 00488 } 00489 } // End for xj 00490 } // End for yj 00491 } // End if not 0 displacment 00492 if( this->progressForward( 1 ) ) /// @todo what is the step here ?? 00493 return; 00494 } // End for xi (displacment) 00495 } // End for yi (displacment) 00496 } // End for zi (displacment) 00497 } 00498 00499 } 00500 } 00501 }
1.7.5.1 
