contrib/brl/bseg/sdet/sdet_denoise_mrf_bp.cxx
Go to the documentation of this file.
00001 // This is brl/bseg/sdet/sdet_denoise_mrf_bp.cxx
00002 #include "sdet_denoise_mrf_bp.h"
00003 #include <vcl_cstdlib.h>
00004 #include <vul/vul_timer.h>
00005 #include <brip/brip_vil_float_ops.h>
00006 #include <brip/brip_line_generator.h>
00007 
00008 //---------------------------------------------------------------
00009 // Constructors
00010 //----------------------------------------------------------------
00011 
00012 // constructor from a parameter block (the only way)
00013 sdet_denoise_mrf_bp::sdet_denoise_mrf_bp(sdet_denoise_mrf_bp_params& dmp)
00014   : sdet_denoise_mrf_bp_params(dmp), output_valid_(false), use_var_(false)
00015 {}
00016 
00017 // Default Destructor
00018 sdet_denoise_mrf_bp::~sdet_denoise_mrf_bp()
00019 {
00020 }
00021 
00022 void sdet_denoise_mrf_bp::
00023 set_image(vil_image_resource_sptr const& resource)
00024 {
00025   vil_image_view_base_sptr view = resource->get_view();
00026   pyr_in_ = vil_pyramid_image_view<float>(view, pyramid_levels_);
00027 }
00028 
00029 void sdet_denoise_mrf_bp::set_variance(vil_image_resource_sptr const& var_resc)
00030 {
00031   vil_image_view_base_sptr var_view = var_resc->get_view();
00032   pyr_var_ = vil_pyramid_image_view<float>(var_view, pyramid_levels_);
00033   use_var_ = true;
00034 }
00035 
00036 
00037 bool sdet_denoise_mrf_bp::denoise()
00038 {
00039   int lev = static_cast<int>(pyramid_levels_)-1;
00040   vil_image_view<float> depth = pyr_in_(lev);
00041   if (use_var_) {
00042     vil_image_view<float> var = pyr_var_(lev);
00043     mrf_ = new sdet_mrf_bp(depth, var, n_labels_, discontinuity_cost_,
00044                            truncation_cost_, kappa_, lambda_);
00045   }
00046   else
00047     mrf_ = new sdet_mrf_bp(depth, n_labels_, discontinuity_cost_,
00048                            truncation_cost_, kappa_, lambda_);
00049 
00050   vul_timer t;
00051   for (unsigned it = 0; it<n_iter_; ++it) {
00052     mrf_->send_messages_optimized();
00053     vcl_cout << '.'<<vcl_flush;
00054   }
00055   vcl_cout << "completed belief propagation at top level in "
00056            << t.real()/1000.0 << " seconds\n";
00057 
00058   for (--lev ; lev>=0; --lev) {
00059     unsigned pre_lev = static_cast<unsigned>(lev+1);
00060     mrf_ = pyramid_upsample(mrf_, pre_lev);
00061     t.mark();
00062     for (unsigned it = 0; it<n_iter_; ++it) {
00063       mrf_->send_messages_optimized();
00064       vcl_cout << '.'<<vcl_flush;
00065     }
00066     vcl_cout << "completed belief propagation at pyramid level " << lev
00067              << " in " << t.real()/1000.0 << " seconds\n";
00068   }
00069   output_valid_ = true;
00070     out_resc_ = mrf_->belief_image();
00071 
00072   return true;
00073 }
00074 
00075 sdet_mrf_bp_sptr sdet_denoise_mrf_bp::
00076 pyramid_upsample(sdet_mrf_bp_sptr const& in_mrf, unsigned level)
00077 {
00078   if (level==0)
00079     return 0;
00080   unsigned nj = in_mrf->nj(), ni = in_mrf->ni();
00081   if (!ni || !nj) return 0;
00082   //initialize a mrf at the next resolution level
00083   vil_image_view<float> in_view = pyr_in_(level-1);
00084   unsigned njd = in_view.nj(), nid = in_view.ni();
00085    sdet_mrf_bp_sptr out_mrf;
00086   if (use_var_) {
00087     vil_image_view<float> var_view = pyr_var_(level-1);
00088     out_mrf= new sdet_mrf_bp(in_view, var_view, n_labels_,
00089                              discontinuity_cost_,
00090                              truncation_cost_,
00091                              kappa_, lambda_);
00092   }
00093   else
00094     out_mrf = new sdet_mrf_bp(in_view, n_labels_,
00095                               discontinuity_cost_,
00096                               truncation_cost_,
00097                               kappa_, lambda_);
00098   for (unsigned n = 0; n<4; ++n)
00099     for (unsigned j = 0; j<nj; ++j) {
00100       unsigned jn = 2*j;
00101       if (jn+1>=njd) continue;
00102       for (unsigned i = 0; i<ni; ++i)
00103       {
00104         unsigned in = 2*i;
00105         if (in+1>=nid) continue;
00106         vcl_vector<float> msg = in_mrf->prior_message(i, j, n);
00107         out_mrf->set_prior_message(in, jn, n, msg);
00108         out_mrf->set_prior_message(in+1, jn, n, msg);
00109         out_mrf->set_prior_message(in, jn+1, n, msg);
00110         out_mrf->set_prior_message(in+1, jn+1, n, msg);
00111       }
00112     }
00113   return out_mrf;
00114 }