Utilities Namespace Reference

More...

Namespaces
namespace	Geometry
	Namespace for geometrical functions mostly having to do with spherical coordinates.
namespace	IO
	namespace for input/output utilities

Classes
class	D2Matrix
	2 dimensional matrix More...
class	D3Matrix
	3 dimensional matrix, fixed size More...
class	DataFrame
	class for impoting data from a csv file and allowing label string lookup like a data frame. More...
class	HilbertCurve
	Class for calculating the Hilbert curve distance in two dimensions. More...
class	Interpolator
	Bilinear interpolation class for interpolating from a 2D uniform grid. More...
class	MixedVector
	A container that can hold mixed objects all derived from a base class and retains the ability to access derived class functions/members. More...
class	MixedVector< BaseT * >
	A MixedVector for pointers. More...
class	RandomNumbers_NR
	This is a class for generating random numbers. It simplifies and fool proofs initialization and allows for multiple independent series of numbers. More...
class	ShuffledIndex
	Gives a randomized sequence of numbers from 0 to N-1. More...
class	SUMMER
	class for adding large amounts of numbers with less error than the simple sum More...
class	SymmetricMatrix
	Symetric matrix. More...

Functions
void	print_date ()
template<class T >
void	fill_linear (std::vector< T > &v, size_t n, T min, T max)
template<class T >
void	fill_logarithmic (std::vector< T > &v, size_t n, T min, T max)
template<class T >
long	locate (const std::vector< T > &v, const T x)
template<class T , class F >
long	locate (const std::vector< T > &v, F x, std::function< bool(F, const T &)> less_than)
template<class T >
void	locate (T *v, unsigned long n, T x, unsigned long *index)
template<typename T >
size_t	locate (const std::vector< T > &v, const std::vector< size_t > &sorted_index, T value, size_t &rank)
	Finds the element of a vector given a sorted index of the vector.
template<class T >
size_t	closest (const std::vector< T > &v, const T x)
	Returns the index of the element of v that is closest to x. v must be sorted.
PosType	InterpolateYvec (std::vector< PosType > &x, std::vector< PosType > &y, PosType xi)
PosType	arctanh (PosType x)
PosType	fmini (PosType a, PosType b)
PosType	fmaxi (PosType a, PosType b)
template<typename T >
T	median (std::vector< T > vec)
	find median of vector
template<typename T >
void	range (std::vector< T > vec, T &max, T &min)
	find the maximum and minimum of a vector
template<typename T >
void	polintT (T xa[], T ya[], int n, T x, T *y, T *dy)
	interpolation
template<typename FunctorType , typename T = double>
T	trapz (FunctorType &func, T a, T b, int n, T *s2)
	used in trapizoidal integral
template<typename FunctorType , typename T = double>
T	nintegrateF (FunctorType func, T a, T b, T tols)
long	lmin (long a, long b)
long	lmax (long a, long b)
void	rotation (float *xout, float *xin, PosType theta)
	Rotates 2 dimensional point without changing input point.
void	rotation (PosType *xout, PosType const *xin, PosType theta)
	Rotates 2 dimensional point without changing input point.
PosType	RandomFromTable (PosType *table, unsigned long Ntable, long *seed)
	Generates a random deviates drawn fEinstein_rom appEinstein_roximately the same as the values of table.
template<typename T , typename P >
Point_2d	subtract (T &p1, P &p2)
template<typename P >
double	crossD (P &O, P &A, P &B)
template<typename T >
size_t	RemoveIntersections (std::vector< T > &curve)
	removes the intersections of the curve
std::vector< Point_2d >	TighterHull (const std::vector< Point_2d > &v)
std::vector< Point_2d >	TightestHull (const std::vector< Point_2d > &v)
	Finds a concave envolope for an arbitrary closed curve. This is done by gridding and then finding points that are withing a sertain distance of a segment of the curve. The outer bounding curve is found and then the cuve is shrunck to the closest point on a segment. This should be fool proof, but is relatively slow and might clip some points.
Point_2d	RandomInTriangle (const Point_2d &x1, const Point_2d &x2, const Point_2d &x3, Utilities::RandomNumbers_NR &ran)
	returns random point within a trinagle
Point_2d	RandomInConvexPoly (const std::vector< Point_2d > &pp, Utilities::RandomNumbers_NR &ran)
	return a point within a convex polygon
std::vector< Point_2d >	RandomInConvexPoly (const std::vector< Point_2d > &pp, int N, Utilities::RandomNumbers_NR &ran)
Point_2d	RandomInPoly (std::vector< Point_2d > &pp, Utilities::RandomNumbers_NR &ran)
	return a point within a polygon that doesn't need to be convex
std::vector< Point_2d >	RandomInPoly (std::vector< Point_2d > &pp, int N, Utilities::RandomNumbers_NR &ran)
Point_2d	RandomNearPoly (std::vector< Point_2d > &pp, double R, Utilities::RandomNumbers_NR &ran)
	return a point within distance R of a polygon, i.e. the center of a cicle that intersects the interior of a polygon
std::vector< Point_2d >	RandomNearPoly (std::vector< Point_2d > &pp, int N, double R, Utilities::RandomNumbers_NR &ran)
template<typename P >
void	find_boundaries (std::vector< bool > &bitmap, long nx, std::vector< std::vector< P > > &points, std::vector< bool > &hits_edge, bool add_to_vector=false, bool outer_only=false)
	finds ordered boundaries to regions where bitmap == true
void	find_islands (std::vector< bool > &bitmap, long nx, std::vector< std::vector< long > > &indexes, std::vector< bool > &hits_edge, bool add_to_vector=false)
	find the indexes of areas with bitmap[]=true broken up into diconnected islands
double	interior_mass (const std::vector< Point_2d > &alpha, const std::vector< Point_2d > &x)
	this returns area within the curve x average kappa iwithin the curve
double	interior_mass (const std::vector< RAY > &rays)
	this returns area within the curve x average kappa iwithin the curve
template<typename T >
std::vector< T >	convex_hull (const std::vector< T > &PP)
	Returns a vector of points on the convex hull in counter-clockwise order.
template<typename T >
void	convex_hull (const std::vector< T > &P, std::vector< size_t > &hull)
template<typename T >
void	concave (std::vector< T > &init_points, std::vector< T > &hull_out, double scale)
	Creates the concave hull of a group of 2 dimensional points by the shrink-wrap algorithm.
template<typename T >
std::vector< T >	concave2 (std::vector< T > &init_points, double scale)
template<typename Ptype >
double	distance_to_segment (const Ptype &P, const Ptype &S1, const Ptype &S2)
template<typename Ptype >
double	distance_to_segment (const Ptype &P, const Ptype &S1, const Ptype &S2, Ptype &closest_point)
template<typename Ptype >
bool	segments_cross (const Ptype &a1, const Ptype &a2, const Ptype &b1, const Ptype &b2)
template<typename Ptype >
bool	inCurve (const Ptype &x, const std::vector< Ptype > &H)
template<typename Ptype >
bool	inhull (PosType x[], const std::vector< Ptype > &H)
template<>
bool	inhull< Point * > (PosType x[], const std::vector< Point * > &H)
template<>
bool	inhull< RAY > (PosType x[], const std::vector< RAY > &H)
	finds in x is within the curve discribed by the H[].x points ie image points
template<>
bool	inhull< PosType * > (PosType x[], const std::vector< PosType * > &H)
template<typename Ptype >
std::vector< Ptype >	concaveK (std::vector< Ptype > &points, int &k, bool check=true)
template<typename Ptype >
void	testconcaveK ()
Point_2d	line_intersection (const Point_2d &v1, const Point_2d &v2, const Point_2d &w1, const Point_2d &w2)
	return the interesetion point of line defined by v1,v2 and w1,w2
bool	circleIntersetsCurve (const Point_2d &x, double r, const std::vector< Point_2d > &v)
	returns true if a circle of radius r around the point x intersects with the curve v. Does not include case where one compleatly encloses the other.
bool	circleOverlapsCurve (const Point_2d &x, double r, const std::vector< Point_2d > &v)
	Returns true if there is any overlap between a curve and a circle. This includea cases where one compleatly encloses the other.
std::vector< Point_2d >	envelope (const std::vector< Point_2d > &v, const std::vector< Point_2d > &w)
	Find a curve that is made up of segments from v and w that surrounds them and does not self intersect.
std::vector< Point_2d >	envelope2 (const std::vector< Point_2d > &v, const std::vector< Point_2d > &w)
	Find a curve that is made up of segments from v and w that surrounds them and does not self intersect.
std::vector< std::vector< Point_2d > >	thicken_poly (const std::vector< Point_2d > &v, double R)
	return the boundaries of the region that is within R of the curve v
void	ReadFileNames (std::string dir, const std::string filespec, std::vector< std::string > &filenames, const std::string file_non_spec=" ", bool verbose=false)
	Reads all the fits files in a directory into a vector of PixelMaps.
PosType	sepSQR (PosType *xx, PosType *yy)
	Separation squared between two positions in 2 dimensions.
void	double_sort (unsigned long n, PosType *arr, unsigned long *brr)
void	double_sort_points (unsigned long n, PosType *arr, Point *brr)
	Sorts points in a point array.
void	quicksortPoints (Point *pointarray, PosType *arr, unsigned long N)
	Sorts points from smallest to largest according to the value of arr[]. Sorts arr[] and pointarray[] simultaneously.
void	quicksortPoints (Point *pointarray, double(*func)(Point &), unsigned long N)
template<typename D >
void	quicksort (unsigned long *particles, D *arr, unsigned long N)
template<typename T >
void	quickPartition (T pivotvalue, unsigned long *pivotindex, unsigned long *particles, T *arr, unsigned long N)
template<typename T >
void	quickPartitionPoints (T pivotvalue, unsigned long *pivotindex, Point *pointarray, T *arr, unsigned long N)
template<typename T >
void	quickPartitionPoints (T pivotvalue, unsigned long *pivotindex, Point *pointarray, T(*func)(Point &p), unsigned long N)
void	log_polar_grid (Point *i_points, PosType rmax, PosType rmin, PosType *center, long Ngrid)
void	findarea (ImageInfo *imageinfo)
void	writeCurves (int m, ImageInfo *critical, int Ncrit, int index)
PosType	cross (const Point *O, const Point *A, const Point *B)
bool	xorder (Point *p1, Point *p2)
bool	yorder (Point *p1, Point *p2)
PosType	crossD (const double *O, const double *A, const double *B)
PosType	crossD (const Point_2d &O, const Point_2d &A, const Point_2d &B)
std::vector< Point * >	convex_hull (std::vector< Point * > &P)
	Returns a vector of points on the convex hull in counter-clockwise order.
std::vector< double * >	convex_hull (std::vector< double * > &P)
	Returns a vector of points on the convex hull in counter-clockwise order.
std::vector< Point * >	concave_hull (std::vector< Point * > &P, int k, bool test=false)
	Returns a vector of points on the convcave hull in counter-clockwise order.
