Skip to content
Snippets Groups Projects
Commit 6786c048 authored by Damien Thisse's avatar Damien Thisse
Browse files

Added a RANSAC method for ACTAR and some functions in LinearCluster3D

parent c5ff41e1
No related branches found
No related tags found
1 merge request!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...
Hide whitespace changes
Inline Side-by-side
NPLib/TrackReconstruction/NPLinearRansac3D.cxx
+ 198
− 2
View file @ 6786c048
...@@ -17,8 +17,14 @@ ...@@ -17,8 +17,14 @@
*****************************************************************************/ *****************************************************************************/
#include "NPLinearRansac3D.h" #include "NPLinearRansac3D.h"
#include<iostream>
using namespace std; using namespace std;
using namespace NPL; using namespace NPL;
ClassImp(LinearCluster3D);
ClassImp(LinearRansac3D);
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void LinearCluster3D::LinearFit() { void LinearCluster3D::LinearFit() {
...@@ -36,6 +42,7 @@ void LinearCluster3D::LinearFit() { ...@@ -36,6 +42,7 @@ void LinearCluster3D::LinearFit() {
double meanZ = 0; double meanZ = 0;
double w; double w;
for (unsigned int i = 0; i < sizeI; i++) { for (unsigned int i = 0; i < sizeI; i++) {
TVector3 Pi((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]); TVector3 Pi((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]);
for (unsigned int j = i + 1; j < sizeI; j++) { for (unsigned int j = i + 1; j < sizeI; j++) {
...@@ -69,11 +76,200 @@ void LinearCluster3D::LinearFit() { ...@@ -69,11 +76,200 @@ void LinearCluster3D::LinearFit() {
m_P0 = TVector3(meanX, meanY, meanZ); m_P0 = TVector3(meanX, meanY, meanZ);
m_Dir = (TVector3(meanDX, meanDY, meanDZ)).Unit(); m_Dir = (TVector3(meanDX, meanDY, meanDZ)).Unit();
m_isFitted = true;
// Cluster dir is arbitrary. we choose to always have positive Z dir // Cluster dir is arbitrary. we choose to always have positive Z dir
if (m_Dir.Z() < 0) if (m_Dir.Z() < 0)
m_Dir = -m_Dir; m_Dir = -m_Dir;
}; };
void LinearCluster3D::FastLinearFit() {
double xx=0, yy=0, zz=0, xx2=0, yy2=0, zz2=0, xy=0, xz=0, yz=0, qq=0;
unsigned int sizeI = m_index.size();
for(unsigned int i = 0; i < sizeI; i++){
xx+=(*m_X)[m_index[i]];
yy+=(*m_Y)[m_index[i]];
zz+=(*m_Z)[m_index[i]];
xx2+=(*m_X)[m_index[i]]*(*m_X)[m_index[i]];
yy2+=(*m_Y)[m_index[i]]*(*m_Y)[m_index[i]];
zz2+=(*m_Z)[m_index[i]]*(*m_Z)[m_index[i]];
xy+=(*m_X)[m_index[i]]*(*m_Y)[m_index[i]];
xz+=(*m_X)[m_index[i]]*(*m_Z)[m_index[i]];
yz+=(*m_Y)[m_index[i]]*(*m_Z)[m_index[i]];
qq+=1;
}
double Xm = xx/qq, Ym = yy/qq, Zm = zz/qq;
double Sxx = -Xm*Xm+xx2/qq, Syy = -Ym*Ym+yy2/qq, Szz = -Zm*Zm+zz2/qq;
double Sxy = -Xm*Ym+xy/qq, Sxz = -Xm*Zm+xz/qq, Syz = -Ym*Zm+yz/qq;
double theta = 0.5*atan(2*Sxy/(Sxx-Syy));
double ct = cos(theta), st = sin(theta), ct2 = ct*ct, st2 = st*st, cst = cos(theta)*sin(theta);
double K11 = (Syy+Szz)*ct2+(Sxx+Szz)*st2-2*Sxy*cst;
double K22 = (Syy+Szz)*st2+(Sxx+Szz)*ct2+2*Sxy*cst;
//double K12 = -Sxy*(ct2-st2)+(Sxx-Syy)*cst; //=0 by construction
double K10 = Sxz*ct+Syz*st;
double K01 = -Sxz*st+Syz*ct;
double K00 = Sxx+Syy;
double C2 = -K00-K11-K22;
double C1 = K00*K11+K00*K22+K11*K22-K01*K01-K10*K10;
double C0 = K01*K01*K11+K10*K10*K22-K00*K11*K22;
double p = C1-C2*C2/3, q = 2*C2*C2*C2/27-C1*C2/3+C0;
double R = q*q/4+p*p*p/27;
double dm2 = 0;
if(R>0)
{
dm2 = -C2/3+pow(-q/2+sqrt(R), 1./3.)+pow(-q/2-sqrt(R), 1./3.);
}
else
{
double rho = sqrt(-p*p*p/27), phi = acos(-q/(2*rho)); double R = q*q/4+p*p*p/27;
double sqrt3rho = pow(rho, 1./3.);
double d1 = -C2/3+2*sqrt3rho*cos(phi/3);
double d2 = -C2/3+2*sqrt3rho*cos((phi+2*M_PI)/3);
double d3 = -C2/3+2*sqrt3rho*cos((phi+4*M_PI)/3);
dm2 = std::min(std::min(d1, d2), d3);
}
double a = -K10*ct/(K11-dm2)+K01*st/(K22-dm2);
double b = -K10*st/(K11-dm2)-K01*ct/(K22-dm2);
double denom = 1./(1+a*a+b*b);
double u = ((1+b*b)*Xm-a*b*Ym+a*Zm)*denom;
double v = (-a*b*Xm+(1+a*a)*Ym+b*Zm)*denom;
double w = (a*Xm+b*Ym+(a*a+b*b)*Zm)*denom;
m_P0 = TVector3(Xm, Ym, Zm);
m_Dir = (TVector3(Xm-u, Ym-v, Zm-w)).Unit();
}
void LinearCluster3D::FastLinearFitShort() {
double xx=0, yy=0, zz=0, xx2=0, yy2=0, zz2=0, xy=0, xz=0, yz=0, qq=0;
unsigned int sizeI = m_index_short.size();
for(unsigned int i = 0; i < sizeI; i++){
xx+=(*m_X)[m_index_short[i]];
yy+=(*m_Y)[m_index_short[i]];
zz+=(*m_Z)[m_index_short[i]];
xx2+=(*m_X)[m_index_short[i]]*(*m_X)[m_index_short[i]];
yy2+=(*m_Y)[m_index_short[i]]*(*m_Y)[m_index_short[i]];
zz2+=(*m_Z)[m_index_short[i]]*(*m_Z)[m_index_short[i]];
xy+=(*m_X)[m_index_short[i]]*(*m_Y)[m_index_short[i]];
xz+=(*m_X)[m_index_short[i]]*(*m_Z)[m_index_short[i]];
yz+=(*m_Y)[m_index_short[i]]*(*m_Z)[m_index_short[i]];
qq+=1;
}
double Xm = xx/qq, Ym = yy/qq, Zm = zz/qq;
double Sxx = -Xm*Xm+xx2/qq, Syy = -Ym*Ym+yy2/qq, Szz = -Zm*Zm+zz2/qq;
double Sxy = -Xm*Ym+xy/qq, Sxz = -Xm*Zm+xz/qq, Syz = -Ym*Zm+yz/qq;
double theta = 0.5*atan(2*Sxy/(Sxx-Syy));
double ct = cos(theta), st = sin(theta), ct2 = ct*ct, st2 = st*st, cst = cos(theta)*sin(theta);
double K11 = (Syy+Szz)*ct2+(Sxx+Szz)*st2-2*Sxy*cst;
double K22 = (Syy+Szz)*st2+(Sxx+Szz)*ct2+2*Sxy*cst;
//double K12 = -Sxy*(ct2-st2)+(Sxx-Syy)*cst; //=0 by construction
double K10 = Sxz*ct+Syz*st;
double K01 = -Sxz*st+Syz*ct;
double K00 = Sxx+Syy;
double C2 = -K00-K11-K22;
double C1 = K00*K11+K00*K22+K11*K22-K01*K01-K10*K10;
double C0 = K01*K01*K11+K10*K10*K22-K00*K11*K22;
double p = C1-C2*C2/3, q = 2*C2*C2*C2/27-C1*C2/3+C0;
double R = q*q/4+p*p*p/27;
double dm2 = 0;
if(R>0)
{
dm2 = -C2/3+pow(-q/2+sqrt(R), 1./3.)+pow(-q/2-sqrt(R), 1./3.);
}
else
{
double rho = sqrt(-p*p*p/27), phi = acos(-q/(2*rho)); double R = q*q/4+p*p*p/27;
double sqrt3rho = pow(rho, 1./3.);
double d1 = -C2/3+2*sqrt3rho*cos(phi/3);
double d2 = -C2/3+2*sqrt3rho*cos((phi+2*M_PI)/3);
double d3 = -C2/3+2*sqrt3rho*cos((phi+4*M_PI)/3);
dm2 = std::min(std::min(d1, d2), d3);
}
double a = -K10*ct/(K11-dm2)+K01*st/(K22-dm2);
double b = -K10*st/(K11-dm2)-K01*ct/(K22-dm2);
double denom = 1./(1+a*a+b*b);
double u = ((1+b*b)*Xm-a*b*Ym+a*Zm)*denom;
double v = (-a*b*Xm+(1+a*a)*Ym+b*Zm)*denom;
double w = (a*Xm+b*Ym+(a*a+b*b)*Zm)*denom;
m_P0 = TVector3(Xm, Ym, Zm);
m_Dir = (TVector3(Xm-u, Ym-v, Zm-w)).Unit();
}
void LinearCluster3D::CheckDirection(const TVector3 & vertex){
if(m_FurthestPoint == TVector3(0, 0, 0)){
double max_distance = 0;
for(int i = 0; i < m_index.size(); i++)
{
TVector3 aPoint((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]);
double distance = (aPoint-vertex).Mag();
if(distance > max_distance){
max_distance = distance;
m_FurthestPoint = aPoint;
}
}
}
if((m_FurthestPoint-vertex).X()*m_Dir.X() < 0) this->Inverse();
}
void LinearCluster3D::CalculateTrackLength(const TVector3 & vertex, const double & percentRange){
TH1D* Hrange = new TH1D("Range", "range", 256, 0, 512);
TH1D* Hrangecumule = new TH1D("RangeCumulated", "rangecumule", 256, 0, 512);
double total_charge = 0;
double cumulated_charge = 0;
double max_distance = 0;
for(int i = 0; i < m_index.size(); i++)
{
TVector3 aPoint((*m_X)[m_index[i]], (*m_Y)[m_index[i]], (*m_Z)[m_index[i]]);
double distance = (aPoint-vertex).Mag();
total_charge+=(*m_Q)[m_index[i]];
Hrange->Fill(distance, (*m_Q)[m_index[i]]);
if(distance > max_distance){
max_distance = distance;
m_FurthestPoint = aPoint;
}
if(distance < m_ShortTrackMaxLength){
m_index_short.push_back(m_index[i]);
}
}
m_Charge = total_charge;
//Inverse the track if it is not going in the good direction
if((m_FurthestPoint-vertex).X()*m_Dir.X() < 0) this->Inverse();
for(int bin = 1; bin < Hrange->GetNbinsX(); bin++)
{
cumulated_charge+=Hrange->GetBinContent(bin);
Hrangecumule->SetBinContent(bin, cumulated_charge/total_charge);
}
Int_t nbin = 1;
while(Hrangecumule->GetBinContent(nbin) < percentRange && nbin < Hrangecumule->GetNbinsX()-1)
{
nbin++;
}
Double_t delta, deltabinap, deltabinav;
delta = Hrangecumule->GetBinContent(nbin)-Hrangecumule->GetBinContent(nbin-1);
deltabinap = Hrangecumule->GetBinContent(nbin)-percentRange;
deltabinav = percentRange-Hrangecumule->GetBinContent(nbin-1);
m_TrackLength = deltabinap*Hrangecumule->GetBinCenter(nbin-1)/delta+deltabinav*Hrangecumule->GetBinCenter(nbin)/delta;
delete Hrange;
delete Hrangecumule;
}
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<double>& Y, vector<double>& Z) { vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<double>& Y, vector<double>& Z) {
cluster_id.clear(); cluster_id.clear();
...@@ -92,7 +288,7 @@ vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<dou ...@@ -92,7 +288,7 @@ vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<dou
m_iteration_d = 0; m_iteration_d = 0;
m_iteration = 0; m_iteration = 0;
while (m_iteration++ < m_max_iteration && m_assigned.size() < sizeX) { while (m_iteration++ < m_max_iteration/* && m_assigned.size() < sizeX*/) {
LinearCluster3D cluster(&X, &Y, &Z); LinearCluster3D cluster(&X, &Y, &Z);
// take 2 distant point randomly that has not been match before // take 2 distant point randomly that has not been match before
d = 0; d = 0;
...@@ -123,7 +319,7 @@ vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<dou ...@@ -123,7 +319,7 @@ vector<LinearCluster3D> LinearRansac3D::TreatEvent(vector<double>& X, vector<dou
} }
} // while } // while
// loop over the cluster starting with the biggest // loop over the cluster starting with the biggest
unsigned int current_cluster = 0; unsigned int current_cluster = 0;
unsigned int index; unsigned int index;
......
NPLib/TrackReconstruction/NPLinearRansac3D.h
+ 45
− 5
View file @ 6786c048
...@@ -24,15 +24,25 @@ ...@@ -24,15 +24,25 @@
//root //root
#include "TVector3.h" #include "TVector3.h"
#include "TRandom2.h" #include "TRandom2.h"
#include "TH1.h"
// nptool // nptool
#include "NPTrackingUtility.h" #include "NPTrackingUtility.h"
namespace NPL{ namespace NPL{
///////////////////////////////// /////////////////////////////////
class LinearCluster3D{ class LinearCluster3D{
public: public:
LinearCluster3D(){}; LinearCluster3D(){
m_isFitted = false;
m_FurthestPoint = TVector3(0, 0, 0);
};
LinearCluster3D(std::vector<double>* X, std::vector<double>* Y, std::vector<double>* Z){ LinearCluster3D(std::vector<double>* X, std::vector<double>* Y, std::vector<double>* Z){
m_X=X; m_Y=Y; m_Z=Z; m_X=X; m_Y=Y; m_Z=Z; m_isFitted = false;
m_FurthestPoint = TVector3(0, 0, 0);
}
LinearCluster3D(std::vector<double>* X, std::vector<double>* Y, std::vector<double>* Z, std::vector<double>* Q){
m_X=X; m_Y=Y; m_Z=Z; m_Q=Q; m_isFitted = false;
m_FurthestPoint = TVector3(0, 0, 0);
} }
~LinearCluster3D(){}; ~LinearCluster3D(){};
...@@ -40,22 +50,48 @@ namespace NPL{ ...@@ -40,22 +50,48 @@ namespace NPL{
void AddIndex(const unsigned int& i) {m_index.push_back(i);}; void AddIndex(const unsigned int& i) {m_index.push_back(i);};
unsigned int size() const {return m_index.size();}; unsigned int size() const {return m_index.size();};
unsigned int GetIndex (const unsigned int& i) const {return m_index[i];}; unsigned int GetIndex (const unsigned int& i) const {return m_index[i];};
void SetShortTrackMaxLength(const double & l) {m_ShortTrackMaxLength = l;}
//Inverse the direction
void Inverse() {m_Dir = -m_Dir;}; //!
// Fit the cluster // Fit the cluster
void LinearFit(); void LinearFit(); //!
void FastLinearFit(); //!
void FastLinearFitShort(); //!
//Calculate the track length from the interaction vertex, and for a given percent of charge deposited along the way
//It also inverses the track if needed, and regroups the points below a certain distance from the vertex to refit it without the deviation
void CalculateTrackLength(const TVector3 & vertex, const double & percentRange = 0.95); //!
void CheckDirection(const TVector3 & vertex);
private: private:
std::vector<double>* m_X;//! std::vector<double>* m_X;//!
std::vector<double>* m_Y;//! std::vector<double>* m_Y;//!
std::vector<double>* m_Z;//! std::vector<double>* m_Z;//!
std::vector<double>* m_Q;//!
std::vector<unsigned int> m_index;//! std::vector<unsigned int> m_index;//!
TVector3 m_P0;// a point on the line std::vector<unsigned int> m_index_short;//!
TVector3 m_Dir;// direction of the line TVector3 m_P0; // a point on the line
TVector3 m_Dir; // direction of the line
double m_TrackLength; //!
double m_Charge; //!
// check if the track has been fitted with the "LinearFit" function
bool m_isFitted; //!
//Furthest point compared to vertex
TVector3 m_FurthestPoint; //!
double m_ShortTrackMaxLength; //!
public: public:
std::vector<unsigned int> GetIndex() const {return m_index;}; std::vector<unsigned int> GetIndex() const {return m_index;};
TVector3 GetP0()const {return m_P0;}; TVector3 GetP0()const {return m_P0;};
TVector3 GetDir()const {return m_Dir;}; TVector3 GetDir()const {return m_Dir;};
double GetTrackLength()const {return m_TrackLength;}
bool IsFitted() const {return m_isFitted;}
TVector3 GetFurthestPoint()const {return m_FurthestPoint;}
double GetCharge()const {return m_Charge;}
public: //operator public: //operator
// reversed logic so the cluster are ordered from largest to smallest // reversed logic so the cluster are ordered from largest to smallest
...@@ -69,6 +105,8 @@ namespace NPL{ ...@@ -69,6 +105,8 @@ namespace NPL{
friend bool operator<(const LinearCluster3D& A, const LinearCluster3D& B){ friend bool operator<(const LinearCluster3D& A, const LinearCluster3D& B){
return (A.size()>B.size()); return (A.size()>B.size());
} }
ClassDef(LinearCluster3D, 0);
}; };
///////////////////////////////// /////////////////////////////////
class LinearRansac3D{ class LinearRansac3D{
...@@ -109,6 +147,8 @@ namespace NPL{ ...@@ -109,6 +147,8 @@ namespace NPL{
public: public:
std::vector<LinearCluster3D> TreatEvent(std::vector<double>& X, std::vector<double>&Y, std::vector<double>&Z); std::vector<LinearCluster3D> TreatEvent(std::vector<double>& X, std::vector<double>&Y, std::vector<double>&Z);
ClassDef(LinearRansac3D, 0);
}; };
......
NPLib/TrackReconstruction/NPRansacACTAR.cxx 0 → 100644
+ 184
− 0
View file @ 6786c048
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
* *
* Original Author : Damien THISSE contact address: damien.thisse@cea.fr *
* *
* Creation Date : October 2024 *
* Last update : October 2024 *
*---------------------------------------------------------------------------*
* Decription: *
* This class deal with finding all the track event by event *
*****************************************************************************/
#include "NPRansacACTAR.h"
#include "NPLinearRansac3D.h"
using namespace NPL;
using namespace std;
/////////////////////////////////////////////////
RansacACTAR::RansacACTAR()
{
fRANSACMaxIteration = 1000;
fRANSACThreshold = 16;//100;
fRANSACPointThreshold = 0.05;//0.07;
fRANSACChargeThreshold = 0.05;//0.07;
fRANSACDistance = 7;//7;
fMAXBarycenterDistance = 10;
fAngleMax = 2;//degree
fNumberOfPadsX = 128;
fNumberOfPadsY = 128;
Rand = new TRandom3();
}
vector<LinearCluster3D> RansacACTAR::SimpleRansac(vector<double> & vX, vector<double> & vY, vector<double> & vZ, vector<double> & vQ){
fTotalCharge = 0;
fOriginalCloudSize = vX.size();
vector<int> index;
vector<LinearCluster3D> all_clusters;
for(int QQ = 0; QQ < vQ.size(); QQ++){
fTotalCharge+= vQ.at(QQ);
index.push_back(QQ);
}
double RemainingCharge = fTotalCharge;
int loop_counter = 0;
// cout << "Total charge is " << fTotalCharge << " ---> Threshold to end RANSAC is thus " << fTotalCharge*fRANSACChargeThreshold << endl;
// cout << "Total voxels is " << index.size() << " ---> Threshold to end RANSAC is thus " << fOriginalCloudSize*fRANSACPointThreshold << endl;
do{
loop_counter++;
//cout << "Loop " << loop_counter << ": remaining charge is " << RemainingCharge << "; number of points remaining : " << index.size() << endl;
LinearCluster3D foundCluster(&vX, &vY, &vZ, &vQ);
int highest_inliers = 0;
vector<int> index_to_remove;
for(int i = 0; i < fRANSACMaxIteration; i++){
LinearCluster3D aCluster(&vX, &vY, &vZ, &vQ);
vector<int> pos_index;
UInt_t rand1 = Rand->Integer(index.size());
UInt_t rand2;
do{
rand2 = Rand->Integer(index.size());
} while(rand2 == rand1);
TVector3 pt1(vX[index[rand1]], vY[index[rand1]], vZ[index[rand1]]);
TVector3 pt2(vX[index[rand2]], vY[index[rand2]], vZ[index[rand2]]);
for(int j = 0; j < index.size(); j++){
TVector3 pt(vX[index[j]], vY[index[j]], vZ[index[j]]);
if(MinimumDistancePointLine(pt1, pt2, pt) < fRANSACDistance){
aCluster.AddIndex(index[j]);
pos_index.push_back(j);
}
}
if(aCluster.size() > highest_inliers){
foundCluster = aCluster;
highest_inliers = aCluster.size();
index_to_remove = pos_index;
}
}
if(highest_inliers > fRANSACThreshold){
all_clusters.push_back(foundCluster);
for(int i = index_to_remove.size()-1; i >-1; i--){
RemainingCharge-= vQ.at(index[index_to_remove[i]]);
index.erase(index.begin()+index_to_remove.at(i));
}
// RemainingCharge = 0;
// for(int QQ = 0; QQ < index.size(); QQ++){
// RemainingCharge+= vQ.at(index.at(QQ));
// }
}
else break;
} while(RemainingCharge > fTotalCharge*fRANSACChargeThreshold && index.size() > fOriginalCloudSize*fRANSACPointThreshold);
return all_clusters;
}
void RansacACTAR::ReadParameterValue(string filename)
{
bool ReadingStatus = false;
ifstream RansacParamFile;
RansacParamFile.open(filename.c_str());
if(!RansacParamFile.is_open()){
cout << "No Paramter File found for Ransac parameters" << endl;
cout << "Using the default parameter:" << endl;
cout << "RANSAC MaxIteration= " << fRANSACMaxIteration << endl;
cout << "RANSAC Threshold=" << fRANSACThreshold << endl;
cout << "RANSAC ChargeThreshold= " << fRANSACChargeThreshold << endl;
cout << "RANSAC Distance= " << fRANSACDistance << endl;
cout << "RANSAC PointThreshold= " << fRANSACPointThreshold << endl;
cout << "MAX Barycenter Distance= " << fMAXBarycenterDistance << endl;
cout << "Angle Max to merge tracks= " << fAngleMax << endl;
}
else{
string LineBuffer, whatToDo, DataBuffer;
while(!RansacParamFile.eof()){
getline(RansacParamFile,LineBuffer);
string name = "ConfigRansac";
if(LineBuffer.compare(0, name.length(), name) == 0){
cout << endl;
cout << "**** ConfigRansac found" << endl;
ReadingStatus = true;
}
while(ReadingStatus){
whatToDo="";
RansacParamFile >> whatToDo;
// Search for comment symbol (%)
if (whatToDo.compare(0, 1, "%") == 0) {
RansacParamFile.ignore(numeric_limits<streamsize>::max(), '\n' );
}
else if (whatToDo.compare(0,19,"RANSACMaxIteration=") == 0) {
RansacParamFile >> DataBuffer;
fRANSACMaxIteration = atoi(DataBuffer.c_str());
cout << "/// RANSAC MaxIteration= " << " " << fRANSACMaxIteration << " ///" << endl;
}
else if (whatToDo.compare(0,15,"RANSACDistance=") == 0) {
RansacParamFile >> DataBuffer;
fRANSACDistance = atoi(DataBuffer.c_str());
cout << "/// RANSAC Distance= " << " " << fRANSACDistance << " ///" << endl;
}
else if (whatToDo.compare(0,22,"RANSACChargeThreshold=") == 0) {
RansacParamFile >> DataBuffer;
fRANSACChargeThreshold = atof(DataBuffer.c_str());
cout << "/// RANSAC ChargeThreshold= " << " " << fRANSACChargeThreshold << " ///" << endl;
}
else if (whatToDo.compare(0,21,"RANSACPointThreshold=") == 0) {
RansacParamFile >> DataBuffer;
fRANSACPointThreshold = atof(DataBuffer.c_str());
cout << "/// RANSAC PointThreshold= " << " " << fRANSACPointThreshold << " ///" << endl;
}
else if (whatToDo.compare(0,16,"RANSACThreshold=") == 0) {
RansacParamFile >> DataBuffer;
fRANSACThreshold = atoi(DataBuffer.c_str());
cout << "/// RANSAC Threshold= " << " " << fRANSACThreshold << " ///" << endl;
}
else if (whatToDo.compare(0,22,"MAXBarycenterDistance=") == 0) {
RansacParamFile >> DataBuffer;
fMAXBarycenterDistance = atoi(DataBuffer.c_str());
cout << "/// Max Barycenter Distance= " << " " << fMAXBarycenterDistance << " ///" << endl;
}
else if (whatToDo.compare(0,16,"AngleMaxToMerge=") == 0) {
RansacParamFile >> DataBuffer;
fAngleMax = atoi(DataBuffer.c_str());
cout << "/// Angle Max to merge tracks= " << " " << fAngleMax << " ///" << endl;
}
else{
ReadingStatus = false;
}
}
}
}
}
\ No newline at end of file
NPLib/TrackReconstruction/NPRansacACTAR.h 0 → 100644
+ 79
− 0
View file @ 6786c048
#ifndef __RANSACACTAR__
#define __RANSACACTAR__
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
* *
* Original Author : Damien THISSE contact address: damien.thisse@cea.fr *
* *
* Creation Date : October 2024 *
* Last update : October 2024 *
*---------------------------------------------------------------------------*
* Decription: *
* This class deal with finding all the track event by event *
*****************************************************************************/
//C++ Header
#include <stdio.h>
#include <vector>
#include <sstream>
#include <iostream>
#include <fstream>
#include <cmath>
#include <stdlib.h>
#include <limits>
#include <string>
#include <algorithm>
//NPL
#include "NPTrack.h"
#include "NPLinearRansac3D.h"
using namespace NPL;
// ROOT Headers
#include <TH2F.h>
#include <TMath.h>
#include <TVector3.h>
#include <TRandom3.h>
#include <TVector.h>
using namespace std;
namespace NPL{
class RansacACTAR
{
public:
RansacACTAR();
~RansacACTAR() {};
public:
void ReadParameterValue(string filename);
vector<NPL::LinearCluster3D> SimpleRansac(vector<double> & vX, vector<double> & vY, vector<double> & vZ, vector<double> & vQ);
private:
TRandom3* Rand;
private:
float fRANSACThreshold;
float fRANSACPointThreshold;
float fRANSACChargeThreshold;
float fRANSACDistance;
float fRANSACMaxIteration;
float fMAXBarycenterDistance;
float fAngleMax;
int fNumberOfTracksMax;
int fOriginalCloudSize;
double fTotalCharge;
int fNumberOfPadsX;
int fNumberOfPadsY;
};
}
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment