Program Listing for File PrepKeypointData.cxx¶
↰ Return to documentation for file (larflow/KeyPoints/PrepKeypointData.cxx)
#include "PrepKeypointData.h"
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/ndarrayobject.h>
#include <sstream>
#include <algorithm>
#include <queue>
#include "TH2D.h"
#include "ublarcvapp/MCTools/MCPixelPGraph.h"
#include "ublarcvapp/MCTools/crossingPointsAnaMethods.h"
#include "larcv/core/DataFormat/IOManager.h"
#include "larcv/core/DataFormat/Image2D.h"
#include "larcv/core/DataFormat/EventImage2D.h"
// larlite
#include "LArUtil/SpaceChargeMicroBooNE.h"
#include "LArUtil/LArProperties.h"
#include "LArUtil/Geometry.h"
#include "DataFormat/storage_manager.h"
#include "DataFormat/mctrack.h"
#include "DataFormat/mcshower.h"
#include "DataFormat/mctruth.h"
namespace larflow {
namespace keypoints {
bool PrepKeypointData::_setup_numpy = false;
PrepKeypointData::PrepKeypointData()
: _adc_image_treename("wire"),
_label_tree(nullptr)
{
_nclose = 0;
_nfar = 0;
hdist[0] = new TH1F("hdist_x","",2002,-500,500.0);
hdist[1] = new TH1F("hdist_y","",2002,-500,500.0);
hdist[2] = new TH1F("hdist_z","",2002,-500,500.0);
hdpix[0] = new TH1F("hdpix_dt","",1001,-500,500);
hdpix[1] = new TH1F("hdpix_du","",1001,-500,500);
hdpix[2] = new TH1F("hdpix_dv","",1001,-500,500);
hdpix[3] = new TH1F("hdpix_dy","",1001,-500,500);
}
PrepKeypointData::~PrepKeypointData()
{
for (int v=0; v<3; v++ )
if ( hdist[v] ) delete hdist[v];
for (int v=0; v<4; v++ )
if ( hdpix[v] ) delete hdpix[v];
if ( _label_tree ) delete _label_tree;
}
void PrepKeypointData::process( larcv::IOManager& iolcv, larlite::storage_manager& ioll )
{
auto ev_adc = (larcv::EventImage2D*)iolcv.get_data(larcv::kProductImage2D,_adc_image_treename);
auto ev_segment = (larcv::EventImage2D*)iolcv.get_data(larcv::kProductImage2D,"segment");
auto ev_instance = (larcv::EventImage2D*)iolcv.get_data(larcv::kProductImage2D,"instance");
auto ev_ancestor = (larcv::EventImage2D*)iolcv.get_data(larcv::kProductImage2D,"ancestor");
auto ev_mctrack = (larlite::event_mctrack*)ioll.get_data( larlite::data::kMCTrack, "mcreco" );
auto ev_mcshower = (larlite::event_mcshower*)ioll.get_data( larlite::data::kMCShower, "mcreco" );
auto ev_mctruth = (larlite::event_mctruth*)ioll.get_data( larlite::data::kMCTruth, "generator" );
std::vector<larcv::Image2D> badch_v;
for ( auto const& img : ev_adc->Image2DArray() ) {
larcv::Image2D blank(img.meta());
blank.paint(0.0);
badch_v.emplace_back( std::move(blank) );
}
std::cout << "[PrepKeypointData Inputs]" << std::endl;
std::cout << " adc images: " << ev_adc->Image2DArray().size() << std::endl;
std::cout << " badch images: " << badch_v.size() << std::endl;
std::cout << " segment images: " << ev_segment->Image2DArray().size() << std::endl;
std::cout << " instance images: " << ev_instance->Image2DArray().size() << std::endl;
std::cout << " ancestor images: " << ev_ancestor->Image2DArray().size() << std::endl;
std::cout << " mctracks: " << ev_mctrack->size() << std::endl;
std::cout << " mcshowers: " << ev_mcshower->size() << std::endl;
std::cout << " mctruths: " << ev_mctruth->size() << std::endl;
_run = iolcv.event_id().run();
_subrun = iolcv.event_id().subrun();
_event = iolcv.event_id().event();
process( ev_adc->Image2DArray(),
badch_v,
ev_segment->Image2DArray(),
ev_instance->Image2DArray(),
ev_ancestor->Image2DArray(),
*ev_mctrack,
*ev_mcshower,
*ev_mctruth );
}
void PrepKeypointData::process( const std::vector<larcv::Image2D>& adc_v,
const std::vector<larcv::Image2D>& badch_v,
const std::vector<larcv::Image2D>& segment_v,
const std::vector<larcv::Image2D>& instance_v,
const std::vector<larcv::Image2D>& ancestor_v,
const larlite::event_mctrack& mctrack_v,
const larlite::event_mcshower& mcshower_v,
const larlite::event_mctruth& mctruth_v ) {
// allocate space charge class
larutil::SpaceChargeMicroBooNE sce;
// make particle graph
ublarcvapp::mctools::MCPixelPGraph mcpg;
mcpg.buildgraph( adc_v, segment_v, instance_v, ancestor_v,
mcshower_v, mctrack_v, mctruth_v );
//mcpg.printGraph();
// build key-points
_kpd_v.clear();
for (int i=0; i<(int)larflow::kNumKeyPoints; i++)
_kppos_v[i].clear();
// build crossing points for muon track primaries
std::vector<KPdata> track_kpd
= getMuonEndpoints( mcpg, adc_v, mctrack_v, &sce );
std::cout << "[Track Endpoint Results]" << std::endl;
for ( auto const& kpd : track_kpd ) {
std::cout << " " << kpd.str() << std::endl;
_kpd_v.emplace_back( std::move(kpd) );
}
// add points for shower starts
std::vector<KPdata> shower_kpd
= getShowerStarts( mcpg, adc_v, mcshower_v, &sce );
std::cout << "[Shower Endpoint Results]" << std::endl;
for ( auto const& kpd : shower_kpd ) {
std::cout << " " << kpd.str() << std::endl;
_kpd_v.emplace_back( std::move(kpd) );
}
// we change the kptype to neutrino vertex for those on it
_label_nu_keypoints( mctruth_v, adc_v, &sce, _kpd_v );
// filter duplicates
//filter_duplicates();
// copy positions of keypoints into flat vector for storage
for ( auto const& kpd : _kpd_v ) {
if ( kpd.kptype>=0 && kpd.kptype<larflow::kNumKeyPoints ) {
_kppos_v[ kpd.kptype ].push_back( kpd.keypt );
}
else {
throw std::runtime_error("unrecognized keypoint type");
}
}
}
std::vector<KPdata>
PrepKeypointData::getMuonEndpoints( ublarcvapp::mctools::MCPixelPGraph& mcpg,
const std::vector<larcv::Image2D>& adc_v,
const larlite::event_mctrack& mctrack_v,
larutil::SpaceChargeMicroBooNE* psce )
{
bool verbose = false;
// get list of primaries
std::vector<ublarcvapp::mctools::MCPixelPGraph::Node_t*> primaries
= mcpg.getPrimaryParticles();
// output vector of keypoint data
std::vector<KPdata> kpd_v;
for ( auto const& pnode : primaries ) {
if ( abs(pnode->pid)!=13
&& abs(pnode->pid)!=2212
&& abs(pnode->pid)!=211 )
continue;
auto const& mctrk = mctrack_v.at( pnode->vidx );
int crossingtype =
ublarcvapp::mctools::CrossingPointsAnaMethods::
doesTrackCrossImageBoundary( mctrk,
adc_v.front().meta(),
4050.0,
psce );
if ( crossingtype>=0 ) {
if ( crossingtype>=0 && crossingtype<=2) {
KPdata kpd;
kpd.crossingtype = crossingtype;
kpd.trackid = pnode->tid;
kpd.pid = pnode->pid;
kpd.vid = pnode->vidx;
kpd.is_shower = 0;
kpd.origin = pnode->origin;
kpd.kptype = larflow::kTrackEnds;
kpd.imgcoord =
ublarcvapp::mctools::CrossingPointsAnaMethods::getFirstStepPosInsideImage( mctrk, adc_v.front().meta(),
4050.0, true, 0.3, 0.1,
kpd.keypt, psce, verbose );
if ( kpd.imgcoord.size()>0 ) {
kpd_v.emplace_back( std::move(kpd) );
}
}
if ( crossingtype>=0 && crossingtype<=2 ) {
KPdata kpd;
kpd.crossingtype = crossingtype;
kpd.trackid = pnode->tid;
kpd.pid = pnode->pid;
kpd.vid = pnode->vidx;
kpd.is_shower = 0;
kpd.origin = pnode->origin;
kpd.kptype = larflow::kTrackEnds;
kpd.imgcoord =
ublarcvapp::mctools::CrossingPointsAnaMethods::getFirstStepPosInsideImage( mctrk, adc_v.front().meta(),
4050.0, false, 0.3, 0.1,
kpd.keypt, psce, verbose );
if ( kpd.imgcoord.size()>0 ) {
kpd_v.emplace_back( std::move(kpd) );
}
}
}//if track in image
}//end of primary loop
return kpd_v;
}
std::vector<KPdata>
PrepKeypointData::getShowerStarts( ublarcvapp::mctools::MCPixelPGraph& mcpg,
const std::vector<larcv::Image2D>& adc_v,
const larlite::event_mcshower& mcshower_v,
larutil::SpaceChargeMicroBooNE* psce )
{
// output vector of keypoint data
std::vector<KPdata> kpd_v;
// loop over nodes, look for electron/gamma pixels
for ( auto const& pnode : mcpg.node_v ) {
if ( abs(pnode.pid)!=11
&& abs(pnode.pid)!=22 )
continue;
int max_plane_pixels = 0;
for (auto const& pix_v : pnode.pix_vv ) {
if ( max_plane_pixels<pix_v.size() )
max_plane_pixels = pix_v.size();
}
if ( max_plane_pixels<10 )
continue;
// start: pnode.start; //should be in apparent position already
KPdata kpd;
kpd.crossingtype = 2;
kpd.trackid = pnode.tid;
kpd.pid = pnode.pid;
kpd.vid = pnode.vidx;
kpd.origin = pnode.origin;
kpd.is_shower = 1;
kpd.kptype = larflow::kShowerStart;
kpd.keypt.resize(3,0);
for (int i=0; i<3; i++)
kpd.keypt[i] = pnode.imgpos4[i];
std::vector<double> dpos(3,0);
for (int i=0; i<3; i++ ) dpos[i] = pnode.imgpos4[i];
kpd.imgcoord.resize(4,0);
try {
for (int p=0; p<3; p++)
kpd.imgcoord[1+p] = (int)larutil::Geometry::GetME()->NearestWire( dpos, p );
}
catch (...) {
continue;
}
float tick = pnode.imgpos4[3];
if ( tick>adc_v[0].meta().min_y() && tick<adc_v[0].meta().max_y() ) {
kpd.imgcoord[0] = adc_v[0].meta().row( tick );
}
else {
continue;
}
kpd_v.emplace_back( std::move(kpd) );
}//end of node loop
return kpd_v;
}
void PrepKeypointData::_label_nu_keypoints( const larlite::event_mctruth& mctruth_v,
const std::vector<larcv::Image2D>& img_v,
larutil::SpaceChargeMicroBooNE* psce,
std::vector<KPdata>& kpdata_v )
{
// loop over all interactions
int inu = -1;
for ( auto const& mct : mctruth_v ) {
inu++;
auto const& nu = mct.GetNeutrino();
auto const& nutraj = nu.Nu().Trajectory();
if (nutraj.size()>0) {
// get the space-charge corrected neutrino vertex
std::vector<double> nupos(3,0);
for (int i=0; i<3; i++ )
nupos[i] = nutraj.front().Position()[i];
std::vector<double> offsets = psce->GetPosOffsets( nupos[0], nupos[1], nupos[2] );
nupos[0] = nupos[0] - offsets[0] + 0.7;
nupos[1] += offsets[1];
nupos[2] += offsets[2];
// make a neutrino keypoint
KPdata kpd;
kpd.crossingtype = 0;
kpd.trackid = 0;
kpd.pid = 12;
kpd.vid = inu;
kpd.is_shower = 0;
kpd.origin = 1;
kpd.kptype = larflow::kNuVertex;
kpd.keypt.resize(3,0);
for (int i=0; i<3; i++) kpd.keypt[i] = nupos[i];
kpd.imgcoord.resize(4,0);
try {
for (int p=0; p<3; p++)
kpd.imgcoord[1+p] = (int)larutil::Geometry::GetME()->NearestWire( nupos, p );
}
catch (...) {
continue;
}
float tick = 3200 + nupos[0]/larutil::LArProperties::GetME()->DriftVelocity()/0.5;
if ( tick>img_v[0].meta().min_y() && tick<img_v[0].meta().max_y() ) {
kpd.imgcoord[0] = img_v[0].meta().row( tick );
}
else {
continue;
}
kpdata_v.push_back( kpd );
}//end of if neutrino truth object has trajectory point for vertex
}//end of loop over mctruth elements
}
void PrepKeypointData::filter_duplicates()
{
// first count the number of unique points
std::set< std::vector<int> > unique_coords;
std::vector< std::vector<int> > kpd_index;
int npts = 0;
for ( size_t ikpd=0; ikpd<_kpd_v.size(); ikpd++ ) {
auto const& kpd = _kpd_v[ikpd];
if (kpd.imgcoord.size()>0) {
if ( unique_coords.find( kpd.imgcoord )==unique_coords.end() ) {
kpd_index.push_back( std::vector<int>{(int)ikpd,0} );
unique_coords.insert( kpd.imgcoord );
npts++;
}
}
}
std::vector<KPdata> kpd_v;
for ( auto const& kpdidx : kpd_index ) {
kpd_v.emplace_back( std::move( _kpd_v[kpdidx[0]] ) );
}
std::swap(kpd_v,_kpd_v);
}
PyObject* PrepKeypointData::get_keypoint_array( int iclass ) const
{
if ( !PrepKeypointData::_setup_numpy ) {
import_array1(0);
PrepKeypointData::_setup_numpy = true;
}
// first count the number of unique points
std::set< std::vector<int> > unique_coords;
std::vector< std::vector<int> > kpd_index;
int npts = 0;
for ( size_t ikpd=0; ikpd<_kpd_v.size(); ikpd++ ) {
auto const& kpd = _kpd_v[ikpd];
if ( kpd.kptype!=(larflow::KeyPoint_t)iclass ) continue;
if (kpd.imgcoord.size()>0) {
if ( unique_coords.find( kpd.imgcoord )==unique_coords.end() ) {
kpd_index.push_back( std::vector<int>{(int)ikpd,0} );
unique_coords.insert( kpd.imgcoord );
npts++;
}
}
}
int nd = 2;
npy_intp dims[] = { npts, 10 };
PyArrayObject* array = (PyArrayObject*)PyArray_SimpleNew( nd, dims, NPY_FLOAT );
size_t ipt = 0;
for ( auto& kpdidx : kpd_index ) {
auto const& kpd = _kpd_v[kpdidx[0]];
if ( kpdidx[1]==0 ) {
// start img coordinates
for ( size_t i=0; i<4; i++ )
*((float*)PyArray_GETPTR2(array,ipt,i)) = (float)kpd.imgcoord[i];
// 3D point
for ( size_t i=0; i<3; i++ )
*((float*)PyArray_GETPTR2(array,ipt,4+i)) = (float)kpd.keypt[i];
// is shower
*((float*)PyArray_GETPTR2(array,ipt,7)) = (float)kpd.is_shower;
// origin
*((float*)PyArray_GETPTR2(array,ipt,8)) = (float)kpd.origin;
// PID
*((float*)PyArray_GETPTR2(array,ipt,9)) = (float)kpd.pid;
ipt++;
}
}// end of loop over keypointdata structs
return (PyObject*)array;
}
PyObject* PrepKeypointData::get_triplet_score_array( float sig ) const
{
if ( !PrepKeypointData::_setup_numpy ) {
import_array1(0);
PrepKeypointData::_setup_numpy = true;
}
// get label info for each triplet proposal
int npts = (int)_match_proposal_labels_v[0].size();
for (size_t iclass=0; iclass<larflow::kNumKeyPoints; iclass++) {
if ( _match_proposal_labels_v[iclass].size()!=npts ) {
throw std::runtime_error("number of triplet labels/scores for each class does not match!");
}
}
int nd = 2;
npy_intp dims[] = { npts, larflow::kNumKeyPoints };
PyArrayObject* array = (PyArrayObject*)PyArray_SimpleNew( nd, dims, NPY_FLOAT );
size_t ipt = 0;
for ( size_t ipt=0; ipt<npts; ipt++ ) {
for (size_t iclass=0; iclass<larflow::kNumKeyPoints; iclass++) {
auto const& label_v = _match_proposal_labels_v[iclass][ipt];
if ( label_v[0]==0.0 ) {
*((float*)PyArray_GETPTR2(array,ipt,iclass)) = 0.0;
}
else {
float dist = 0.;
for (int i=0; i<3; i++) dist += label_v[1+i]*label_v[1+i];
*((float*)PyArray_GETPTR2(array,ipt,iclass)) = exp( -0.5*dist/(sig*sig) );
}
}
}// end of loop over keypointdata structs
return (PyObject*)array;
}
void PrepKeypointData::make_proposal_labels( const larflow::prep::PrepMatchTriplets& match_proposals )
{
for (int i=0; i<larflow::kNumKeyPoints; i++) {
_match_proposal_labels_v[i].clear();
_match_proposal_labels_v[i].reserve(match_proposals._triplet_v.size());
}
for (int imatch=0; imatch<match_proposals._triplet_v.size(); imatch++ ) {
// triplet index (index of the sparse matrix coordinates)
const std::vector<int>& triplet = match_proposals._triplet_v[imatch];
// 3D position formed by intersection of the wires
const std::vector<float>& pos = match_proposals._pos_v[imatch];
// std::cout << "[match " << imatch << "] "
// << "testpt=(" << pos[0] << "," << pos[1] << "," << pos[2] << ") "
// << std::endl;
// make a score for each class
for (int ikpclass=0; ikpclass<(int)larflow::kNumKeyPoints; ikpclass++) {
// store label values
// [0]: has a match to a true keypoint
// [1,2,3]: (dx,dy,dz) to closest keypoint
// [4,5,6,7]: (dr,du,dv,dy) shift in row and columns
std::vector<float> label_v(10,0);
float dist = 1.0e9;
// dumb assignment, loops over all truth keypoints
// seems fast enough
const KPdata* kpd = nullptr;
std::vector<float> leafpos(3,0);
for (auto const& testkpd : _kpd_v ) {
// ignore those not in class
if ( testkpd.kptype!=(larflow::KeyPoint_t)ikpclass )
continue;
float testdist = 0.;
for ( int v=0; v<3; v++ )
testdist += (testkpd.keypt[v]-pos[v])*(testkpd.keypt[v]-pos[v]);
testdist = sqrt(testdist);
if ( dist>testdist ) {
// update the leadpos and kpd pointer
dist = testdist;
for (int v=0; v<3; v++ )
leafpos[v] = testkpd.keypt[v];
kpd = &testkpd;
}
}//end of loop over true points
// make label vector
// within 50 pixels/15 cm
if ( dist<0.3*50 ) {
label_v[0] = 1.0;
_nclose++;
}
else {
label_v[0] = 0.0;
_nfar++;
}
// make shift in 3D label
if ( dist<0.3*50 ) {
for (int i=0; i<3; i++ ) {
label_v[1+i] = leafpos[i]-pos[i];
hdist[i]->Fill(label_v[1+i]);
}
// shift in imgcoords
std::vector<int> imgcoords(4,0);
imgcoords[0] = match_proposals._sparseimg_vv[0][triplet[0]].row;
for (int p=0; p<3; p++ ) {
imgcoords[1+p] = match_proposals._sparseimg_vv[p][triplet[p]].col;
}
for (int i=0; i<4; i++) {
label_v[4+i] = imgcoords[i]-kpd->imgcoord[i];
hdpix[i]->Fill( label_v[4+i] );
}
}
_match_proposal_labels_v[ikpclass].push_back(label_v);
}//end of keypoint class loop
}//end of match proposal loop
}
void PrepKeypointData::writeHists()
{
for (int i=0; i<3; i++ ) {
hdist[i]->Write();
}
for (int i=0; i<4; i++ ) {
hdpix[i]->Write();
}
}
void PrepKeypointData::defineAnaTree()
{
_label_tree = new TTree("keypointlabels","Key point Training Labels");
_label_tree->Branch("run",&_run,"run/I");
_label_tree->Branch("subrun",&_subrun,"subrun/I");
_label_tree->Branch("event",&_event,"event/I");
_label_tree->Branch("kplabel_nuvertex", &_match_proposal_labels_v[0]);
_label_tree->Branch("kplabel_trackends", &_match_proposal_labels_v[1]);
_label_tree->Branch("kplabel_showerstart", &_match_proposal_labels_v[2]);
_label_tree->Branch("kppos_nuvertex", &_kppos_v[0]);
_label_tree->Branch("kppos_trackends", &_kppos_v[1]);
_label_tree->Branch("kppos_showerstart", &_kppos_v[2]);
}
void PrepKeypointData::writeAnaTree()
{
if (_label_tree)
_label_tree->Write();
}
void PrepKeypointData::printKeypoints() const
{
std::cout << "[PrepKeypointData::printKeypoints] -----------------" << std::endl;
for ( size_t i=0; i<_kpd_v.size(); i++ ) {
auto const& kpd = _kpd_v[i];
std::cout << " [" << i << "] "
<< "(" << kpd.keypt[0] << "," << kpd.keypt[1] << "," << kpd.keypt[2] << ")"
<< std::endl;
}
std::cout << "----------------------------------------------------" << std::endl;
}
std::vector<TH2D> PrepKeypointData::makeScoreImage( const int ikpclass, const float sigma,
const std::string histname,
const larflow::prep::PrepMatchTriplets& tripmaker,
const std::vector<larcv::Image2D>& adc_v ) const
{
std::vector<TH2D> hist_v;
for ( size_t p=0; p<adc_v.size(); p++ ) {
std::stringstream ss;
ss << histname << "_p" << (int)p;
TH2D hist( ss.str().c_str(), ss.str().c_str(),
adc_v[p].meta().cols(), adc_v[p].meta().min_x(), adc_v[p].meta().max_x(),
adc_v[p].meta().rows(), adc_v[p].meta().min_y(), adc_v[p].meta().max_y() );
for (size_t ipt=0; ipt<tripmaker._triplet_v.size(); ipt++) {
int r = tripmaker._sparseimg_vv[p][ tripmaker._triplet_v[ipt][p] ].row;
int c = tripmaker._sparseimg_vv[p][ tripmaker._triplet_v[ipt][p] ].col;
auto const& label_v = _match_proposal_labels_v[ikpclass][ipt];
if ( label_v[0]==0.0 && hist.GetBinContent(c+1,r+1)<0.01 ) {
hist.SetBinContent( c+1, r+1, 0.01 );
}
else if (label_v[0]>0.0) {
float dist = 0.;
for (int i=0; i<3; i++) dist += label_v[1+i]*label_v[1+i];
float score = exp( -0.5*dist/(sigma*sigma) );
if ( hist.GetBinContent(c+1,r+1)<score )
hist.SetBinContent( c+1, r+1, score );
}
}
hist_v.emplace_back( std::move(hist) );
}
return hist_v;
}
}
}