std::vector< double * >	concave_hull (std::vector< double * > &P, int k)
	Returns a vector of points on the convcave hull in counter-clockwise order.
void	contour_ellipse (std::vector< Point_2d > &P, Point_2d center, unsigned long Npoints, std::vector< Point_2d > &C, double *ellipticity, double *ellipse_area)
	Returns axis ratio, area and points of an ellipse engulfed by some contour (e.g. a contour of same convergence calculated with find_contour).
void	contour_ellipse (std::vector< RAY > &P, Point_2d center, unsigned long Npoints, std::vector< Point_2d > &C, double *ellipticity, double *ellipse_area)
Point_2d	contour_center (std::vector< Point_2d > &P, unsigned long Npoints)
	Returns the center of a contour defined as the midpoint between the two points in the contour that are farthest apart from one another.
Point_2d	contour_center (std::vector< RAY > &P, unsigned long Npoints)
long	IndexFromPosition (PosType *x, long Npixels, PosType range, const PosType *center)
void	PositionFromIndex (unsigned long i, PosType *x, long Npixels, PosType range, PosType const *center)
	This should work for square regions.
long	IndexFromPosition (PosType x, long Npixels, PosType range, PosType center)
PosType	TwoDInterpolator (PosType *x,int Npixels,PosType range,PosType *center,PosType *map,bool init)
	bilinear interpolation from a map.
PosType	TwoDInterpolator (PosType *map)
template<typename T >
int	cutbox (const T *center, T *p1, T *p2, float rmax)
template<int lev>
void	quicksortPoints_multithread (Point *pointarray, PosType *arr, unsigned long N, int level=0)
	Multi-threaded quicksort. The maximum number of threads used is 2^lev. The last parameter should be left out when calling so that it takes the default value.
template<int lev>
void	quicksortPoints_multithread (Point *pointarray, double(*func)(Point &), unsigned long N, int level=0)
	Multi-threaded quicksort. The maximum number of threads used is 2^lev. The function func takes a point and returns the value that is should be sorted by. The last parameter should be left out when calling so that it takes the default value.
template<typename T , int lev>
void	quicksort_multithread (T *array, double(*func)(T &), unsigned long N, int level=0)
	Multi-threaded quicksort. The maximum number of threads used is 2^lev. The function func takes a T type and returns the value that is should be sorted by. The last parameter should be left out when calling so that it takes the default value. std::swap() is used to swap elements of the array.
float	isLeft (Point *p0, Point *p1, PosType *x)
float	isLeft (const Point_2d &p0, const Point_2d &p1, const Point_2d &x)
unsigned long	prevpower (unsigned long k)
int	windings (PosType *x,Point *points,unsigned long Npoints,PosType *area,short image)
	windings(): winding number test for a point in a polygon Returns: Number of times a curves winds around the point x.
int	windings (PosType *x, Point **points, unsigned long Npoints, PosType *area, short image=0)
int	windings (const Point_2d &x, const std::vector< Point_2d > &point, PosType *area)
int	windings (Point_2d &x, std::vector< RAY > &point, PosType *area)
int	windings (PosType *x, Kist< Point > *kist, PosType *area, short image=0)
int	windings2 (PosType *x,Point *points_original,unsigned long Npoints,PosType *area,short image)
	determines whether a point is inside a curve, that has been stretched 1.2 times returns the area of the stretched curve
unsigned long	order_curve4 (Point *curve, long Npoints)
	returns 1 if it is in the curve and 0 if it is out. Borders count as in.
unsigned long	order_curve4 (Kist< Point > *curve)
	Overloads and is dependent on version that takes a point array. Returns number of points that have been ordered. Remaining, unordered points are left at the end of the kist.
unsigned long	order_curve5 (Kist< Point > *curve)
	For odering the curve by the convex hull method. Warning: Does not work very well.
void	ordered_convexhull (Kist< Point > *curve)
	Replaces curve->imagekist with its convex hull. The number of points will change.
void	ordered_shrink_wrap (Kist< Point > *curve)
void	ordered_concavehull (Kist< Point > *curve)
PosType	ConvexHullArea (Kist< Point > *curve)
	Replaces curve with its convex hull. The number of points will change.
template<class T >
T	to_numeric (const std::string &str)
	convert a string to a numerical value of various types
template<>
long	to_numeric< long > (const std::string &str)
template<>
int	to_numeric< int > (const std::string &str)
template<>
float	to_numeric< float > (const std::string &str)
template<>
double	to_numeric< double > (const std::string &str)
template<typename T >
bool	AlwaysTrue (T t)
template<typename T >
bool	AlwaysFalse (T t)
template<typename T >
T	vec_sum (const std::vector< T > &v)
template<class T >
void	Matrix (T **matrix, long rows, long cols)
template<class T >
void	free_Matrix (T **matrix, long rows, long)
PosType **	PosTypeMatrix (size_t rows, size_t cols)
void	free_PosTypeMatrix (PosType **matrix, size_t rows, size_t cols)
PosType **	PosTypeMatrix (long rows1, long rows2, long cols1, long cols2)
void	free_PosTypeMatrix (PosType **matrix, long rows1, long rows2, long cols1, long cols2)
template<class BaseT >
std::size_t	lower_bound (std::vector< BaseT * > &items, PosType target)
template<typename Container >
void	delete_container (Container &c)
	delete the objects that are pointed to in a container of pointers
template<typename T >
T	between (const T &x, const T &l, const T &u)
template<typename T , typename R >
void	shuffle (std::vector< T > &vec, R &ran)
	Shuffles a vector into a random order.
