PointManager.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////
2 /// \brief Class to help Slicer4D manage the collection of hits
3 /// \authors mbaird42@indiana.edu
4 /// \date 2013/05/01
5 /////////////////////////////////////////////////////////////////////////
6 
7 #include "Slicer/PointManager.h"
8 
9 #include <cmath>
10 #include <cstdlib>
11 #include <iostream>
12 
13 namespace slicer {
14 
15  //......................................................................
17  { }
18 
19  //......................................................................
21  { }
22 
23  //......................................................................
24  void PointManager::Initalize(double eps, double pepen,
25  double distpen, double opvplpen,
26  bool usepepen)
27  {
28  fEps = eps;
29  fPEPen = pepen;
30  fDistPenP = distpen/6.2;
31  fDistPenC = distpen/4.1;
32  fOpVPlPen = opvplpen;
33  fUsePEPen = usepepen;
34  }
35 
36  //......................................................................
37  void PointManager::Reset(unsigned int n, bool skip)
38  {
39  // initialize cluster IDs to -1 (-1 = "unassigned")
40  fClustIDs.clear();
41  for(unsigned int i = 0; i < n; ++i) fClustIDs.push_back(-1);
42 
43  fNhits = n;
44 
45  fTime .clear();
46  fTres .clear();
47  fPlane.clear();
48  fCell .clear();
49  fPE .clear();
50  fView .clear();
51 
52  fSkip = skip;
53 
54  NScount = 0;
55  RQcount = 0;
56  }
57 
58  //......................................................................
60  {
61 
62  fNeighborList.clear();
63  fNeighborList.resize(fNhits);
64 
65  // First add each point to its own list (the DBSCAN algorithm needs
66  // this to work correctly.)
67  for(unsigned int i = 0; i < fNhits; ++i) fNeighborList[i].push_back(i);
68 
69  // use i+1 < fNhits to protect against the case where fNhits = 0
70  for(unsigned int i = 0; i+1 < fNhits; ++i) {
71  for(unsigned int j = i+1; j < fNhits; ++j) {
72 
73  // If hits are too far apart in time, then we don't need to determine
74  // the neighbor score. Since the hits are time sorted, this will
75  // also be true of all hits later in the list.
76  if( fSkip && ((fTime[j]-fTime[i]) > 1000.0) ) break;
77 
78  double score = this->NeighborScore(i,j);
79  if(score <= fEps) {
80  // i and j are neighbors, so they go into each others lists
81  fNeighborList[i].push_back(j);
82  fNeighborList[j].push_back(i);
83  } // end if
84 
85  } // end for j
86  } // end for i
87 
88  }
89 
90  //......................................................................
92  {
93  // An optional function to clean up the neighbor assignments (to
94  // be run AFTER the FindNeighbors function.)
95  //
96  // Currently, this function will remove points from neighborhoods
97  // iff all of their neighbors are in the opposite view.
98 
99  for(unsigned int i = 0; i < fNhits; ++i) {
100 
101  if(fNeighborList[i].size() <= 1) continue;
102 
103  bool goodneighbors = false;
104 
105  // The first point in the neighbor list is the original point,
106  // so start with the second point in the list.
107  for(unsigned int j = 1; j < fNeighborList[i].size(); ++j) {
108  if(fView[i] == fView[fNeighborList[i][j]]) {
109  goodneighbors = true;
110  break;
111  }
112  } // end for j (loop over all neighbors)
113 
114  if(!goodneighbors) {
115 
116  // erase point i from the neighbor list of all of point i's neighbors
117  for(unsigned int j = 1; j < fNeighborList[i].size(); ++j) {
118 
119  unsigned int n = fNeighborList[i][j];
120  std::vector<unsigned int> nList = fNeighborList[n];
121 
122  fNeighborList[n].clear();
123 
124  for(unsigned int k = 0; k < nList.size(); ++k) {
125 
126  if(nList[k] != i) fNeighborList[n].push_back(nList[k]);
127 
128  } // end for k (loop through the neighbor list of the j-th neighbor of i)
129 
130  } // end for j (loop over all neighbors of i)
131 
132  // clear point i's neighbor list and add point i back to that list
133  fNeighborList[i].clear();
134  fNeighborList[i].push_back(i);
135 
136  } // end if no good neighbors
137 
138  } // end for i (loop over all hits)
139 
140  }
141 
142  //......................................................................
143  std::vector<unsigned int> PointManager::regionQuery(unsigned int point)
144  {
145  RQcount++;
146  return fNeighborList[point];
147  }
148 
149  //......................................................................
150  void PointManager::changeClustID(unsigned int point, unsigned int id)
151  {
152  fClustIDs[point] = id;
153  }
154 
155  //......................................................................
156  void PointManager::AddHit(double t, double tres,
157  unsigned int p, unsigned int c, double pe,
158  geo::View_t v)
159  {
160  fTime .push_back(t);
161  fTres .push_back(tres);
162  fPlane.push_back(p);
163  fCell .push_back(c);
164  fPE .push_back(pe);
165  fView .push_back(v);
166  }
167 
168  //......................................................................
169  double PointManager::NeighborScore(unsigned int i, unsigned int pt)
170  {
171  NScount++;
172 
173  double score = 0.0;
174  double T = 0.0, R = 0.0, Z = 0.0, XY = 0.0, PE = 0.0;
175 
176  double dT = std::abs(fTime[i] - fTime[pt]);
177  double dC = std::max(fCell[i], fCell[pt]) - std::min(fCell[i], fCell[pt]);
178  double dP = std::max(fPlane[i], fPlane[pt]) - std::min(fPlane[i], fPlane[pt]);
179 
180  ///\ TODO: perhaps some of the below calculations can be made more
181  // efficient...
182  if(fView[i] == fView[pt]) {
183  XY = dC/fDistPenC;
184  Z = dP/fDistPenP;
185  R = sqrt(dC*dC*16.81 + dP*dP*38.44);
186  if(fUsePEPen) {
187  double denom = sqrt(fPE[i]*fPE[i]+fPE[pt]*fPE[pt]);
188  if(denom == 0.0) denom = 1.0e-6;
189  PE = pow(fPEPen/denom, 5);
190  }
191  }
192  else {
193  XY = 0.0;
194  Z = dP/fOpVPlPen;
195  R = dP*6.2;
196  if(fUsePEPen) {
197  double denom = sqrt(fPE[i]*fPE[i]+fPE[pt]*fPE[pt]);
198  if(denom == 0.0) denom = 1.0e-6;
199  PE = pow(fPEPen/denom, 5);
200  }
201  }
202 
203  T = (dT - R/(30.0))/sqrt(fTres[i]*fTres[i]+fTres[pt]*fTres[pt]);
204 
205  score = T*T + Z*Z + XY*XY + PE;
206 
207  return score;
208  }
209 
210  //......................................................................
212  {
213  std::cout << "\n\n==========\nDebug Print:\n\n";
214  std::cout << "Number of hits = " << fNhits << "\n";
215  std::cout << "Number of Neighbor Score calculations = " << NScount << "\n";
216  std::cout << "Number of Region Queries = " << RQcount << "\n";
217  std::cout << "\n==========\n\n";
218  }
219 
220  //......................................................................
221  void PointManager::PreparekDistInfo(unsigned int kmax)
222  {
223  // This function calculates the nearest neighbor distances for all
224  // hits. This function is NOT used for clustering, it is only used to
225  // help determine the correct input parameters for Slicer4D.
226 
227  // BTW, this is NOT maximally efficient but it doesn't have to be since
228  // it is only intended to be used to help tune the Slicer4D input
229  // parameters.
230 
231  fkDistList.clear();
232  fkDistList.resize(fNhits);
233 
234  for(unsigned int i = 0; i < fNhits; ++i) {
235 
236  std::vector<double> list;
237  for(unsigned int k = 0; k < kmax; ++k) list.push_back(1.0e9);
238 
239  for(unsigned int j = 0; j < fNhits; ++j) {
240 
241  double dist = this->NeighborScore(i,j);
242 
243  if(dist < list[kmax-1]) {
244  // add this dist to the list, resort, and resize
245  list.push_back(dist);
246 
247  bool done = false;
248  while(!done) {
249  done = true;
250  for(unsigned int k = 1; k <= kmax; ++k) {
251  if(list[k-1] > list[k]) {
252  double temp = list[k];
253  list[k] = list[k-1];
254  list[k-1] = temp;
255  done = false;
256  } // end if
257  } // end for k
258  } // end while !done
259 
260  list.resize(kmax);
261 
262  } // end if
263  } // end for j
264 
265  fkDistList[i] = list;
266 
267  } // end for i
268 
269  }
270 
271  //......................................................................
272  double PointManager::kDist(unsigned int k, unsigned int pt)
273  {
274  // This function returns the k-dist value (distance to the k-th nearest
275  // neighbor) for a point. This function is NOT used for clustering, it
276  // is only used to help determine the correct input parameters for Slicer4D.
277 
278  std::vector<double> list = fkDistList[pt];
279 
280  return list[k];
281  }
282 
283 } // end namespace slicer
284 /////////////////////////////////////////////////////////////////////////
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
T max(const caf::Proxy< T > &a, T b)
double fPEPen
PE penalty.
Definition: PointManager.h:60
void Initalize(double eps, double pepen, double distpen, double opvplpen, bool usepepen)
std::vector< unsigned int > fPE
PE for each hit.
Definition: PointManager.h:52
double NeighborScore(unsigned int i, unsigned int pt)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void Reset(unsigned int n, bool skip=true)
std::vector< double > fTime
Time for each hit.
Definition: PointManager.h:48
unsigned int fNhits
Total number of hits.
Definition: PointManager.h:55
void PreparekDistInfo(unsigned int kmax)
float abs(float number)
Definition: d0nt_math.hpp:39
double dist
Definition: runWimpSim.h:113
Float_t Z
Definition: plot.C:38
bool fUsePEPen
Use PE penalty term?
Definition: PointManager.h:59
std::vector< unsigned int > fPlane
Plane number for each hit.
Definition: PointManager.h:50
base_types push_back(int_type())
double fDistPenP
distance penalty converted to number of planes
Definition: PointManager.h:61
#define R(x)
const double j
Definition: BetheBloch.cxx:29
std::vector< unsigned int > regionQuery(unsigned int point)
void AddHit(double t, double tres, unsigned int p, unsigned int c, double pe, geo::View_t v)
OStream cout
Definition: OStream.cxx:6
double fDistPenC
distance penalty converted to number of cells
Definition: PointManager.h:62
double fOpVPlPen
opposite view plane penalty
Definition: PointManager.h:63
std::vector< double > fTres
T resolution for each hit.
Definition: PointManager.h:49
std::vector< std::vector< double > > fkDistList
A list of the distances to the closest neighbors for each hit (only used for analysis) ...
Definition: PointManager.h:45
double fEps
Neighbor Score cut off value.
Definition: PointManager.h:56
Class to help Slicer4D manage the collection of hits.
double T
Definition: Xdiff_gwt.C:5
T min(const caf::Proxy< T > &a, T b)
unsigned int RQcount
Count of the # of times regionQuery is called (for debugging)
Definition: PointManager.h:66
std::vector< int > fClustIDs
Assigned cluster ID for each hit.
Definition: PointManager.h:43
bool fSkip
Skip calculating the full neighbor score if hits are too far apart in time.
Definition: PointManager.h:57
std::vector< unsigned int > fCell
Cell number for each hit.
Definition: PointManager.h:51
std::vector< std::vector< unsigned int > > fNeighborList
A list of neighbors for each hit.
Definition: PointManager.h:44
std::vector< geo::View_t > fView
View for each hit.
Definition: PointManager.h:53
double kDist(unsigned int k, unsigned int pt)
unsigned int NScount
Count of the # of times NeighborScore is called (for debugging)
Definition: PointManager.h:65
void changeClustID(unsigned int point, unsigned int id)