.. _program_listing_file_larflow_Reco_cluster_functions.cxx:

Program Listing for File cluster_functions.cxx
==============================================

|exhale_lsh| :ref:`Return to documentation for file <file_larflow_Reco_cluster_functions.cxx>` (``larflow/Reco/cluster_functions.cxx``)

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

.. code-block:: cpp

   #include "cluster_functions.h"
   
   #include <fstream>
   
   #include "ublarcvapp/dbscan/DBScan.h"  
   #include "ublarcvapp/dbscan/sDBScan.h" 
   #include "ublarcvapp/ContourTools/ContourClusterAlgo.h"
   #include <cilantro/principal_component_analysis.hpp>
   
   namespace larflow {
   namespace reco {
   
     ClusterFunctions::ClusterFunctions() {}
     
     void cluster_larflow3dhits( const std::vector<larlite::larflow3dhit>& hit_v,
                                 std::vector< cluster_t >& cluster_v,
                                 const float maxdist, const int minsize, const int maxkd )
     {
       
       // convert points into list of floats
       std::vector< std::vector<float> > points_v;
       points_v.reserve( hit_v.size() );
       
       for ( auto const& lfhit : hit_v ) {
         std::vector<float > hit = { (float)lfhit[0], (float)lfhit[1], (float)lfhit[2] };
         points_v.push_back( hit );
       }
   
   
       clock_t start = clock();  
       std::vector< ublarcvapp::dbscan::dbCluster > dbcluster_v = ublarcvapp::dbscan::DBScan::makeCluster3f( maxdist, minsize, maxkd, points_v );
   
       
       for (int ic=0; ic<(int)dbcluster_v.size()-1;ic++) {
         // skip the last cluster, which are noise points
         auto const& cluster = dbcluster_v[ic];
         cluster_t c;
         c.points_v.reserve(cluster.size());
         c.imgcoord_v.reserve(cluster.size());
         c.hitidx_v.reserve(cluster.size());
         for ( auto const& hitidx : cluster ) {
           // store 3d position and 2D image coordinates
           c.points_v.push_back( points_v.at(hitidx) );
           std::vector<int> coord(4,0);
           coord[3] = hit_v[hitidx].tick;
           coord[0] = (int)hit_v[hitidx].targetwire[0]; // U-plane
           coord[1] = (int)hit_v[hitidx].targetwire[1]; // V-plane
           coord[2] = (int)hit_v[hitidx].srcwire;
           c.imgcoord_v.push_back( coord );
           c.hitidx_v.push_back(hitidx);
         }
         cluster_v.emplace_back(std::move(c));
       }
       clock_t end = clock();
       double elapsed = double(end-start)/CLOCKS_PER_SEC;
       
       std::cout << "[cluster_larflow3dhit] made clusters: " << dbcluster_v.size() << " elpased=" << elapsed << " secs" << std::endl;    
     }
   
     void cluster_spacepoint_v( const std::vector< std::vector<float> >& points_v,
                                std::vector< cluster_t >& cluster_v,
                                const float maxdist, const int minsize, const int maxkd )
     {
   
       clock_t start = clock();  
       std::vector< ublarcvapp::dbscan::dbCluster > dbcluster_v
         = ublarcvapp::dbscan::DBScan::makeCluster3f( maxdist, minsize, maxkd, points_v );
   
       
       for (int ic=0; ic<(int)dbcluster_v.size()-1;ic++) {
         // skip the last cluster, which are noise points
         auto const& cluster = dbcluster_v[ic];
         cluster_t c;
         c.points_v.reserve(cluster.size());
         //c.imgcoord_v.reserve(cluster.size());
         c.hitidx_v.reserve(cluster.size());
         for ( auto const& hitidx : cluster ) {
           // store 3d position and 2D image coordinates
           c.points_v.push_back( points_v.at(hitidx) );
           c.hitidx_v.push_back(hitidx);
         }
         cluster_v.emplace_back(std::move(c));
       }
       clock_t end = clock();
       double elapsed = double(end-start)/CLOCKS_PER_SEC;
       
       std::cout << "[cluster_spacepoint] made clusters: " << dbcluster_v.size() << " elpased=" << elapsed << " secs" << std::endl;    
         
     }
   
     
     void cluster_sdbscan_larflow3dhits( const std::vector<larlite::larflow3dhit>& hit_v,
                                         std::vector< cluster_t >& cluster_v,
                                         const float maxdist, const int minsize, const int maxkd ) {
       
       clock_t start = clock();
       
       // convert points into list of floats
       std::vector< std::vector<float> > points_v;
       points_v.reserve( hit_v.size() );
       
       for ( auto const& lfhit : hit_v ) {
         std::vector<float> hit = { lfhit[0], lfhit[1], lfhit[2] };
         points_v.push_back( hit );
       }
       
       auto sdbscan = ublarcvapp::dbscan::SDBSCAN< std::vector<float>, float >();
       sdbscan.Run( &points_v, 3, maxdist, minsize );
       
       auto noise = sdbscan.Noise;
       auto dbcluster_v = sdbscan.Clusters;
       
       for (int ic=0; ic<(int)dbcluster_v.size();ic++) {
         // skip the last cluster, which are noise points
         auto const& cluster = dbcluster_v[ic];
         cluster_t c;
         c.points_v.reserve(cluster.size());
         c.imgcoord_v.reserve(cluster.size());
         c.hitidx_v.reserve(cluster.size());
         for ( auto const& hitidx : cluster ) {
           // store 3d position and 2D image coordinates
           c.points_v.push_back( points_v.at(hitidx) );
           std::vector<int> coord(4,0);
           coord[3] = hit_v[hitidx].tick;
           coord[0] = (int)hit_v[hitidx].targetwire[0]; // U-plane
           coord[1] = (int)hit_v[hitidx].targetwire[1]; // V-plane
           coord[2] = (int)hit_v[hitidx].targetwire[2]; // Y-plane
           c.imgcoord_v.push_back( coord );
           c.hitidx_v.push_back(hitidx);
         }
         cluster_v.emplace_back(std::move(c));
       }
       clock_t end = clock();
       double elapsed = double(end-start)/CLOCKS_PER_SEC;
       
       std::cout << "[cluster_simple_larflow3dhits] made clusters: " << dbcluster_v.size() << " elpased=" << elapsed << " secs" << std::endl;    
     }
   
     void cluster_sdbscan_spacepoints( const std::vector< std::vector<float> >& hit_v,
                                       std::vector< cluster_t >& cluster_v,
                                       const float maxdist, const int minsize, const int maxkd )
     {
       
       clock_t start = clock();
       
       // convert points into list of floats
       std::vector< std::vector<float> > points_v;
       points_v.reserve( hit_v.size() );
       
       for ( auto const& lfhit : hit_v ) {
         std::vector<float> hit = { lfhit[0], lfhit[1], lfhit[2] };
         points_v.push_back( hit );
       }
       
       auto sdbscan = ublarcvapp::dbscan::SDBSCAN< std::vector<float>, float >();
       sdbscan.Run( &points_v, 3, maxdist, minsize );
       
       auto noise = sdbscan.Noise;
       auto dbcluster_v = sdbscan.Clusters;
       
       for (int ic=0; ic<(int)dbcluster_v.size();ic++) {
         // skip the last cluster, which are noise points
         auto const& cluster = dbcluster_v[ic];
         cluster_t c;
         c.points_v.reserve(cluster.size());
         c.imgcoord_v.reserve(cluster.size());
         c.hitidx_v.reserve(cluster.size());
         for ( auto const& hitidx : cluster ) {
           // store 3d position and 2D image coordinates
           c.points_v.push_back( points_v.at(hitidx) );
           c.hitidx_v.push_back(hitidx);        
         }
         cluster_v.emplace_back(std::move(c));
       }
       clock_t end = clock();
       double elapsed = double(end-start)/CLOCKS_PER_SEC;
       
       std::cout << "[cluster_sdbscan_spacepoints] made clusters: " << dbcluster_v.size() << " elpased=" << elapsed << " secs" << std::endl;    
     }
     
     void cluster_pca( cluster_t& cluster ) {
   
       //std::cout << "[cluster_pca]" << std::endl;
       cluster.bbox_v.resize(3);
       for (size_t i=0; i<3; i++ ) {
         cluster.bbox_v[i].resize(2,0);
         cluster.bbox_v[i][0] = 1.0e9;  // initial min value
         cluster.bbox_v[i][1] = -1.0e9; // initial max value
       }
       
       std::vector< Eigen::Vector3f > eigen_v;
       eigen_v.reserve( cluster.points_v.size() );
       for ( auto const& hit : cluster.points_v ) {
         // store values for PCA 
         eigen_v.push_back( Eigen::Vector3f( hit[0], hit[1], hit[2] ) );
         // also determine axis-aligned bounding box values
         for (int i=0; i<3; i++ ) {
           if ( hit[i]<cluster.bbox_v[i][0] ) cluster.bbox_v[i][0] = hit[i];
           if ( hit[i]>cluster.bbox_v[i][1] ) cluster.bbox_v[i][1] = hit[i];
         }
       }
   
       
       cilantro::PrincipalComponentAnalysis3f pca( eigen_v );
       cluster.pca_center.resize(3,0);
       cluster.pca_eigenvalues.resize(3,0);
       cluster.pca_axis_v.clear();
       for (int i=0; i<3; i++) {
         cluster.pca_center[i] = pca.getDataMean()(i);
         cluster.pca_eigenvalues[i] = pca.getEigenValues()(i);
         std::vector<float> e_v = { pca.getEigenVectors()(0,i),
                                    pca.getEigenVectors()(1,i),
                                    pca.getEigenVectors()(2,i) };
         cluster.pca_axis_v.push_back( e_v );
       }
       // we provide a way to order the points
       struct pca_coord {
         int idx;
         float s;
         bool operator<( const pca_coord& rhs ) const {
           if ( s<rhs.s) return true;
           return false;
         };
       };
   
       // normalize first axis (should be though)
       for ( int ipca=0; ipca<3; ipca++) {
         float lenpc = 0.;
         for (int i=0; i<3; i++ ) {
           lenpc += cluster.pca_axis_v[ipca][i]*cluster.pca_axis_v[ipca][i];
         }
         if ( lenpc>0 ) {
           lenpc = sqrt(lenpc);
           for ( int i=0; i<3; i++ )
             cluster.pca_axis_v[ipca][i] /= lenpc;
         }
       }
       
       std::vector< pca_coord > ordered_v;
       ordered_v.reserve( cluster.points_v.size() );
       for ( size_t idx=0; idx<cluster.points_v.size(); idx++ ) {
         float s = 0.;
         for (int i=0; i<3; i++ ) {
           s += (cluster.points_v[idx][i]-cluster.pca_center[i])*cluster.pca_axis_v[0][i];
         }
         pca_coord coord;
         coord.idx = idx;
         coord.s   = s;
         ordered_v.push_back( coord );
       }
       std::sort( ordered_v.begin(), ordered_v.end() );
       cluster.ordered_idx_v.resize( ordered_v.size() );
       cluster.pca_proj_v.resize( ordered_v.size() );
       for ( size_t i=0; i<ordered_v.size(); i++ ) {
         cluster.ordered_idx_v[i] = ordered_v[i].idx;
         cluster.pca_proj_v[i]    = ordered_v[i].s;
       }
   
       // define ends
       cluster.pca_ends_v.resize(2);
       cluster.pca_ends_v[0].resize(3,0);
       cluster.pca_ends_v[1].resize(3,0);
       for (int i=0; i<3; i++) {
         cluster.pca_ends_v[0][i] = cluster.pca_center[i] + cluster.pca_proj_v.front()*cluster.pca_axis_v[0][i];
         cluster.pca_ends_v[1][i] = cluster.pca_center[i] + cluster.pca_proj_v.back()*cluster.pca_axis_v[0][i];
       }
   
       float len3 = 0;
       std::vector<float> d3(3);    
       for (int i=0; i<3; i++ ) {
         d3[i] = cluster.pca_ends_v[1][i] - cluster.pca_ends_v[0][i];
         len3 += d3[i]*d3[i];      
       }
       len3 = sqrt(len3);    
       cluster.pca_len = len3;
         
       // define radius of each point to the pca-axis
       float max_r = 0.;
       float ave_r2 = 0.;
       cluster.pca_radius_v.resize( ordered_v.size(), 0 );
       for ( size_t i=0; i<ordered_v.size(); i++ ) {
         // get distance of point from pca-axis
         // http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
         std::vector<float>& pt = cluster.points_v[ ordered_v[i].idx ];
   
         std::vector<float> d1(3);
         std::vector<float> d2(3);
   
         float len1 = 0.;
         for (int i=0; i<3; i++ ) {
           d1[i] = pt[i] - cluster.pca_ends_v[0][i];
           d2[i] = pt[i] - cluster.pca_ends_v[1][i];
           len1 += d1[i]*d1[i];
         }
         len1 = sqrt(len1);
   
         if ( len3<1.0e-4 ) {
           cluster.pca_radius_v[i] = len1;
         }
         else {
   
           // cross-product
           std::vector<float> d1xd2(3);
           d1xd2[0] =  d1[1]*d2[2] - d1[2]*d2[1];
           d1xd2[1] = -d1[0]*d2[2] + d1[2]*d2[0];
           d1xd2[2] =  d1[0]*d2[1] - d1[1]*d2[0];
           float len1x2 = 0.;
           for ( int i=0; i<3; i++ ) {
             len1x2 += d1xd2[i]*d1xd2[i];
           }
           len1x2 = sqrt(len1x2);
           cluster.pca_radius_v[i]  = len1x2/len3; // distance of point from PCA-axis
         }
         ave_r2 += cluster.pca_radius_v[i]*cluster.pca_radius_v[i];
         if ( max_r < cluster.pca_radius_v[i] )
           max_r = cluster.pca_radius_v[i];
         
       }//end of loop over ordered points
   
       ave_r2 = sqrt( ave_r2/cluster.pca_radius_v.size() );
   
       cluster.pca_max_r  = max_r;
       cluster.pca_ave_r2 = ave_r2;
       
       
     }
   
     void cluster_runpca( std::vector<cluster_t>& cluster_v ) {
       clock_t start = clock();
       for ( auto& cluster : cluster_v ) {
         cluster_pca( cluster );
       }
       clock_t end = clock();
       double elapsed = double(end-start)/CLOCKS_PER_SEC;
       std::cout << "[ cluster_runpca ] elapsed=" << elapsed << " secs" << std::endl;
     }
   
     void cluster_dump2jsonfile( const std::vector<cluster_t>& cluster_v, std::string outfilename ) {
       nlohmann::json j;
       std::vector< nlohmann::json > jcluster_v;
       j["clusters"] = jcluster_v;
   
       for ( auto const& cluster : cluster_v ) {
         //std::cout << "cluster: nhits=" << cluster.points_v.size() << std::endl;
         nlohmann::json jcluster = cluster_json(cluster);
         j["clusters"].emplace_back( std::move(jcluster) );      
       }
       
       std::ofstream o(outfilename.c_str());
       j >> o;
       o.close();
     }
   
     nlohmann::json cluster_json( const cluster_t& cluster ) {
   
       nlohmann::json jcluster;
       jcluster["hits"] = cluster.points_v;
   
       // save start, center, end of pca to make line
       std::vector< std::vector<float> > pca_points(3);
       for (int i=0; i<3; i++ )
         pca_points[i].resize(3,0);
   
       for (int i=0; i<3; i++ ) {
         pca_points[0][i] = cluster.pca_ends_v[0][i];// - 5.0*cluster.pca_axis_v[0][i];
         pca_points[1][i] = cluster.pca_center[i];
         pca_points[2][i] = cluster.pca_ends_v[1][i];// + 5.0*cluster.pca_axis_v[0][i];
       }
       jcluster["pca"] = pca_points;
       return jcluster;
     }
       
   
     void cluster_splitbytrackshower( const std::vector<larlite::larflow3dhit>& hit_v,
                                      const std::vector<larcv::Image2D>& ssnettrack_image_v,
                                      std::vector<larlite::larflow3dhit>& track_hit_v,
                                      std::vector<larlite::larflow3dhit>& shower_hit_v,
                                      float min_larmatch_score ) {
   
       clock_t begin = clock();
       
       track_hit_v.clear();
       shower_hit_v.clear();
       track_hit_v.reserve(  hit_v.size() );
       shower_hit_v.reserve( hit_v.size() );
   
       std::vector< const larcv::ImageMeta* > meta_v( ssnettrack_image_v.size(),0);
       for ( size_t p=0; p<ssnettrack_image_v.size(); p++ )
         meta_v[p] = &(ssnettrack_image_v[p].meta());
   
       int below_threshold = 0.;
       
       for ( auto const & hit : hit_v ) {
   
         //std::cout << "hit[9]=" << hit[9] << std::endl;
         if ( min_larmatch_score>0 && hit.size()>=10 && hit[9]<min_larmatch_score ) {
           below_threshold++;
           continue;
         }
         
         std::vector<float> scores(3,0);
         scores[0] = ssnettrack_image_v[0].pixel( meta_v[0]->row( hit.tick, __FILE__, __LINE__ ), hit.targetwire[0], __FILE__, __LINE__ );
         scores[1] = ssnettrack_image_v[1].pixel( meta_v[1]->row( hit.tick, __FILE__, __LINE__ ), hit.targetwire[1], __FILE__, __LINE__ );
         scores[2] = ssnettrack_image_v[2].pixel( meta_v[2]->row( hit.tick, __FILE__, __LINE__ ), hit.srcwire,       __FILE__, __LINE__ );
   
         // condition ... gather metrics
         int n_w_score = 0;
         float tot_score = 0.;
         float max_score = 0.;
         float min_non_zero = 1.;
         for ( auto s : scores ) {
           if ( s>0 ) n_w_score++;
           tot_score += s;
           if ( max_score<s )
             max_score = s;
           if ( s>1 && s<min_non_zero )
             min_non_zero = 0;
         }
   
         if ( n_w_score>0 && tot_score/float(n_w_score)>0.5 ) {
           track_hit_v.push_back( hit );
         }
         else
           shower_hit_v.push_back( hit );
       }
   
       clock_t end = clock();
       double elapsed = double(end-begin)/CLOCKS_PER_SEC;
       
       std::cout << "[cluster_split_by_trackshower]: "
                 << "original=" << hit_v.size()
                 << " into track=" << track_hit_v.size()
                 << " and shower=" << shower_hit_v.size()
                 << " below-threshold=" << below_threshold
                 << " elasped=" << elapsed << " secs"
                 << std::endl;
     }
   
     void cluster_imageprojection( const cluster_t& cluster,
                                   std::vector<larcv::Image2D>& clust2d_images_v ) {
   
       // zero the images
       for ( auto& img : clust2d_images_v )
         img.paint(0.0);
       
       for ( auto const& imgcoord : cluster.imgcoord_v ) {
         for ( auto& img : clust2d_images_v ) {
           img.set_pixel( img.meta().row( imgcoord[3], __FILE__, __LINE__ ),
                          img.meta().col( imgcoord[img.meta().plane()], __FILE__, __LINE__ ),
                          50.0 /* dummy value */ );
         }
       }
       
     }
   
     void cluster_getcontours( std::vector<larcv::Image2D>& clust2d_images_v ) {
       ublarcvapp::ContourClusterAlgo contour_algo;
       contour_algo.analyzeImages( clust2d_images_v, 10.0, 2, 5, 10, 10, 2 );
   
       // contours made. now what.
       for ( size_t p=0; p<3; p++ ) {
         int ncontours = contour_algo.m_plane_atomics_v[p].size();
         std::cout << "[cluster_getcontours] ncontours plane[" << p << "] = " << ncontours << std::endl;
         for (int ictr=0; ictr<ncontours; ictr++ ) {
   
           if ( !contour_algo.m_plane_atomicmeta_v[p][ictr].hasValidPCA() ) {
             std::cout << "  icontour[" << ictr << "] pca not valid" << std::endl;
             continue;
           }
           
           // dist from start-end of pca-axis
           float pca_dist = 0;
           float dx = 0;
           for ( size_t i=0; i<2; i++ ) {
             dx = ( contour_algo.m_plane_atomicmeta_v[p][ictr].getPCAEndPos()[i] -
                    contour_algo.m_plane_atomicmeta_v[p][ictr].getPCAStartPos()[i] );
             pca_dist += dx*dx;
           }
           pca_dist = sqrt(pca_dist);
           std::cout << "  icontour[" << ictr << "] "
                     << " pca-dist=" << pca_dist
                     << " eigenvals [0]=" << contour_algo.m_plane_atomicmeta_v[p][ictr].getPCAeigenvalue(0)
                     << " [1]=" << contour_algo.m_plane_atomicmeta_v[p][ictr].getPCAeigenvalue(1)
                     << std::endl;
         }
       }
       
     }
   
     float cluster_closest_endpt_dist( const cluster_t& clust_a, const cluster_t& clust_b,
                                       std::vector< std::vector<float> >& endpts ) {
       // 2x2 calculation
   
       float dist[4] = {0};
       for (int i=0; i<2; i++) {
         for (int j=0; j<2; j++) {
           for (int k=0; k<3; k++) {
             dist[ 2*i+j ] += (clust_a.pca_ends_v[i][k] - clust_b.pca_ends_v[j][k])*(clust_a.pca_ends_v[i][k] - clust_b.pca_ends_v[j][k]);
           }
         }
       }
   
       float mindist = -1.0;
       int mini = 0;
       int minj = 0;
       for (int k=0; k<4; k++ ) {
         if ( dist[k]<mindist || mindist<0 ) {
           mindist = dist[k];
           mini = k/2;
           minj = k%2;
         }
       }
   
       endpts.resize(2);
       endpts[0].resize(3,0);
       endpts[1].resize(3,0);
       for (int k=0; k<3; k++ ) {
         endpts[0][k] = clust_a.pca_ends_v[mini][k];
         endpts[1][k] = clust_b.pca_ends_v[minj][k];
       }
       
       if ( mindist>0 )
         mindist = sqrt(mindist);
       
       return mindist;
     }
   
     bool cluster_endpt_in_bbox( const cluster_t& clust_a, const cluster_t& clust_b ) {
   
       bool a_in_b = false;
       bool b_in_a = false;
       for (int iend=0; iend<2; iend++) {
         int dims_inside = 0;
         for (int i=0; i<3; i++ ) {
           if ( clust_a.pca_ends_v[iend][i]>=clust_b.bbox_v[i][0] && clust_a.pca_ends_v[iend][i]<=clust_b.bbox_v[i][1] )
             dims_inside++;
         }
         if ( dims_inside==3 )
           a_in_b = true;
         
         dims_inside = 0;
         for (int i=0; i<3; i++ ) {
           if ( clust_b.pca_ends_v[iend][i]>=clust_a.bbox_v[i][0] && clust_b.pca_ends_v[iend][i]<=clust_a.bbox_v[i][1] )
             dims_inside++;
         }
         if ( dims_inside==3 )
           b_in_a = true;
       }
   
       return a_in_b || b_in_a;
     }
   
     float cluster_cospca( const cluster_t& clust_a, const cluster_t& clust_b ) {
       float cospca = 0.;
       float lena = 0.;
       float lenb = 0.;
       for ( int k=0; k<3; k++ ) {
         cospca += clust_a.pca_axis_v[0][k]*clust_b.pca_axis_v[0][k];
       }
       return cospca;
     }
   
     cluster_t cluster_merge( const cluster_t& clust_a, const cluster_t& clust_b ) {
       // copy first
       cluster_t merge = clust_a;
   
       merge.points_v.reserve(   merge.points_v.size()+clust_b.points_v.size() );
       merge.imgcoord_v.reserve( merge.imgcoord_v.size()+clust_b.imgcoord_v.size() );
       merge.hitidx_v.reserve( merge.imgcoord_v.size()+clust_b.imgcoord_v.size() );
       for ( size_t i=0; i<clust_b.points_v.size(); i++ ) {
         merge.points_v.push_back(   clust_b.points_v[i] );
         merge.imgcoord_v.push_back( clust_b.imgcoord_v[i] );
         merge.hitidx_v.push_back( clust_b.hitidx_v[i] );
       }
   
       cluster_pca( merge );
       return merge;
     }
   
     larlite::pcaxis cluster_make_pcaxis( const cluster_t& c, int cidx ) {
   
       // pca-axis
       larlite::pcaxis::EigenVectors e_v;
       // just std::vector< std::vector<double> >
       // we store axes (3) and then the 1st axis end points. So five vectors.
       for ( auto const& a_v : c.pca_axis_v ) {
         std::vector<double> da_v = { (double)a_v[0], (double)a_v[1], (double) a_v[2] };
           e_v.push_back( da_v );
       }
       // start and end points
       for ( auto const& p_v : c.pca_ends_v ) {
         std::vector<double> dp_v = { (double)p_v[0], (double)p_v[1], (double)p_v[2] };
         e_v.push_back( dp_v );
       }
       double eigenval[3] = { c.pca_eigenvalues[0], c.pca_eigenvalues[1], c.pca_eigenvalues[2] };
       double centroid[3] = { c.pca_center[0], c.pca_center[1], c.pca_center[2] };
       larlite::pcaxis llpca( true, c.points_v.size(), eigenval, e_v, centroid, 0, cidx );
   
       return llpca;
     }
   
   
     cluster_t cluster_from_larflowcluster( const larlite::larflowcluster& lfcluster ) {
   
       std::cout << "[cluster_from_larflowcluster] input cluster size=" << lfcluster.size() << std::endl;
       cluster_t c;
       c.points_v.reserve( 2*lfcluster.size() );
       c.imgcoord_v.reserve( 2*lfcluster.size() );
       c.hitidx_v.reserve( 2*lfcluster.size() );
       for ( auto const& lfhit : lfcluster ) {
         //if ( lfhit.size()<3 ) continue;
         //if ( lfhit.targetwire.size()<3 ) continue;
         // std::cout << "lfhit: " << lfhit[0] << "," << lfhit[1] << "," << lfhit[2] << std::endl;
         // std::cout << "   coord=(" << lfhit.targetwire[0] << "," << lfhit.targetwire[1] << "," << lfhit.targetwire[2] << ")" << std::endl;
         // std::cout << "   tick=" << lfhit.tick << std::endl;
         std::vector<float> pos   = { lfhit[0], lfhit[1], lfhit[2] };
         std::vector<int>   coord = { lfhit.targetwire[0],
                                      lfhit.targetwire[1],
                                      lfhit.targetwire[2],
                                      lfhit.tick };
         c.points_v.push_back( pos );
         c.imgcoord_v.push_back( coord );
         c.hitidx_v.push_back( (int)c.points_v.size()-1 );
       }
   
       cluster_pca( c );
   
       return c;
     }
   
     /*
      * @brief project hits into the ADC image and get the total pixel sum for each plane
      *
      * @param[in] cluster Cluster whose hits we project. (Note we use cluster_t::imgcoord).
      * @param[in] img_v   Wire plane image to get the charge sum from.
      * @return charge sum on each plane.
      *
      */
     std::vector<float> cluster_pixelsum( const cluster_t& cluster,
                                          const std::vector<larcv::Image2D>& img_v  ) {
       
       std::vector<float> cluster_adc(3,0.0);
   
       std::vector<larcv::Image2D> blank_v;
       for (auto const& img : img_v ) {
         larcv::Image2D blank(img.meta());
         blank.paint(0);
         blank_v.emplace_back( std::move(blank) );
       }
   
       for ( size_t ihit=0; ihit<cluster.imgcoord_v.size(); ihit++ ) {
         int row = img_v[0].meta().row( cluster.imgcoord_v[ihit][3], __FILE__, __LINE__ );
         for ( size_t p=0; p<img_v.size(); p++ ) {
           int col = img_v[p].meta().col( cluster.imgcoord_v[ihit][p], __FILE__, __LINE__ );
           float used   = blank_v[p].pixel(row,col);
           if ( used<1.0 ) {
             float pixval = img_v[p].pixel(row,col);
             cluster_adc[p] += pixval;
             // mark it to not be reused
             blank_v[p].set_pixel(row,col,10.0);
           }
         }
       }
       return cluster_adc;
     }
   
     float cluster_dist_from_pcaline( const cluster_t& cluster, const std::vector<float>& pt ) {
   
       // 
       // http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
       
       std::vector<float> d1(3);
       std::vector<float> d2(3);
   
       float len1 = 0.;
       for (int i=0; i<3; i++ ) {
         d1[i] = pt[i] - cluster.pca_ends_v[0][i];
         d2[i] = pt[i] - cluster.pca_ends_v[1][i];
         len1 += d1[i]*d1[i];
       }
       len1 = sqrt(len1);
   
       if ( cluster.pca_len<1.0e-4 ) {
         // short cluster, use distance to end point
         return len1;
       }
   
       // cross-product
       std::vector<float> d1xd2(3);
       d1xd2[0] =  d1[1]*d2[2] - d1[2]*d2[1];
       d1xd2[1] = -d1[0]*d2[2] + d1[2]*d2[0];
       d1xd2[2] =  d1[0]*d2[1] - d1[1]*d2[0];
       float len1x2 = 0.;
       for ( int i=0; i<3; i++ ) {
         len1x2 += d1xd2[i]*d1xd2[i];
       }
       len1x2 = sqrt(len1x2);
       float r = len1x2/cluster.pca_len;
       return r;
     }
     
   }
   }