template<typename T >
void	sort_indexes (const std::vector< T > &v, std::vector< size_t > &index)
	Find the indexes that sort a vector in asending order.
template<typename T >
void	sort_indexes (const T *v, std::vector< size_t > &index, size_t N)
template<typename T >
void	sort_indexes_decending (const std::vector< T > &v, std::vector< size_t > &index)
	Find the indexes that sort a vector in descending order.
template<typename T >
void	apply_permutation (T *vec, const std::vector< std::size_t > &p)
template<typename T >
void	apply_permutation (std::vector< T > &vec, const std::vector< std::size_t > &p)
void	powerspectrum2d (std::valarray< double > const &aa, std::valarray< double > const &bb, long nx, long ny, double boxlx, double boxly, std::vector< double > &ll, std::vector< double > &Pl, double zeropaddingfactor)
	Calculates power spectrum from a 2d map or the cross-power spectrum between two 2d maps.
void	powerspectrum2d (std::valarray< double > &aa, long nx, long ny, double boxlx, double boxly, std::vector< double > &ll, std::vector< double > &Pl)
void	powerspectrum2dprebin (std::valarray< double > &aa, long nx, long ny, double boxlx, double boxly, const std::vector< double > &ll, std::vector< double > &Pl, std::vector< double > &llave)
void	powerspectrum2d (std::valarray< float > const &aa, std::valarray< float > const &bb, long nx, long ny, double boxlx, double boxly, std::vector< double > &ll, std::vector< double > &Pl, double zeropaddingfactor)
void	powerspectrum2d (std::valarray< float > &aa, long nx, long ny, double boxlx, double boxly, std::vector< double > &ll, std::vector< double > &Pl)
void	powerspectrum2dprebin (std::valarray< float > &aa, int nx, int ny, double boxlx, double boxly, const std::vector< double > &ll, std::vector< double > &Pl, std::vector< double > &llave)
int	GetNThreads ()
	returns the compiler variable N_THREADS that is maximum number of threads to be used.
void	splitstring (std::string &line, std::vector< std::string > &vec, const std::string &delimiter)
	split string into vector of seporate strings that were seporated by
template<typename T >
double	KleinSum (std::vector< T > &input)
template<typename T , typename F >
double	PairWiseSum (T *begin, T *end, F value)
template<typename T >
double	PairWiseSum (T *begin, T *end)
	Does a parwise sumation of a vector which incresses the precision of the sumation to ~O(epsilon log(n) )
template<typename T >
double	PairWiseSum (std::vector< T > &v)
	Does a parwise sumation of a vector which incresses the precision of the sumation to ~O(epsilon log(n) )
template<typename T , typename F >
double	PairWiseSum (std::vector< T > &v, F value)
	This version allows you to specify a function( value(T * p) ), that returns the value of a pointer to a T type that is to be summed. This is useful for summing the squares or summing some a particular variable within a list of objects. A lambda function is particularly useful here.
template<typename T = double>
double	hypergeometric (T a, T b, T c, T x)
template<typename T >
std::valarray< T >	AdaptiveSmooth (const std::valarray< T > &map_in, size_t Nx, size_t Ny, T value)
	Smooth a 2 dimensional map stored in a valarray with a density dependent kernel.
bool	xorderD (double *p1, double *p2)
bool	yorderD (double *p1, double *p2)
long	IndexFromPosition (PosType *x, long Nx, long Ny, PosType range, const PosType *center)
	this is the nonsquare version of the function, it will return -1 is outside of region
void	PositionFromIndex (unsigned long i, PosType *x, long Nx, long Ny, PosType range, PosType const *center)
	This should work for square or rectangular regions as long as Npixels and range are the x-axis values and the pixels are square.

Variables
const double	nXbin =64.

Detailed Description

The Utilities namespace contains functions for wide use in many classes that perform generic tasks.

structure to construct, hold and destruct image and source points

• \ingEinstein_roup ChangeLens
• \ingEinstein_roup ChangeLens
• orders points in a curve, separates disconnected curves curves[0...Maxcurves] must be allocated before

uses neighbors-of-neighbors to split into curves and then uses sorts by angle and then walks the ist sorting

can break down for crescent curves

The two functions below are inverses of each other for converting between a 1d array index and a square grid of positions Npixels in the number of point is 1 dimension index is between 0 and Npixels*Npixels-1 If x is outside of the region -1 is returned.

Function Documentation

AdaptiveSmooth()

template<typename T >

std::valarray< T > Utilities::AdaptiveSmooth	(	const std::valarray< T > &	map_in,
		size_t	Nx,
		size_t	Ny,
		T	value )

Smooth a 2 dimensional map stored in a valarray with a density dependent kernel.

The smoothing is done by finding the circle around each point whose total pixel values are larger than value. In the case of a density map made from particles if value = (mass of particle)*(number of neighbours) an approximate N nearest neighbour smoothing is done.

concave()

template<typename T >

void Utilities::concave	(	std::vector< T > &	init_points,
		std::vector< T > &	hull_out,
		double	scale )

Creates the concave hull of a group of 2 dimensional points by the shrink-wrap algorithm.

The type of the input vector points must have an operator []. If the input vector is the same as the output vector it will be replaced, and the function will still work.

It is guaranteed that the resulting hull will surround the all the points. Any edge that is greater than scale will be refined until it is either smaller than scale or it cannot be refined further. As a result some edges might be larger than scale and some smaller.

This should be a NlogN algorithm.

