Belle II Software development
TriangularInterpolation Class Reference

The TriangularInterpolation class. More...

Public Member Functions

const vector< xy_t > & getPoints () const
 returns list of verticies
 
const vector< triangle_t > & getTriangles () const
 returns list of triangles
 
 TriangularInterpolation ()=default
 Default constructor.
 
 TriangularInterpolation (vector< xy_t > &pc, vector< triangle_t > &ts, double d)
 More complex constructor.
 
 ~TriangularInterpolation ()=default
 Destructor.
 
void init (vector< xy_t > &points, vector< triangle_t > &triangles, double d)
 Calculate extents of a triangular mesh and build spatial index.
 
xy_t triangleCenter (const triangle_t &p) const
 Calculate triangle center.
 
void makeIndex (int nx, double xmin, double xmax, int ny, double ymin, double ymax)
 Make spatial index.
 
short int sideCross (short int prev, short int curr, const xy_t &r, const xy_t &v0) const
 Determine which triangle side is crossed by a line segment defined by r and v0 points.
 
void weights (short int i, const xy_t &r, double &w0, double &w1, double &w2) const
 Calculate barycentric coordinates of a point inside triangle.
 
short int findTriangle (const xy_t &r0) const
 Find the triangle which contain the point.
 

Protected Attributes

vector< triangle_tm_triangles
 Triangle list.
 
vector< xy_tm_points
 Vertex list.
 
vector< xy_tm_triangleCenters
 Triangle centers.
 
vector< double > m_triangleAreas
 Triangle areas.
 
vector< short int > m_spatialIndex
 Spatial index.
 
double m_xmin
 Border of the region where the spatial index is constructed.
 
double m_xmax
 Border of the region where the spatial index is constructed.
 
double m_ymin
 Border of the region where the spatial index is constructed.
 
double m_ymax
 Border of the region where the spatial index is constructed.
 
unsigned int m_nx
 Spatial index grid size.
 
unsigned int m_ny
 Spatial index grid size.
 
double m_ixnorm {1}
 Reciprocals to speedup the index calculation.
 
double m_iynorm {1}
 Reciprocals to speedup the index calculation.
 

Detailed Description

The TriangularInterpolation class.

This class travers triangular meshes which satisfies the Delaunay condition. In other meshes it may give a wrong result. The mesh is represented by a list of triangles coupled to a list of xy-points. To speed up the traverse a sort of spatial index is constructed where on a regular cartesian grid the closeses triangle center is memorized.

Definition at line 61 of file BFieldComponentBeamline.cc.

Constructor & Destructor Documentation

◆ TriangularInterpolation()

TriangularInterpolation ( vector< xy_t > &  pc,
vector< triangle_t > &  ts,
double  d 
)
inline

More complex constructor.

Definition at line 72 of file BFieldComponentBeamline.cc.

73 {
74 init(pc, ts, d);
75 }
void init(vector< xy_t > &points, vector< triangle_t > &triangles, double d)
Calculate extents of a triangular mesh and build spatial index.

Member Function Documentation

◆ findTriangle()

short int findTriangle ( const xy_t r0) const
inline

Find the triangle which contain the point.

If not returns the closest one. First using the spatial index locate triangle close to the point and then traverse the mesh using Delaunay triangulation properties.

Parameters
r02d cartesian point
Returns
Triangle index in the list

Definition at line 232 of file BFieldComponentBeamline.cc.

233 {
234 unsigned int ix = lrint((r0.x - m_xmin) * m_ixnorm), iy = lrint((r0.y - m_ymin) * m_iynorm);
235 short int curr = (ix < m_nx && iy < m_ny) ? m_spatialIndex[iy + m_ny * ix] : 0;
236 xy_t r = m_triangleCenters[curr];
237 const short int end = m_triangles.size();
238 short int prev = end;
239 do {
240 short int next = sideCross(prev, curr, r, r0);
241 if (next == end) break;
242 prev = curr;
243 curr = next;
244 } while (true);
245 return curr;
246 }
vector< short int > m_spatialIndex
Spatial index.
double m_xmin
Border of the region where the spatial index is constructed.
unsigned int m_ny
Spatial index grid size.
double m_ixnorm
Reciprocals to speedup the index calculation.
vector< triangle_t > m_triangles
Triangle list.
double m_ymin
Border of the region where the spatial index is constructed.
vector< xy_t > m_triangleCenters
Triangle centers.
short int sideCross(short int prev, short int curr, const xy_t &r, const xy_t &v0) const
Determine which triangle side is crossed by a line segment defined by r and v0 points.
unsigned int m_nx
Spatial index grid size.
double m_iynorm
Reciprocals to speedup the index calculation.

◆ getPoints()

const vector< xy_t > & getPoints ( ) const
inline

returns list of verticies

Definition at line 64 of file BFieldComponentBeamline.cc.

64{ return m_points;}

◆ getTriangles()

const vector< triangle_t > & getTriangles ( ) const
inline

returns list of triangles

Definition at line 66 of file BFieldComponentBeamline.cc.

66{ return m_triangles;}

◆ init()

void init ( vector< xy_t > &  points,
vector< triangle_t > &  triangles,
double  d 
)
inline

Calculate extents of a triangular mesh and build spatial index.

Moves vector contents inside the class.

Parameters
pointsList of verticies
trianglesList of triangles
dHint how close spatial index should be built

Definition at line 88 of file BFieldComponentBeamline.cc.

89 {
90 std::swap(points, m_points);
91 std::swap(triangles, m_triangles);
92 const double inf = numeric_limits<double>::infinity();
93 double xmin = inf, xmax = -inf;
94 double ymin = inf, ymax = -inf;
95 auto limit = [&xmin, &xmax, &ymin, &ymax](const xy_t & p) -> void {
96 xmin = std::min(xmin, p.x); xmax = std::max(xmax, p.x);
97 ymin = std::min(ymin, p.y); ymax = std::max(ymax, p.y);
98 };
99 for (const auto& t : m_triangles) {
100 limit(m_points[t.j0]);
101 limit(m_points[t.j1]);
102 limit(m_points[t.j2]);
103 }
104 double xw = xmax - xmin;
105 double yw = ymax - ymin;
106
107 xmin -= xw * (1. / 256), xmax += xw * (1. / 256);
108 ymin -= yw * (1. / 256), ymax += yw * (1. / 256);
109 int nx = lrint(xw * (1 + 1. / 128) / d), ny = lrint(yw * (1 + 1. / 128) / d);
110 nx = std::max(nx, 2);
111 ny = std::max(ny, 2);
112 makeIndex(nx, xmin, xmax, ny, ymin, ymax);
113 }
void makeIndex(int nx, double xmin, double xmax, int ny, double ymin, double ymax)
Make spatial index.

◆ makeIndex()

void makeIndex ( int  nx,
double  xmin,
double  xmax,
int  ny,
double  ymin,
double  ymax 
)
inline

Make spatial index.

Parameters
nxGrid size in X axis
xminLow border on in X
xmaxUpper border on in X
nyGrid size in Y axis
yminLow border on in Y
ymaxUpper border on in Y

Definition at line 136 of file BFieldComponentBeamline.cc.

137 {
138 m_nx = nx, m_ny = ny, m_xmin = xmin, m_xmax = xmax, m_ymin = ymin, m_ymax = ymax;
139 m_ixnorm = (nx - 1) / (xmax - xmin), m_iynorm = (ny - 1) / (ymax - ymin);
140 double dx = (xmax - xmin) / (nx - 1), dy = (ymax - ymin) / (ny - 1);
141 // cout<<nx<<" "<<xmin<<" "<<xmax<<" "<<ny<<" "<<ymin<<" "<<ymax<<endl;
142 m_spatialIndex.resize(nx * ny);
143 m_triangleCenters.resize(m_triangles.size());
144 for (unsigned int i = 0; i < m_triangles.size(); i++) m_triangleCenters[i] = triangleCenter(m_triangles[i]);
145
146 m_triangleAreas.resize(m_triangles.size());
147 auto getTriangleArea = [this](const triangle_t& p) ->double{
148 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
149 double d21x = p2.x - p1.x, d21y = p2.y - p1.y;
150 double d01x = p0.x - p1.x, d01y = p0.y - p1.y;
151 return 1 / (d21x * d01y - d21y * d01x);
152 };
153 for (unsigned int i = 0; i < m_triangles.size(); i++) m_triangleAreas[i] = getTriangleArea(m_triangles[i]);
154
155
156 for (int ix = 0; ix < nx; ix++) {
157 double x = xmin + ix * dx;
158 for (int iy = 0; iy < ny; iy++) {
159 double y = ymin + iy * dy;
160 int imin = -1; double dmin = 1e100;
161 for (unsigned int i = 0; i < m_triangleCenters.size(); i++) {
162 const xy_t& p = m_triangleCenters[i];
163 double d = pow(p.x - x, 2) + pow(p.y - y, 2);
164 if (d < dmin) {imin = i; dmin = d;}
165 }
166 int k = iy + ny * ix;
167 m_spatialIndex[k] = imin;
168 }
169 }
170 }
xy_t triangleCenter(const triangle_t &p) const
Calculate triangle center.
double m_xmax
Border of the region where the spatial index is constructed.
vector< double > m_triangleAreas
Triangle areas.
double m_ymax
Border of the region where the spatial index is constructed.

◆ sideCross()

short int sideCross ( short int  prev,
short int  curr,
const xy_t r,
const xy_t v0 
) const
inline

