1 // This is brl/bseg/bstm/cpp/pro/processes/bstm_cpp_box_roc_process.cxx
2 #include <iostream>
3 #include <fstream>
4 #include <bprb/bprb_func_process.h>
5 //:
6 // \file
7 // \brief
8 //
9 // \author Ali Osman Ulusoy
10 // \date Nov 1, 2013
11 
12 #ifdef _MSC_VER
13 #  include "vcl_msvc_warnings.h"
14 #endif
15 #include <bstm/io/bstm_cache.h>
16 #include <bstm/io/bstm_lru_cache.h>
17 #include <bstm/bstm_scene.h>
18 #include <bstm/bstm_block.h>
19 #include <bstm/bstm_data_base.h>
20 //brdb stuff
21 #include <brdb/brdb_value.h>
22 
23 #include "vgl/vgl_box_3d.h"
24 #include "vgl/vgl_intersection.h"
25 
26 #include <bstm/bstm_util.h>
27 #include <bbas_pro/bbas_1d_array_float.h>
28 #include <boct/boct_bit_tree.h>
29 
30 
31 namespace bstm_cpp_box_roc_process_globals
32 {
33   constexpr unsigned n_inputs_ = 14;
34   constexpr unsigned n_outputs_ = 2;
35 
36   typedef unsigned char uchar;
37   typedef unsigned short ushort;
38   typedef vnl_vector_fixed<uchar, 16> uchar16;
39   typedef vnl_vector_fixed<uchar, 8> uchar8;
40   typedef vnl_vector_fixed<ushort, 4> ushort4;
41 
42 
43 }
44 
bstm_cpp_box_roc_process_cons(bprb_func_process & pro)45 bool bstm_cpp_box_roc_process_cons(bprb_func_process& pro)
46 {
47   using namespace bstm_cpp_box_roc_process_globals;
48 
49   //process takes 1 input
50   std::vector<std::string> input_types_(n_inputs_);
51 
52   input_types_[0] = "bstm_scene_sptr";
53   input_types_[1] = "bstm_cache_sptr";
54   //change region
55   input_types_[2] = "float"; //center x
56   input_types_[3] = "float"; //center y
57   input_types_[4] = "float"; //center z
58   input_types_[5] = "float"; //len x
59   input_types_[6] = "float"; //len y
60   input_types_[7] = "float"; //len z
61 
62   //box to ignore outside
63   input_types_[8] = "float"; //center x
64   input_types_[9] = "float"; //center y
65   input_types_[10] = "float"; //center z
66   input_types_[11] = "float"; //len x
67   input_types_[12] = "float"; //len y
68   input_types_[13] = "float"; //len z
69 
70   // process has 1 output:
71   // output[0]: scene sptr
72   std::vector<std::string> output_types;
73 
74   output_types.emplace_back("bbas_1d_array_float_sptr");  // tpr
75   output_types.emplace_back("bbas_1d_array_float_sptr");  // fpr
76   bool good = pro.set_input_types(input_types_) && pro.set_output_types(output_types);
77   return good;
78 }
79 
bstm_cpp_box_roc_process(bprb_func_process & pro)80 bool bstm_cpp_box_roc_process(bprb_func_process& pro)
81 {
82   using namespace bstm_cpp_box_roc_process_globals;
83 
84   if ( pro.n_inputs() < n_inputs_ ) {
85     std::cout << pro.name() << ": The input number should be " << n_inputs_<< std::endl;
86     return false;
87   }
88 
89   //get the inputs
90   unsigned i = 0;
91   bstm_scene_sptr scene =pro.get_input<bstm_scene_sptr>(i++);
92   bstm_cache_sptr cache= pro.get_input<bstm_cache_sptr>(i++);
93   auto change_region_center_x = pro.get_input<float>(i++);
94   auto change_region_center_y = pro.get_input<float>(i++);
95   auto change_region_center_z = pro.get_input<float>(i++);
96   auto change_region_len_x = pro.get_input<float>(i++);
97   auto change_region_len_y = pro.get_input<float>(i++);
98   auto change_region_len_z = pro.get_input<float>(i++);
99 
100   auto center_x = pro.get_input<float>(i++);
101   auto center_y = pro.get_input<float>(i++);
102   auto center_z = pro.get_input<float>(i++);
103   auto len_x = pro.get_input<float>(i++);
104   auto len_y = pro.get_input<float>(i++);
105   auto len_z = pro.get_input<float>(i++);
106   //create vgl box,
107   const vgl_point_3d<double> change_region_center(change_region_center_x,change_region_center_y,change_region_center_z);
108   vgl_box_3d<double> change_box(change_region_center,change_region_len_x,change_region_len_y,change_region_len_z, vgl_box_3d<double>::centre);
109   vgl_box_3d<double> change_box_insides(change_region_center,change_region_len_x,change_region_len_y,change_region_len_z, vgl_box_3d<double>::centre);
110   change_box_insides.scale_about_centroid(0.85);
111 
112   //create vgl box,
113   const vgl_point_3d<double> center(center_x,center_y,center_z);
114   vgl_box_3d<double> box(center,len_x,len_y,len_z, vgl_box_3d<double>::centre);
115 
116 
117   //init roc
118   unsigned const numPoints = 1000; //num points on ROC
119   auto * tp=new bbas_1d_array_float(numPoints);
120   auto * tn=new bbas_1d_array_float(numPoints);
121   auto * fp=new bbas_1d_array_float(numPoints);
122   auto * fn=new bbas_1d_array_float(numPoints);
123   for (unsigned int pnt=0; pnt<numPoints; ++pnt) {
124     tp->data_array[pnt]=0.0f;
125     fp->data_array[pnt]=0.0f;
126     tn->data_array[pnt]=0.0f;
127     fn->data_array[pnt]=0.0f;
128   }
129 
130 
131   //iterate over each block/metadata to check if bbox intersects the input bbox
132   std::map<bstm_block_id, bstm_block_metadata> blocks = scene->blocks();
133   std::map<bstm_block_id, bstm_block_metadata> ::const_iterator bstm_iter = blocks.begin();
134   for(; bstm_iter != blocks.end() ; bstm_iter++)
135   {
136     bstm_block_id bstm_id = bstm_iter->first;
137     bstm_block_metadata bstm_metadata = bstm_iter->second;
138 
139 
140     bstm_block* blk = cache->get_block(bstm_metadata.id_);
141     bstm_time_block* blk_t = cache->get_time_block(bstm_metadata.id_);
142 
143     bstm_data_base *change_buffer = cache->get_data_base(bstm_metadata.id_,bstm_data_traits<BSTM_CHANGE>::prefix()  );
144     auto* change_array = (bstm_data_traits<BSTM_CHANGE>::datatype*) change_buffer->data_buffer();
145 
146     //iterate through each tree
147     boxm2_array_3d<uchar16>&  trees = blk->trees();
148     for (unsigned int x = 0; x < trees.get_row1_count(); ++x) {
149       for (unsigned int y = 0; y < trees.get_row2_count(); ++y) {
150        for (unsigned int z = 0; z < trees.get_row3_count(); ++z) {
151          //load current block/tree
152          uchar16 tree = trees(x, y, z);
153          boct_bit_tree bit_tree((unsigned char*) tree.data_block());
154 
155          //first check if the tree box is contained in the box,
156          vgl_point_3d<double> tree_min_pt(bstm_metadata.local_origin_.x() + x*bstm_metadata.sub_block_dim_.x(),
157                                            bstm_metadata.local_origin_.y() + y*bstm_metadata.sub_block_dim_.y(),
158                                            bstm_metadata.local_origin_.z() + z*bstm_metadata.sub_block_dim_.z());
159          vgl_box_3d<double> tree_box(tree_min_pt,bstm_metadata.sub_block_dim_.x(),
160                                                   bstm_metadata.sub_block_dim_.y(),
161                                                   bstm_metadata.sub_block_dim_.z(),
162                                                   vgl_box_3d<double>::min_pos);
163          if(!vgl_intersection<double>(tree_box,box).is_empty() )
164          {
165 
166 
167            //iterate through leaves of the tree
168            std::vector<int> leafBits = bit_tree.get_leaf_bits();
169            std::vector<int>::iterator iter;
170            for (iter = leafBits.begin(); iter != leafBits.end(); ++iter) {
171              int curr_depth = bit_tree.depth_at((*iter));
172              double side_len = 1.0 / (double) (1<<curr_depth);
173              int data_offset = bit_tree.get_data_index(*iter);
174 
175              vgl_point_3d<double> localCenter = bit_tree.cell_center(*iter);
176              vgl_point_3d<double> cellCenter(localCenter.x() + x, localCenter.y() + y, localCenter.z() + z);
177              vgl_point_3d<double> cellCenter_global(   float(cellCenter.x()*bstm_metadata.sub_block_dim_.x() + bstm_metadata.local_origin_.x()),
178                                                         float(cellCenter.y()*bstm_metadata.sub_block_dim_.y() + bstm_metadata.local_origin_.y()),
179                                                         float(cellCenter.z()*bstm_metadata.sub_block_dim_.z() + bstm_metadata.local_origin_.z()));
180 
181              if(change_box_insides.contains(cellCenter_global))
182                continue;
183 
184              for (unsigned int pnt=0; pnt<numPoints; ++pnt) {
185                float detection_threshold = 1.0f-(float)pnt/(float)numPoints;
186                if(change_array[data_offset] >= detection_threshold)
187                  if(change_box.contains(cellCenter_global))
188                    tp->data_array[pnt]+= side_len;
189                  else
190                    fp->data_array[pnt]+= side_len;
191                else
192                  if(change_box.contains(cellCenter_global))
193                    fn->data_array[pnt]+= side_len;
194                  else
195                    tn->data_array[pnt]+= side_len;
196              }
197 
198            }
199 
200          }
201 
202        }
203       }
204     }
205   }
206 
207 
208   auto * tpr=new bbas_1d_array_float(numPoints);
209   auto * fpr=new bbas_1d_array_float(numPoints);
210 
211   for (unsigned int pnt=0; pnt<numPoints; ++pnt) {
212     tpr->data_array[pnt]= tp->data_array[pnt] / (tp->data_array[pnt] + fn->data_array[pnt]);
213     fpr->data_array[pnt]= fp->data_array[pnt] / (fp->data_array[pnt] + tn->data_array[pnt]);
214     std::cout << "tpr: " <<tpr->data_array[pnt] << " and fpr: " << fpr->data_array[pnt] << std::endl;
215   }
216 
217   pro.set_output_val<bbas_1d_array_float_sptr>(0, tpr);
218   pro.set_output_val<bbas_1d_array_float_sptr>(1, fpr);
219   return true;
220 
221 }
222