The algorithm: 1) The convex hull is found. 2) The longest edge is found 3) all the points that are not in the hull are tested to see if they are within the rays extending from the end point perpendicular to the edge. 4) Of the points that are the one that makes the smallest area triangle with the end points is chosen and added 5) go back to 3 if there are edges that are larger than scale and new points exist to be added 6) remove all intersections in the hull

concave_hull() [1/2]

std::vector< double * > Utilities::concave_hull	(	std::vector< double * > &	P,
		int	k )

Returns a vector of points on the convcave hull in counter-clockwise order.

This uses a K-nearest neighbour adapted from Moreira & Santos (GRAPP 2007 conference proceedings). This is a modified gift wrap algorithm using k neighbours. The value of k will automatically increase when certain special cases are encountered.

This is an overloaded version of the other concave_hull()

concave_hull() [2/2]

std::vector< Point * > Utilities::concave_hull	(	std::vector< Point * > &	P,
		int	k,
		bool	test = false )

Returns a vector of points on the convcave hull in counter-clockwise order.

This uses a K-nearest neighbour adapted from Moreira & Santos (GRAPP 2007 conference proceedings). This is a modified gift wrap algorithm using k neighbours. The value of k will automatically increase when certain special cases are encountered.

concaveK()

template<typename Ptype >

std::vector< Ptype > Utilities::concaveK	(	std::vector< Ptype > &	points,
		int &	k,
		bool	check = true )

will be added back later

contour_center()

Point_2d Utilities::contour_center	(	std::vector< Point_2d > &	P,
		unsigned long	Npoints )

Returns the center of a contour defined as the midpoint between the two points in the contour that are farthest apart from one another.

The performance of the algorithm is ~O(N^2). Less naive methods go like O(N) at best. Most commonly a combined convex hull plus rotating calipers algorithm is used. Since we have the convex_hull already, we only need to implement the latter algorithm.

contour_ellipse()

void Utilities::contour_ellipse	(	std::vector< Point_2d > &	P,
		Point_2d	center,
		unsigned long	Npoints,
		std::vector< Point_2d > &	C,
		double *	ellipticity,
		double *	ellipse_area )

Returns axis ratio, area and points of an ellipse engulfed by some contour (e.g. a contour of same convergence calculated with find_contour).

The axis ratio of the ellipse b/a is equal to the ratio of the distances between center and the nearest contour point (i.e. b) and between center and the farthest contour point (i.e. a). NOTE that the center used to calculate a and b is an input parameter. The definition of the center is crucial to the meaning of above output parameters. The center of the convex_hull produces for even slightly distorted hulls significant offsets resulting in overestimated major axis values (a). The function Utilities::contour_center() calculates the center as the midpoint between the two points in the contour that are farthest apart from one another, which gives already more reliable results. The output vector describing the ellipse is resized to match the size of the contour vector.

convex_hull()

template<typename T >

void Utilities::convex_hull	(	const std::vector< T > &	P,
		std::vector< size_t > &	hull )

Returns a vector of points on the convex hull in counter-clockwise order.

ConvexHullArea()

PosType Utilities::ConvexHullArea

(

Kist< Point > *

curve

)

Replaces curve with its convex hull. The number of points will change.

gives the area within the convex hull of the curve

double_sort_points()

void Utilities::double_sort_points	(	unsigned long	n,
		PosType *	arr,
		Point *	brr )

Sorts points in a point array.

arr array uses NR standard indexing i.e arr[1...n] but brr[0..n-1] if the point array is two-way-coupled to another point array the image pointers of that array will follow sort if the array is not two-way-coupled to another the image pointers in the other array will be untouched

envelope()

std::vector< Point_2d > Utilities::envelope	(	const std::vector< Point_2d > &	v,
		const std::vector< Point_2d > &	w )

Find a curve that is made up of segments from v and w that surrounds them and does not self intersect.

v and w can be non-self intersecting If they do not intersect and one is not inside the other an empty vector is returned.

Unlike Utilities::envelope(), this algorithm is pretty foolproof. It uses the same concept as Utilities::TightestHull(). There may be small segements that are not in either curve, but they should increase the area be a small fraction.

envelope2()

std::vector< Point_2d > Utilities::envelope2	(	const std::vector< Point_2d > &	v,
		const std::vector< Point_2d > &	w )

Find a curve that is made up of segments from v and w that surrounds them and does not self intersect.

v and w must be non-self intersecting If they do not intersect and one is not inside the other an empty vector is returned

fill_linear()

template<class T >

void Utilities::fill_linear	(	std::vector< T > &	v,
		size_t	n,
		T	min,
		T	max )

Fills a vector with equidistant points from [min, max].

fill_logarithmic()

template<class T >

void Utilities::fill_logarithmic	(	std::vector< T > &	v,
		size_t	n,
		T	min,
		T	max )

Fills a vector with logarithmically equidistant points from [min, max].

find_boundaries()

template<typename P >

void Utilities::find_boundaries	(	std::vector< bool > &	bitmap,
		long	nx,
		std::vector< std::vector< P > > &	points,
		std::vector< bool > &	hits_edge,
		bool	add_to_vector = false,
		bool	outer_only = false )

finds ordered boundaries to regions where bitmap == true

This can be used to find critical curves or contours. If the boundary curve touches the edge of the bitmap it will be indicated in hits_boundary as true.

Boundaries will never cross or lead off the grid.
On the edges they will leave the edge pixels out even if they should be in. This is a technical compromise.

Output points are in pixel units with (0,0) being point (0,0)

Parameters

outer_only

finds only the fist boundary which will be the outer one if there are not seporated islands

InterpolateYvec()

PosType Utilities::InterpolateYvec	(	std::vector< PosType > &	x,
		std::vector< PosType > &	y,
		PosType	xi )

Interpolate (cubic interpolation) the value of a function given xi

locate() [1/3]

template<typename T >

size_t Utilities::locate	(	const std::vector< T > &	v,
		const std::vector< size_t > &	sorted_index,
		T	value,
		size_t &	rank )

Finds the element of a vector given a sorted index of the vector.

The vector v does not need to be sorted, but the index does with Utiltites::sort_index(), Utilities::sort_index_decending() or in some other way.

returns the index of v that is the largest value that is <= value

Parameters

v	vector to search, does not need to be sorted
sorted_index	sorted index for v, could be assending or decending
value	value to be matched
rank	index in sorted_index

locate() [2/3]

template<class T >

long Utilities::locate	(	const std::vector< T > &	v,
		const T	x )

Locates the element of the given vector which, together with the following element, brackets the given number. If x is smaller than the smallest entry or larger than the largest, the result is either -1 or n-1.

locate() [3/3]

template<class T , class F >

long Utilities::locate	(	const std::vector< T > &	v,
		F	x,
		std::function< bool(F, const T &)>	less_than )

Locates the element of the given vector which, together with the following element, brackets the given number. If x is smaller than the smallest entry or larger than the largest, the result is either -1 or n-1.

nintegrateF()

template<typename FunctorType , typename T = double>

T Utilities::nintegrateF	(	FunctorType	func,
		T	a,
		T	b,
		T	tols )

Parameters

func	struct or class to be integrated
a	limit of integrations
b	limit of integrations
tols	target fractional error

order_curve4()

unsigned long Utilities::order_curve4	(	Point *	curve,
		long	Npoints )

returns 1 if it is in the curve and 0 if it is out. Borders count as in.

Orders points on a closed curve.

returns 1 if it is in the curve and 0 if it is out. Borders count as in. returns 1 if it is in the curve and 0 if it is out. Borders count as in.

The algorithm first finds the "center" of the curve. It then does a rough ordering according to the angle around this center. It then walks the curve jumping to a neighbor cell each step choosing a neighbor along one of the x or y-axis before taking a diagonal step. If it comes to the point where there is no more neighbors (as may occur after going through a self-intersection and then returning to it) the algorithm backtracks until it finds a point in the ordered list that is also a neighbor to a point in the not yet ordered list and attaches this to the end of the ordered list and continuous to walk. This algorithm works well at finding a closed loop. It can cut off points from the curve that are either in loops or if four cells intersect and are all on the curve as can happen when there is a lot of structure in the curve that is not resolved at the gridsize used. The points that are cut off are at the end of the array in no guaranteed order.

Returns the number of point that have been ordered - total number minus the cuttout points.

This algorithm could be improve by inserting the remaining points, if any, into the existing curve and recursively calling itself.

order_curve5()

unsigned long Utilities::order_curve5

(

Kist< Point > *

curve

)

For odering the curve by the convex hull method. Warning: Does not work very well.

The convex hull is found for the points in the kist. Then each additional point is inserted into the curve where it will increase the length of the curve the least. This method leaves loops where they shouldn't be and probably doesn't handle self-intersections well.

powerspectrum2d() [1/4]

void Utilities::powerspectrum2d	(	std::valarray< double > &	aa,
		long	nx,
		long	ny,
		double	boxlx,
		double	boxly,
		std::vector< double > &	ll,
		std::vector< double > &	Pl )

Parameters

aa	first realspace map to be
nx	number of pixels in x direction
ny	number of pixels in y direction
boxlx	range of image in x direction
boxly	range of image in y direction
ll	output multiplot number of bins
Pl	output binned power spectrum

powerspectrum2d() [2/4]

void Utilities::powerspectrum2d	(	std::valarray< double > const &	aa,
		std::valarray< double > const &	bb,
		long	nx,
		long	ny,
		double	boxlx,
		double	boxly,
		std::vector< double > &	ll,
		std::vector< double > &	Pl,
		double	zeropaddingfactor )

Calculates power spectrum from a 2d map or the cross-power spectrum between two 2d maps.

Adapted from Carlo Giocoli's pl() routine.

Parameters

aa	first realspace map to be
bb	second realspace map, same as aa to get power spectrum
nx	number of pixels in x direction
ny	number of pixels in y direction
boxlx	range of image in x direction
boxly	range of image in y direction
ll	output multiplot number of bins
Pl	output binned power spectrum

powerspectrum2d() [3/4]

void Utilities::powerspectrum2d	(	std::valarray< float > &	aa,
		long	nx,
		long	ny,
		double	boxlx,
		double	boxly,
		std::vector< double > &	ll,
		std::vector< double > &	Pl )

Parameters

aa	first realspace map to be
nx	number of pixels in x direction
ny	number of pixels in y direction
boxlx	range of image in x direction
boxly	range of image in y direction
ll	output multiplot number of bins
Pl	output binned power spectrum

powerspectrum2d() [4/4]

void Utilities::powerspectrum2d	(	std::valarray< float > const &	aa,
		std::valarray< float > const &	bb,
		long	nx,
		long	ny,
		double	boxlx,
		double	boxly,
		std::vector< double > &	ll,
		std::vector< double > &	Pl,
		double	zeropaddingfactor )

Parameters

