The TriangularInterpolation class. More...

Public Member Functions
const vector< xy_t > &	getPoints () const
	returns list of vertices

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_t >	m_triangles
	Triangle list.

vector< xy_t >	m_points
	Vertex list.

vector< xy_t >	m_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 62 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 73 of file BFieldComponentBeamline.cc.

    {
      init(pc, ts, d);
    }

Member Function Documentation

◆ findTriangle()

short int findTriangle ( const xy_t & r0 ) const

inlinenodiscard

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

r0	2d cartesian point

Returns: Triangle index in the list

Definition at line 233 of file BFieldComponentBeamline.cc.

    {
      unsigned int ix = lrint((r0.x - m_xmin) * m_ixnorm), iy = lrint((r0.y - m_ymin) * m_iynorm);
      short int curr = (ix < m_nx && iy < m_ny) ? m_spatialIndex[iy + m_ny * ix] : 0;
      xy_t r = m_triangleCenters[curr];
      const short int end = m_triangles.size();
      short int prev = end;
      do {
        short int next = sideCross(prev, curr, r, r0);
        if (next == end) break;
        prev = curr;
        curr = next;
      } while (true);
      return curr;
    }

◆ getPoints()

const vector< xy_t > & getPoints ( ) const

inlinenodiscard

returns list of vertices

Definition at line 65 of file BFieldComponentBeamline.cc.

65{ return m_points;}

◆ getTriangles()

const vector< triangle_t > & getTriangles ( ) const

inlinenodiscard

returns list of triangles

Definition at line 67 of file BFieldComponentBeamline.cc.

67{ 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

points	List of vertices
triangles	List of triangles
d	Hint how close spatial index should be built

Definition at line 89 of file BFieldComponentBeamline.cc.

    {
      std::swap(points, m_points);
      std::swap(triangles, m_triangles);
      const double inf = numeric_limits<double>::infinity();
      double xmin = inf, xmax = -inf;
      double ymin = inf, ymax = -inf;
      auto limit = [&xmin, &xmax, &ymin, &ymax](const xy_t & p) -> void {
        xmin = std::min(xmin, p.x); xmax = std::max(xmax, p.x);
        ymin = std::min(ymin, p.y); ymax = std::max(ymax, p.y);
      };
      for (const auto& t : m_triangles) {
        limit(m_points[t.j0]);
        limit(m_points[t.j1]);
        limit(m_points[t.j2]);
      }
      double xw = xmax - xmin;
      double yw = ymax - ymin;
 
      xmin -= xw * (1. / 256), xmax += xw * (1. / 256);
      ymin -= yw * (1. / 256), ymax += yw * (1. / 256);
      int nx = lrint(xw * (1 + 1. / 128) / d), ny = lrint(yw * (1 + 1. / 128) / d);
      nx = std::max(nx, 2);
      ny = std::max(ny, 2);
      makeIndex(nx, xmin, xmax, ny, ymin, ymax);
    }

◆ makeIndex()

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

inline

Make spatial index.

Parameters

nx	Grid size in X axis
xmin	Low border on in X
xmax	Upper border on in X
ny	Grid size in Y axis
ymin	Low border on in Y
ymax	Upper border on in Y

Definition at line 137 of file BFieldComponentBeamline.cc.

    {
      m_nx = nx, m_ny = ny, m_xmin = xmin, m_xmax = xmax, m_ymin = ymin, m_ymax = ymax;
      m_ixnorm = (nx - 1) / (xmax - xmin), m_iynorm = (ny - 1) / (ymax - ymin);
      double dx = (xmax - xmin) / (nx - 1), dy = (ymax - ymin) / (ny - 1);
      //    cout<<nx<<" "<<xmin<<" "<<xmax<<" "<<ny<<" "<<ymin<<" "<<ymax<<endl;
      m_spatialIndex.resize(nx * ny);
      m_triangleCenters.resize(m_triangles.size());
      for (unsigned int i = 0; i < m_triangles.size(); i++) m_triangleCenters[i] = triangleCenter(m_triangles[i]);
 
      m_triangleAreas.resize(m_triangles.size());
      auto getTriangleArea = [this](const triangle_t& p) ->double{
        const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
        double d21x = p2.x - p1.x, d21y = p2.y - p1.y;
        double d01x = p0.x - p1.x, d01y = p0.y - p1.y;
        return 1 / (d21x * d01y - d21y * d01x);
      };
      for (unsigned int i = 0; i < m_triangles.size(); i++) m_triangleAreas[i] = getTriangleArea(m_triangles[i]);
 
 
      for (int ix = 0; ix < nx; ix++) {
        double x = xmin + ix * dx;
        for (int iy = 0; iy < ny; iy++) {
          double y = ymin + iy * dy;
          int imin = -1; double dmin = 1e100;
          for (unsigned int i = 0; i < m_triangleCenters.size(); i++) {
            const xy_t& p = m_triangleCenters[i];
            double d = square(p.x - x) + square(p.y - y);
            if (d < dmin) {imin = i; dmin = d;}
          }
          int k = iy + ny * ix;
          m_spatialIndex[k] = imin;
        }
      }
    }

◆ sideCross()

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

inlinenodiscard

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

Parameters

prev	Triangle number which has been already checked
curr	Triangle number which is being checked
r	Starting point of the line segment
v0	Ending 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 182 of file BFieldComponentBeamline.cc.

    {
      const double vx = r.x - v0.x, vy = r.y - v0.y;
      auto isCrossed = [&vx, &vy, &v0](const xy_t & p1, const xy_t & p0) -> bool {
        double u0x = p0.x, u0y = p0.y;
        double ux = p1.x - u0x, uy = p1.y - u0y;
        double dx = u0x - v0.x, dy = u0y - v0.y;
        double D = uy * vx - ux * vy;
        double t = dx * vy - dy * vx;
        double s = dx * uy - dy * ux;
        return ((t < D) != (t < 0)) && ((s < D) != (s < 0));
      };
 
      const triangle_t& p = m_triangles[curr];
      const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
      if (p.n0 != prev && isCrossed(p1, p2)) return p.n0;
      if (p.n1 != prev && isCrossed(p2, p0)) return p.n1;
      if (p.n2 != prev && isCrossed(p0, p1)) return p.n2;
      return m_triangles.size();
    }

◆ triangleCenter()

xy_t triangleCenter ( const triangle_t & p ) const

inlinenodiscard

Calculate triangle center.

Parameters

p Triangle

Definition at line 121 of file BFieldComponentBeamline.cc.

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

◆ 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

i	Triangle index in the list
r	2d cartesian point
w0	Weight of 0 vertex
w1	Weight of 1 vertex
w2	Weight of 2 vertex

Definition at line 212 of file BFieldComponentBeamline.cc.

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