1 // This is brl/bpro/core/brad_pro/processes/brad_estimate_empty_process.cxx
2 #include <iostream>
3 #include <algorithm>
4 #include <bprb/bprb_func_process.h>
5 #include <bpro/core/bbas_pro/bbas_1d_array_float.h>
6 #include <brad/brad_phongs_model_est.h>
7 #ifdef _MSC_VER
8 #  include "vcl_msvc_warnings.h"
9 #endif
10 #include <vnl/algo/vnl_levenberg_marquardt.h>
11 #include "vnl/vnl_math.h"
12 //:
13 // \file
14 namespace brad_estimate_empty_process_globals
15 {
16     constexpr unsigned n_inputs_ = 4;
17     constexpr unsigned n_outputs_ = 1;
18 }
19 
20 
21 //: Constructor
brad_estimate_empty_process_cons(bprb_func_process & pro)22 bool brad_estimate_empty_process_cons(bprb_func_process& pro)
23 {
24     using namespace brad_estimate_empty_process_globals;
25 
26     std::vector<std::string> input_types_(n_inputs_);
27     input_types_[0] = "bbas_1d_array_float_sptr";
28     input_types_[1] = "bbas_1d_array_float_sptr";
29     input_types_[2] = "bbas_1d_array_float_sptr";
30     input_types_[3] = "bbas_1d_array_float_sptr";
31 
32     std::vector<std::string>  output_types_(n_outputs_);
33 
34     output_types_[0] = "float";
35 
36     return pro.set_input_types(input_types_)
37         && pro.set_output_types(output_types_);
38 }
39 
40 
41 //: Execute the process
brad_estimate_empty_process(bprb_func_process & pro)42 bool brad_estimate_empty_process(bprb_func_process& pro)
43 {
44     // Sanity check
45     if (pro.n_inputs()< 4) {
46         std::cout << "brip_extrema_process: The input number should be 6" << std::endl;
47         return false;
48     }
49 
50     // get the inputs
51     unsigned i=0;
52     bbas_1d_array_float_sptr intensities = pro.get_input<bbas_1d_array_float_sptr>(i++);
53     bbas_1d_array_float_sptr visibilities = pro.get_input<bbas_1d_array_float_sptr>(i++);
54     bbas_1d_array_float_sptr camera_elev_array = pro.get_input<bbas_1d_array_float_sptr>(i++);
55     bbas_1d_array_float_sptr camera_azim_array = pro.get_input<bbas_1d_array_float_sptr>(i++);
56 
57     unsigned num_samples=intensities->data_array.size();
58     vnl_vector<double> Iobs(num_samples);
59     vnl_vector<double> vis(num_samples);
60     vnl_vector<double> camera_elev(num_samples);
61     vnl_vector<double> camera_azim(num_samples);
62 
63     float mean_intensities = 0.0f ;
64     float sum_weights = 0.0f ;
65     for (unsigned i=0;i<num_samples;i++)
66     {
67         camera_elev[i]    =camera_elev_array->data_array[i];
68         camera_azim[i]    =camera_azim_array->data_array[i];
69 
70         Iobs[i]        =intensities->data_array[i];
71         vis[i]=visibilities->data_array[i];
72         if (Iobs[i] <0.0 || Iobs[i] > 1.0 )
73             vis[i] = 0.0;
74         mean_intensities += float(vis[i]* Iobs[i]);
75         sum_weights      += float(vis[i]);
76     }
77     std::vector<float> temp_histogram(8,0.125f);
78 
79     float sum = 1.0;
80     for (unsigned i=0;i<Iobs.size();i++)
81     {
82         unsigned index = i + 1;
83         if (i == Iobs.size()-1)
84             index =0;
85         auto gradI = (float)std::fabs(Iobs[i]-Iobs[index]);
86 
87         int bin_index = (int) std::floor(gradI*8);
88         bin_index = bin_index>7 ? 7:bin_index;
89         temp_histogram[bin_index] += (float)std::min(vis[i],vis[index]);
90         sum += (float)std::min(vis[i],vis[index]);
91     }
92     for (unsigned i =0; i < 8;i++) temp_histogram[i] /= sum;
93     float entropy_histo  =0.0;
94     for (unsigned int i = 0; i<8; ++i)
95         entropy_histo += temp_histogram[i]*std::log(temp_histogram[i]);
96 
97     entropy_histo /= float(vnl_math::log2e);
98     entropy_histo = std::exp(-entropy_histo);
99 
100     i = 0;
101     pro.set_output_val<float>(i++, entropy_histo);
102 
103     return true;
104 }
105