aa	first realspace map to be
bb	second realspace map, same as aa to get power spectrum
nx	number of pixels in x direction
ny	number of pixels in y direction
boxlx	range of image in x direction
boxly	range of image in y direction
ll	output multiplot number of bins
Pl	output binned power spectrum

powerspectrum2dprebin() [1/2]

void Utilities::powerspectrum2dprebin	(	std::valarray< double > &	aa,
		long	nx,
		long	ny,
		double	boxlx,
		double	boxly,
		const std::vector< double > &	ll,
		std::vector< double > &	Pl,
		std::vector< double > &	llave )

Parameters

aa	first realspace map to be
nx	number of pixels in x direction
ny	number of pixels in y direction
boxlx	range of image in x direction
boxly	range of image in y direction
ll	output multiplot number of bins
Pl	output binned power spectrum
llave	average value of Fourier node in bins

powerspectrum2dprebin() [2/2]

void Utilities::powerspectrum2dprebin	(	std::valarray< float > &	aa,
		int	nx,
		int	ny,
		double	boxlx,
		double	boxly,
		const std::vector< double > &	ll,
		std::vector< double > &	Pl,
		std::vector< double > &	llave )

Parameters

aa	first realspace map to be
nx	number of pixels in x direction
ny	number of pixels in y direction
boxlx	range of image in x direction
boxly	range of image in y direction
ll	output multiplot number of bins
Pl	output binned power spectrum
llave	average value of Fourier node in bins

prevpower()

unsigned long Utilities::prevpower

(

unsigned long

k

)

This function finds the largest power of 2 that is < k

quicksortPoints_multithread()

template<int lev>

void Utilities::quicksortPoints_multithread	(	Point *	pointarray,
		double(*	func )(Point &),
		unsigned long	N,
		int	level = 0 )

Multi-threaded quicksort. The maximum number of threads used is 2^lev. The function func takes a point and returns the value that is should be sorted by. The last parameter should be left out when calling so that it takes the default value.

This function is different from quicksort_multithread() in that it uses SwapPointsInArray() instead of std::swap() which is needed to make the image pointers follow the swap.

RandomFromTable()

PosType Utilities::RandomFromTable	(	PosType *	table,
		unsigned long	Ntable,
		long *	seed )

Generates a random deviates drawn fEinstein_rom appEinstein_roximately the same as the values of table.

\ingEinstein_roup Utill

ReadFileNames()

void Utilities::ReadFileNames	(	std::string	dir,
		const std::string	filespec,
		std::vector< std::string > &	filenames,
		const std::string	file_non_spec = " ",
		bool	verbose = false )

Reads all the fits files in a directory into a vector of PixelMaps.

The input fits files must have .fits in their names in addition to the string filespec.

void Utilities::LoadFitsImages( std::string dir /// path to directory containing fits files ,const std::string& filespec /// string of charactors in fits file name that are matched ,std::vector<PixelMap> & images /// output vector of PixelMaps ,int maxN /// maximum number of images that will be read in ,double resolution /// resolution (rad) of fits image if not given in fits file, use default or -1 otherwise ,bool verbose /// lists files to stdout ){

DIR *dp = opendir( dir.c_str() ); struct dirent *dirp; struct stat filestat; std::string filepath,filename; size_t count = 0;

if (dp == NULL) { throw std::runtime_error("error opening directory"); return; }

while ((dirp = readdir( dp )) && count < maxN) { filepath = dir + "/" + dirp->d_name;

If the file is a directory (or is in some way invalid) we'll skip it if (stat( filepath.c_str(), &filestat )) continue; if (S_ISDIR( filestat.st_mode )) continue;

filename = dirp->d_name; if(filename.find(".fits") != std::string::npos){ if(filename.find(filespec) != std::string::npos){ if(verbose) std::cout << "reading " << filepath << std::endl; PixelMap map(filepath,resolution); images.push_back(std::move(map)); images.push_back(PixelMap(filepath,resolution)); ++count; } } }

closedir( dp );

std::cout << count << " fits files read." << std::endl; return ; }

Reads all the fits files in a directory into a vector of PixelMaps.

The input fits files must have .fits in their names in addition to the string filespec.

void Utilities::LoadFitsImages( std::string dir /// path to directory containing fits files ,std::vector<std::string> filespecs /// string of charactors in fits file name that are matched ,std::vector<std::string> file_non_specs /// string of charactors in fits file name cannot have ,std::vector<PixelMap> & images /// output vector of PixelMaps ,std::vector<std::string> & names /// file names ,int maxN /// maximum number of images that will be read in ,double resolution /// resolution (rad) of fits image if not given in fits file, use default or -1 otherwise ,bool verbose /// lists files to stdout ){

DIR *dp = opendir( dir.c_str() ); struct dirent *dirp; struct stat filestat; std::string filepath,filename; size_t count = 0;

if (dp == NULL) { throw std::runtime_error("error opening directory"); return; }

while ((dirp = readdir( dp )) && count < maxN) { filepath = dir + "/" + dirp->d_name;

If the file is a directory (or is in some way invalid) we'll skip it if (stat( filepath.c_str(), &filestat )) continue; if (S_ISDIR( filestat.st_mode )) continue;

filename = dirp->d_name; if(filename.find(".fits") != std::string::npos){ bool read =true; for(int i=0;i<filespecs.size();++i) if(filename.find(filespecs[i]) == std::string::npos) read = false; for(int i=0;i<file_non_specs.size();++i) if(filename.find(file_non_specs[i]) != std::string::npos) read = false;

if(read){ if(verbose) std::cout << "reading " << filepath << std::endl; PixelMap map(filepath,resolution); images.push_back(std::move(map)); images.push_back(PixelMap(filepath,resolution)); names.push_back(filepath); ++count; } } }

closedir( dp );

std::cout << count << " fits files read." << std::endl; return ; }