Determine which triangle side is crossed by a line segment defined by r and v0 points.

Parameters
prevTriangle number which has been already checked
currTriangle number which is being checked
rStarting point of the line segment
v0Ending point of the line segment
Returns
Next triangle index in the list if nothing found returns the total number of triangles in the list

Definition at line 181 of file BFieldComponentBeamline.cc.

182 {
183 const double vx = r.x - v0.x, vy = r.y - v0.y;
184 auto isCrossed = [&vx, &vy, &v0](const xy_t & p1, const xy_t & p0) -> bool {
185 double u0x = p0.x, u0y = p0.y;
186 double ux = p1.x - u0x, uy = p1.y - u0y;
187 double dx = u0x - v0.x, dy = u0y - v0.y;
188 double D = uy * vx - ux * vy;
189 double t = dx * vy - dy * vx;
190 double s = dx * uy - dy * ux;
191 return ((t < D) != (t < 0)) && ((s < D) != (s < 0));
192 };
193
194 const triangle_t& p = m_triangles[curr];
195 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
196 if (p.n0 != prev && isCrossed(p1, p2)) return p.n0;
197 if (p.n1 != prev && isCrossed(p2, p0)) return p.n1;
198 if (p.n2 != prev && isCrossed(p0, p1)) return p.n2;
199 return m_triangles.size();
200 }

◆ triangleCenter()

xy_t triangleCenter ( const triangle_t p) const
inline

Calculate triangle center.

Parameters
pTriangle

Definition at line 120 of file BFieldComponentBeamline.cc.

121 {
122 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
123 return {(p0.x + p1.x + p2.x)* (1. / 3), (p0.y + p1.y + p2.y)* (1. / 3)};
124 }

◆ weights()

void weights ( short int  i,
const xy_t r,
double &  w0,
double &  w1,
double &  w2 
) const
inline

Calculate barycentric coordinates of a point inside triangle.

Parameters
iTriangle index in the list
r2d cartesian point
w0Weight of 0 vertex
w1Weight of 1 vertex
w2Weight of 2 vertex

Definition at line 211 of file BFieldComponentBeamline.cc.

212 {
213 const triangle_t& p = m_triangles[i];
214 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
215 double dx2 = p2.x - r.x, dy2 = p2.y - r.y;
216 double d21x = p2.x - p1.x, d21y = p2.y - p1.y;
217 double d02x = p0.x - p2.x, d02y = p0.y - p2.y;
218 w0 = (dx2 * d21y - dy2 * d21x) * m_triangleAreas[i];
219 w1 = (dx2 * d02y - dy2 * d02x) * m_triangleAreas[i];
220 w2 = 1 - (w0 + w1);
221 }

Member Data Documentation

◆ m_ixnorm

double m_ixnorm {1}
protected

Reciprocals to speedup the index calculation.

Definition at line 271 of file BFieldComponentBeamline.cc.

◆ m_iynorm

double m_iynorm {1}
protected

Reciprocals to speedup the index calculation.

Definition at line 273 of file BFieldComponentBeamline.cc.

◆ m_nx

unsigned int m_nx
protected

Spatial index grid size.

Definition at line 267 of file BFieldComponentBeamline.cc.

◆ m_ny

unsigned int m_ny
protected

Spatial index grid size.

Definition at line 269 of file BFieldComponentBeamline.cc.

◆ m_points

vector<xy_t> m_points
protected

Vertex list.

Definition at line 251 of file BFieldComponentBeamline.cc.

◆ m_spatialIndex

vector<short int> m_spatialIndex
protected

Spatial index.

Definition at line 257 of file BFieldComponentBeamline.cc.

◆ m_triangleAreas

vector<double> m_triangleAreas
protected

Triangle areas.

Definition at line 255 of file BFieldComponentBeamline.cc.

◆ m_triangleCenters

vector<xy_t> m_triangleCenters
protected

Triangle centers.

Definition at line 253 of file BFieldComponentBeamline.cc.

◆ m_triangles

vector<triangle_t> m_triangles
protected

Triangle list.

Definition at line 249 of file BFieldComponentBeamline.cc.

◆ m_xmax

double m_xmax
protected

Border of the region where the spatial index is constructed.

Definition at line 261 of file BFieldComponentBeamline.cc.

◆ m_xmin

double m_xmin
protected

Border of the region where the spatial index is constructed.

Definition at line 259 of file BFieldComponentBeamline.cc.

◆ m_ymax

double m_ymax
protected

Border of the region where the spatial index is constructed.

Definition at line 265 of file BFieldComponentBeamline.cc.

◆ m_ymin

double m_ymin
protected

Border of the region where the spatial index is constructed.

Definition at line 263 of file BFieldComponentBeamline.cc.


The documentation for this class was generated from the following file: