.. _program_listing_file_larflow_FlowContourMatch_makecontourflowmatches.cxx:

Program Listing for File makecontourflowmatches.cxx
===================================================

|exhale_lsh| :ref:`Return to documentation for file <file_larflow_FlowContourMatch_makecontourflowmatches.cxx>` (``larflow/FlowContourMatch/makecontourflowmatches.cxx``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "makecontourflowmatches.h"
   
   #include "larcv/core/Base/larcv_logger.h"
   #include "larcv/core/ROOTUtil/ROOTUtils.h"
   
   #include "TGraph.h"
   #include "TH2D.h"
   #include "TBox.h"
   
   #include "opencv2/imgproc.hpp"
   
   namespace larflow {
   
     void createMatchData( const ublarcvapp::ContourClusterAlgo& contour_data,
                           const larcv::Image2D& src_adc_full,
                           const larcv::Image2D& tar_adc_full,
                           const larcv::Image2D& flow_img_crop,
                           const larcv::Image2D& src_adc_crop,
                           const larcv::Image2D& tar_adc_crop,
                           ContourFlowMatchDict_t& matchdict,
                           const float threshold,
                           const float max_dist_to_target_contour,
                           const larcv::msg::Level_t verbosity,
                           const bool visualize )
     {
                           
       larcv::logger log("createMatchData");
       log.set(verbosity);
       
       int src_planeid = src_adc_crop.meta().plane();
       int tar_planeid = tar_adc_crop.meta().plane();
       // if ( log.debug() )
       //   log.send(larcv::msg::kDEBUG,__FUNCTION__,__LINE__)
       //     << " source plane=" << src_planeid << " "
       //     << " target plane=" << tar_planeid
       //     << std::endl;
   
       const larcv::ImageMeta& srcfullmeta = src_adc_full.meta();
       const larcv::ImageMeta& tarfullmeta = tar_adc_full.meta();    
       const larcv::ImageMeta& srcmeta = src_adc_crop.meta();
       const larcv::ImageMeta& tarmeta = tar_adc_crop.meta();
       
       // allocate arrays for image pixel to contour index lookup
       if ( !matchdict.index_map_initialized ) {
         matchdict.src_ctr_index_map.resize( src_adc_full.meta().rows()*src_adc_full.meta().cols(), -1 );
         matchdict.tar_ctr_index_map.resize( tar_adc_full.meta().rows()*tar_adc_full.meta().cols(), -1 );
         matchdict.index_map_initialized = true;
       }
       if ( !matchdict.src_ctr_pixel_v_initialized )
         matchdict.src_ctr_pixel_v.resize( contour_data.m_plane_atomics_v.at(src_planeid).size() );
   
       // if we visualize, we fill points here
       std::vector<float> vis_row[2];
       std::vector<float> vis_col[3];
         
       int nsrcpix_in_ctr = 0; // number of source image pixels inside a contour
       int ntarpix_in_ctr = 0; // number of target image pixels inside a contour
       int nflow_into_ctr = 0; // number of pixels that flow into a contour. unlikely zero, so here for a check
       for ( int r=0; r<(int)srcmeta.rows(); r++) {
         
         // for each row, we find the available contours on the images.
         // saves us search each time
   
         std::set< int >  tar_ctr_ids;     // (crop) index of target contours in this row
         std::vector<int> src_cols_in_ctr; // columns on the source (crop) image that have ctrs
         std::map<int,int> src_cols2ctrid; // map from (crop) src col to src ctr
         //src_cols_in_ctr.reserve(20);
   
         // std::cout << "------------------------------------------" << std::endl;
         // std::cout << "Find row=" << r << " contours" << std::endl;
   
         // Find contours on source image in this row
         // std::cout << "source: ";      
         for ( int c=0; c<(int)srcmeta.cols(); c++) {
       if ( src_adc_crop.pixel(r,c)<threshold )
         continue;
   
           float full_row = srcfullmeta.row( srcmeta.pos_y(r) );
           float full_col = srcfullmeta.col( srcmeta.pos_x(c) );
       cv::Point pt( full_col, full_row ); // point in the full image
   
           // search the contours
           bool found_contour = false;
           float closestdist  = -20000;
           for ( size_t ictr=0; ictr<contour_data.m_plane_atomics_v[src_planeid].size(); ictr++ ) {
             auto const& ctr     = contour_data.m_plane_atomics_v[src_planeid][ictr];
             auto const& ctrmeta = contour_data.m_plane_atomicmeta_v[src_planeid][ictr];
   
             // bbox check
             if ( pt.x < ctrmeta.getMinX() || pt.x > ctrmeta.getMaxX() ) continue;
             if ( pt.y < ctrmeta.getMinY() || pt.y > ctrmeta.getMaxY() ) continue;
             
         double result =  cv::pointPolygonTest( ctr, pt, true );
             if ( closestdist < result )
               closestdist = result;
             
         if ( result>-1 ) {
           src_cols_in_ctr.push_back( c );
           src_cols2ctrid[c] = ictr;
           //std::cout << " source pix (" << r << "," << c << ") "
               //          << "found in ctr=" << ictr << ". nelems=" << src_cols_in_ctr.size() << std::endl;
               
               // store in lookup map for crop
               int pixindex = full_row*srcfullmeta.cols() + full_col;
               matchdict.src_ctr_pixel_v.at( ictr ).push_back( pixindex );
           matchdict.src_ctr_index_map[ pixindex ] = ictr;
               nsrcpix_in_ctr++;
               found_contour = true;
   
               if ( visualize ) {
                 // source points with charge
                 vis_row[0].push_back( srcmeta.pos_y(r) );
                 vis_col[0].push_back( full_col );
               }
               
           break;
   
             }
       }//end of contour loop
           if ( !found_contour && log.debug() ) {
             log.send( larcv::msg::kDEBUG,__FUNCTION__,__LINE__ ) << "did not find source contour for point (col,row)=(" << pt.x << "," << pt.y << ")! "
                                                                  << " closest dist=" << closestdist << std::endl;
           }
         }//end of cropped column loop
         //std::cout << std::endl;
         
         // Find Contours on the target image in this row
         //std::cout << "target: ";      
         for ( int c=0; c<(int)tarmeta.cols(); c++) {
       if ( tar_adc_crop.pixel(r,c)<threshold )
         continue;
   
           float full_row = tarfullmeta.row( tarmeta.pos_y(r) );
           float full_col = tarfullmeta.col( tarmeta.pos_x(c) );
           
       cv::Point pt( full_col, full_row );
           
           // search the contours
           float closest_dist = -1000000;
           bool found_contour = false;
           for ( size_t ictr=0; ictr<contour_data.m_plane_atomics_v[tar_planeid].size(); ictr++ ) {
             
             auto const& ctr     = contour_data.m_plane_atomics_v[tar_planeid][ictr];
             auto const& ctrmeta = contour_data.m_plane_atomicmeta_v[tar_planeid][ictr];
   
             // bbox check
             if ( pt.x < ctrmeta.getMinX() || pt.x > ctrmeta.getMaxX() ) continue;
             if ( pt.y < ctrmeta.getMinY() || pt.y > ctrmeta.getMaxY() ) continue;          
       
         double result =  cv::pointPolygonTest( ctr, pt, true );
             if ( result > closest_dist )
               closest_dist = result;
             
         if ( result>=-1 ) {
           tar_ctr_ids.insert( ictr );
               matchdict.tar_ctr_index_map[ full_row*tarfullmeta.cols() + full_col ] = ictr;
               ntarpix_in_ctr++;
               found_contour = true;
           break;
         }
       }//end of target contour loop
           if ( !found_contour && log.debug() ) {
             log.send( larcv::msg::kDEBUG,__FUNCTION__,__LINE__ )
               << "did not find target contour for point (col,row)=(" << pt.x << "," << pt.y << ")! "
               << " localpix=(" << c << "," << r << ") "
               << " closest dist=" << closest_dist << std::endl;
           }
         }//end of target col loop
         //std::cout << std::endl;
         
         // Nothing in this row, move on to the next row
         if ( src_cols_in_ctr.size()==0 || tar_ctr_ids.size()==0 ) {
           // if ( log.debug() )
           //   log.send(larcv::msg::kDEBUG,__FUNCTION__,__LINE__)
           //     << "nothing to match. "
           //     << " srcwires=" << src_cols_in_ctr.size()
           //     << " tarwires=" << tar_ctr_ids.size()
           //     << std::endl;
           continue;
         }
         else if ( src_cols_in_ctr.size()>0 && tar_ctr_ids.size()>0 ) {
           // if ( log.debug() )
           //   log.send(larcv::msg::kDEBUG,__FUNCTION__,__LINE__)
           //     << " srcwires=" << src_cols_in_ctr.size()
           //     << " tarwires=" << tar_ctr_ids.size()
           //     << std::endl;        
         }
   
         // now loop over source columns in contours and make matches to target contours
         for ( auto const& source_crop_col : src_cols_in_ctr ) {
   
       float flow     = flow_img_crop.pixel(r,source_crop_col);
       int target_crop_col = source_crop_col+flow;
           //std::cout << " source-crop-col flow (" << r << "," << source_crop_col << ") flow=" << flow << " target-crop-col=" << target_crop_col << std::endl;
           if ( target_crop_col>= tarmeta.cols() )
             continue; // flow out of bounds
           if ( target_crop_col<0 )
             continue; // flow out of bounds
           float tar_row_full = tarfullmeta.row( tarmeta.pos_y(r) );
           float tar_col_full = tarfullmeta.col( tarmeta.pos_x(target_crop_col) );
           cv::Point tar_pt( tar_col_full, tar_row_full );
       
       // retrieve the contour we're in
       int src_ctr_id = src_cols2ctrid[source_crop_col];
       
       // loop through the target contours, find the contour closest to the flowed-to target wire
           int closest_contour_idx = -1;
           float closest_dist = -50;
       for ( auto const& ctrid : tar_ctr_ids ) {
             // have to move point back to full image coordinates
             
         float dist = cv::pointPolygonTest( contour_data.m_plane_atomics_v[tar_planeid][ctrid], tar_pt, true );
             // if dist is negatice, its outside the contour
             // if dist is positive, its inside the contour
   
             if ( dist>=0 ) {
               // inside, no need to search further
               closest_contour_idx = ctrid;
               closest_dist        = 0.;
               nflow_into_ctr++;
               break;
             }
             else if ( dist>-1*max_dist_to_target_contour ) {
               // enforce a max distance to a contour
               if ( closest_contour_idx<0 || closest_dist < dist ) {
                 closest_contour_idx = ctrid;
                 closest_dist        = dist;
               }
             }
           }//end of contour loop
   
           // if no contour found, move on
           //std::cout << "source -> tar contour index=" << closest_contour_idx << ". closest dist=" << closest_dist << std::endl;        
           if ( closest_contour_idx<0 ) {
             continue;
           }
   
   
           if ( visualize ) {
             vis_row[1].push_back( tarmeta.pos_y(r) );
             vis_col[1].push_back( tar_col_full );
           }
           
           // else we store/append a src-tar contour match for this source pixel
           
           // store the match data
           // --------------------
   
           // create a pair key
           SrcTarPair_t idpair = { src_ctr_id, closest_contour_idx };
   
           // check if we have one for this pair already
           auto it_indexmap = matchdict.find( idpair );
           if ( it_indexmap==matchdict.end() ) {
         // if the map doesn't have the pair we're looking for, we create the data
             ContourFlowMatch_t x( src_ctr_id,  closest_contour_idx);
             matchdict.insert( std::pair<SrcTarPair_t,ContourFlowMatch_t>(idpair,x) );
             it_indexmap = matchdict.find(idpair);
           }
           
           
           // store the specific pixel flow information
           float src_full_wire = srcmeta.pos_x(source_crop_col);
           float tar_full_wire = tarmeta.pos_x(target_crop_col);
           float in_contour_wire = tar_full_wire;
           
           if ( closest_dist<0 ) {
             // if we missed a contour, we find the offset that moves us into it
             const ublarcvapp::Contour_t& tar_ctr = contour_data.m_plane_atomics_v[tar_planeid][closest_contour_idx];
             bool shift_pos = false;
             float shifted_pos  = 0;
             for (int ioffset=-1; ioffset<=1; ioffset+=2) {
               float tar_col_full_shifted = tar_full_wire + float(ioffset)*fabs(closest_dist);
               float testdist = cv::pointPolygonTest( tar_ctr, cv::Point( tar_col_full_shifted, tar_row_full ), true );
               if ( testdist>=0 || (testdist>-max_dist_to_target_contour && fabs(testdist)<fabs(closest_dist)) ) {
                 shifted_pos = tar_col_full_shifted;
                 //std::cout << "shifted position better: dist=" << closest_dist << " to " << testdist << std::endl;
                 closest_dist = testdist;
                 shift_pos = true;
               }
             }
             // replace, as we did better
             if ( shift_pos )
               in_contour_wire = shifted_pos;
           }
   
           if ( visualize )
             vis_col[2].push_back( in_contour_wire );
             
           ContourFlowMatch_t::FlowPixel_t flowpix;
           flowpix.src_wire = src_full_wire;
           flowpix.tar_wire = in_contour_wire;
           flowpix.tar_orig = tar_full_wire;
           flowpix.row      = tar_row_full;
           flowpix.tick     = tarfullmeta.pos_y( tar_row_full );
           flowpix.pred_miss = std::fabs(closest_dist);
           flowpix.dist2cropcenter = std::fabs(source_crop_col - (float)srcmeta.cols()/2);
           int srcindex = (int)flowpix.row * srcfullmeta.cols() + (int)flowpix.src_wire;
           std::vector<ContourFlowMatch_t::FlowPixel_t>& flowpix_v = it_indexmap->second.getFlowPixelList( srcindex );
           flowpix_v.emplace_back( std::move(flowpix) );
   
           // debug
           if ( log.debug() )
             log.send( larcv::msg::kDEBUG,__FUNCTION__,__LINE__)
               << "creating a src-tar datapoint: "
               << "[src=" << src_full_wire << " -> "
               << "tar(orig)=" << tar_full_wire << " tar(shift)=" << in_contour_wire << "] "
               << "dist2ctr=" << closest_dist << "(orig) "
               << std::endl;
           
         }//end of loop over source cols
       }//end of loop over rows
       if ( log.info() )
         log.send( larcv::msg::kINFO,__FUNCTION__,__LINE__ )
           << "Number of pixels inside crop that flowed into a contour: "
           << nflow_into_ctr
           << std::endl;
   
       if ( visualize ) {
         // we plot all three planes.
         log.send( larcv::msg::kDEBUG,__FUNCTION__,__LINE__ )
           << "number of src pts, " << vis_row[0].size() << ", "
           << " tar pts, " << vis_row[1].size() << "."
           << std::endl;
         log.send( larcv::msg::kDEBUG,__FUNCTION__,__LINE__ )
           << "source meta: " << src_adc_crop.meta().dump()
           << std::endl;
         log.send( larcv::msg::kDEBUG,__FUNCTION__,__LINE__ )
           << "target meta: " << tar_adc_crop.meta().dump()
           << std::endl;
         
         TCanvas c("cflow","FlowContourMatch",1400,600);
         c.Divide(2,1);
         TH2D hadc_source = larcv::rootutils::as_th2d( src_adc_full, "hadc_source" );
         TH2D hadc_target = larcv::rootutils::as_th2d( tar_adc_full, "hadc_target" );
         TGraph gsrcpts( vis_row[0].size(), vis_col[0].data(), vis_row[0].data() );
         TGraph gtarpts( vis_row[1].size(), vis_col[1].data(), vis_row[1].data() );
         TGraph gtarcontourpts( vis_col[2].size(), vis_col[2].data(), vis_row[1].data() );
   
         TBox src_box( src_adc_crop.meta().min_x(), src_adc_crop.meta().min_y(),
                       src_adc_crop.meta().max_x(), src_adc_crop.meta().max_y() );
         TBox tar_box( tar_adc_crop.meta().min_x(), tar_adc_crop.meta().min_y(),
                       tar_adc_crop.meta().max_x(), tar_adc_crop.meta().max_y() );
         src_box.SetLineColor(kRed);
         src_box.SetLineWidth(2);
         src_box.SetFillStyle(0);
         tar_box.SetLineColor(kRed);
         tar_box.SetLineWidth(2);      
         tar_box.SetFillStyle(0);
   
         c.Draw();      
   
         gsrcpts.SetMarkerStyle(24);
   
         gtarpts.SetMarkerStyle(24);
   
         gtarcontourpts.SetLineColor(kRed);      
         gtarcontourpts.SetMarkerColor(kRed);
         gtarcontourpts.SetMarkerStyle(24);
   
         c.cd(1);
         hadc_source.Draw("colz");
         gsrcpts.Draw("p");
         src_box.Draw("same");
         
         c.cd(2);
         hadc_target.Draw("colz");
         gtarpts.Draw("psame");
         gtarcontourpts.Draw("psame");
         tar_box.Draw("same");
   
         c.Update();
   
         TCanvas ccheck("cropcheck","Crop Check",1400,600);
         ccheck.Divide(3,1);
         ccheck.Draw();
         ccheck.cd(1);
         TH2D hsrc_crop = larcv::rootutils::as_th2d( src_adc_crop, "hsrc_crop" );
         hsrc_crop.Draw("colz");
         ccheck.cd(2);
         TH2D htar_crop = larcv::rootutils::as_th2d( tar_adc_crop, "htar_crop" );
         htar_crop.Draw("colz");
         ccheck.cd(3);
         TH2D hflow_crop = larcv::rootutils::as_th2d( flow_img_crop, "hflow_crop" );
         //hflow_crop.GetZaxis()->SetRangeUser(-832,832);      
         hflow_crop.Draw("colz");
         ccheck.Update();
   
         char cname[100];
         sprintf( cname, "flowcheck_src%d_to_tar%d_%d_%d.png",
                  (int)src_adc_crop.meta().plane(),(int)tar_adc_crop.meta().plane(),
                  (int)src_adc_crop.meta().min_x(),(int)src_adc_crop.meta().min_y() );
         c.SaveAs(cname);
   
         //std::cout << "[enter] to continue." << std::endl;
         //std::cin.get();
       }
   
     }
   
   }