Parameters

dir	path to directory containing fits files
filespec	string of charactors in file name that are matched. It can be an empty string.
filenames	output vector of PixelMaps
file_non_spec	string of charactors in file name that file must not have.

shuffle()

template<typename T , typename R >

void Utilities::shuffle	(	std::vector< T > &	vec,
		R &	ran )

Shuffles a vector into a random order.

Parameters

vec	The vector to be shuffled
ran	a random number generator so that ran() gives a number between 0 and 1

sort_indexes() [1/2]

template<typename T >

void Utilities::sort_indexes	(	const std::vector< T > &	v,
		std::vector< size_t > &	index )

Find the indexes that sort a vector in asending order.

Parameters

v	the original data that is not changed
index	vector of indexes that if put into v will sort it

sort_indexes() [2/2]

template<typename T >

void Utilities::sort_indexes	(	const T *	v,
		std::vector< size_t > &	index,
		size_t	N )

Parameters

v	the original data that is not changed
index	vector of indexes that if put into v will sort it

sort_indexes_decending()

template<typename T >

void Utilities::sort_indexes_decending	(	const std::vector< T > &	v,
		std::vector< size_t > &	index )

Find the indexes that sort a vector in descending order.

Parameters

v	the original data that is not changed
index	vector of indexes that if put into v will sort it

thicken_poly()

std::vector< std::vector< Point_2d > > Utilities::thicken_poly	(	const std::vector< Point_2d > &	v,
		double	R )

return the boundaries of the region that is within R of the curve v

If the radius R is small compared to the dimenstions of the polygon an approximation is make to save time in which parallel segments to each side of the polygon are used.

If the radius R is not small a gridding method is used which is more reliable.

TighterHull()

std::vector< Point_2d > Utilities::TighterHull

(

const std::vector< Point_2d > &

v

)

removes the intersections while removing interior loops The input curve needs to be ordered already. No points in the input curve will be outside the output hull. Will fail if there are overlapping segments on the hull.

TwoDInterpolator()

PosType Utilities::TwoDInterpolator	(	PosType *	x,
		int	Npixels,
		PosType	range,
		PosType *	center,
		PosType *	map,
		bool	init )

bilinear interpolation from a map.

Out of bounds points return 0. map is a i dimensional array representing a 2 dimensional map. Don't use init. After it is used once, later calls can use TwoDInterpolator(PosType *map) for the same point in the same coordinate system to save time in calculating the indexes.

bilinear interpolation

bilinear interpolation

windings() [1/5]

int Utilities::windings	(	const Point_2d &	x,
		const std::vector< Point_2d > &	point,
		PosType *	area )

Parameters

x	Point for which the winding number is calculated
point	The points on the border. These must be ordered.
area	returns absolute the area within the curve with oriented border

windings() [2/5]

int Utilities::windings	(	Point_2d &	x,
		std::vector< RAY > &	point,
		PosType *	area )

Parameters

x	Point for which the winding number is calculated
point	The points on the border. Uses image plane points
area	returns absolute the area within the curve with oriented border

windings() [3/5]

int Utilities::windings	(	PosType *	x,
		Kist< Point > *	kist,
		PosType *	area,
		short	image = 0 )

Parameters

x	Point for which the winding number is calculated
kist	Kist of points on the border. These must be ordered.
area	returns absolute the area within the curve with oriented border
image	if == 1 the image of the curve is uses as the curve

windings() [4/5]

int Utilities::windings	(	PosType *	x,
		Point **	points,
		unsigned long	Npoints,
		PosType *	area,
		short	image = 0 )

Parameters

x	Point for which the winding number is calculated
points	The points on the border. These must be ordered.
Npoints	number of points in curve
area	returns absolute the area within the curve with oriented border
image	if == 1 the image of the curve is uses as the curve

windings() [5/5]

int Utilities::windings	(	PosType *	x,
		Point *	points,
		unsigned long	Npoints,
		PosType *	area,
		short	image )

windings(): winding number test for a point in a polygon Returns: Number of times a curves winds around the point x.

The number of times the curve loops around a point is calculated.

The area of a self-intersecting curve will be the area of the regions encircled in a clockwise direction minus the regions encircled in a counterclockwise direction - an infinity symbol has zero area.

Parameters

x	Point for which the winding number is calculated
points	The points on the border. These must be ordered.
Npoints	number of points in curve
area	returns absolute the area within the curve with oriented border
image	if == 1 the image of the curve is uses as the curve

windings2()

int Utilities::windings2	(	PosType *	x,
		Point *	points_original,
		unsigned long	Npoints,
		PosType *	area,
		short	image )

determines whether a point is inside a curve, that has been stretched 1.2 times returns the area of the stretched curve

Parameters

x	Point for which the winding number is calculated
points_original	The points on the border. These must be ordered.
Npoints	number of points in curve
area	returns absolute the area within the curve with oriented border
image	if == 0 the image of the curve is used as the curve

writeCurves()

void Utilities::writeCurves	(	int	m,
		ImageInfo *	critical,
		int	Ncrit,
		int	ind_caustic )

writes in four files the critical curves and the caustics for all the curves found and also for a specified one (ind_causic)

Parameters

m	part of te filename, could be the number/index of the main lens
critical	the crit curve
Ncrit	the number of crit curves
ind_caustic	the index of the cuvre of interest