centerline.h 19.4 KB
``````#ifndef STIM_CENTERLINE_H
#define STIM_CENTERLINE_H

#include <vector>
#include <stim/math/vec3.h>
#include <stim/structures/kdtree.cuh>

namespace stim{

/**	This class stores information about a single fiber represented as a set of geometric points
*	between two branch or end points. This class is used as a fundamental component of the stim::network
*	class to describe an interconnected (often biological) network.
*/
template<typename T>
class centerline : public std::vector< stim::vec3<T> >{

protected:

std::vector<T> L;										//stores the integrated length along the fiber (used for parameterization)

///Return the normalized direction vector at point i (average of the incoming and outgoing directions)
vec3<T> d(size_t i) {
if (size() <= 1) return vec3<T>(0, 0, 0);						//if there is insufficient information to calculate the direction, return a null vector
if (size() == 2) return (at(1) - at(0)).norm();					//if there are only two points, the direction vector at both is the direction of the line segment
if (i == 0) return (at(1) - at(0)).norm();						//the first direction vector is oriented towards the first line segment
if (i == size() - 1) return (at(size() - 1) - at(size() - 2)).norm();	//the last direction vector is oriented towards the last line segment

//all other direction vectors are the average direction of the two joined line segments
vec3<T> a = at(i) - at(i - 1);
vec3<T> b = at(i + 1) - at(i);
vec3<T> ab = a.norm() + b.norm();
return ab.norm();
}

//initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
void update_L(size_t start = 0) {
L.resize(size());									//allocate space for the L array
if (start == 0) {
L[0] = 0;											//initialize the length value for the first point to zero (0)
start++;
}

stim::vec3<T> d;
for (size_t i = start; i < size(); i++) {		//for each line segment in the centerline
d = at(i) - at(i - 1);
L[i] = L[i - 1] + d.len();				//calculate the running length total
}
}

void init() {
if (size() == 0) return;								//return if there aren't any points
update_L();
}

/// Returns a stim::vec representing the point at index i

/// @param i is an index of the desired centerline point
stim::vec<T> get_vec(unsigned i){
return std::vector< stim::vec3<T> >::at(i);
}

///finds the index of the point closest to the length l on the lower bound.
///binary search.
size_t findIdx(T l) {
for (size_t i = 1; i < L.size(); i++) {				//for each point in the centerline
if (L[i] > l) return i - 1;						//if we have passed the desired length value, return i-1
}
return L.size() - 1;
/*size_t i = L.size() / 2;
size_t max = L.size() - 1;
size_t min = 0;
while (i < L.size() - 1){
if (l < L[i]) {
max = i;
i = min + (max - min) / 2;
}
else if (L[i] <= l && L[i + 1] >= l) {
break;
}
else {
min = i;
i = min + (max - min) / 2;
}
}
return i;*/
}

///Returns a position vector at the given length into the fiber (based on the pvalue).
///Interpolates the radius along the line.
///@param l: the location of the in the cylinder.
///@param idx: integer location of the point closest to l but prior to it.
stim::vec3<T> p(T l, int idx) {
T rat = (l - L[idx]) / (L[idx + 1] - L[idx]);
stim::vec3<T> v1 = at(idx);
stim::vec3<T> v2 = at(idx + 1);
return(v1 + (v2 - v1)*rat);
}

public:

using std::vector< stim::vec3<T> >::at;
using std::vector< stim::vec3<T> >::size;

centerline() : std::vector< stim::vec3<T> >() {
init();
}
centerline(size_t n) : std::vector< stim::vec3<T> >(n){
init();
}
centerline(std::vector<stim::vec3<T> > pos) :
std::vector<stim::vec3<T> > (pos)
{
init();
}

//overload the push_back function to update the length vector
void push_back(stim::vec3<T> p) {
std::vector< stim::vec3<T> >::push_back(p);
update_L(size() - 1);
}

///Returns a position vector at the given p-value (p value ranges from 0 to 1).
///interpolates the position along the line.
///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
stim::vec3<T> p(T pvalue) {
if (pvalue <= 0.0) return at(0);			//return the first element
if (pvalue >= 1.0) return back();			//return the last element

T l = pvalue*L[L.size() - 1];
int idx = findIdx(l);
return p(l, idx);
}

///Update centerline internal parameters (currently the L vector)
void update() {
init();
}
///Return the length of the entire centerline
T length() {
return L.back();
}

/// stitch two centerlines
///@param c1, c2: two centerlines
///@param sigma: sample rate
static std::vector< stim::centerline<T> > stitch(stim::centerline<T> c1, stim::centerline<T> c2 = stim::centerline<T>()) {

std::vector< stim::centerline<T> > result;
stim::centerline<T> new_centerline;
stim::vec3<T> new_vertex;

// if only one centerline, stitch itself!
if (c2.size() == 0) {
size_t num = c1.size();
size_t id = 100000;							// store the downsample start position
T threshold;
if (num < 4) {								// if the number of vertex is less than 4, do nothing
result.push_back(c1);
return result;
}
else {
// test geometry start vertex
stim::vec3<T> v1 = c1[1] - c1[0];		// vector from c1[0] to c1[1]
for (size_t p = 2; p < num; p++) {		// 90ยฐ standard???
stim::vec3<T> v2 = c1[p] - c1[0];
float cosine = v2.dot(v1);
if (cosine < 0) {
id = p;
threshold = v2.len();
break;
}
}
if (id != 100000) {						// find a downsample position on the centerline
T* c;
c = (T*)malloc(sizeof(T) * (num - id) * 3);
for (size_t p = id; p < num; p++) {
for (size_t d = 0; d < 3; d++) {
c[p * 3 + d] = c1[p][d];
}
}
stim::kdtree<T, 3> kdt;
kdt.create(c, num - id, 5);			// create tree

T* query = (T*)malloc(sizeof(T) * 3);
for (size_t d = 0; d < 3; d++)
query[d] = c1[0][d];
size_t index;
T dist;

kdt.search(query, 1, &index, &dist);

free(query);
free(c);

if (dist > threshold) {
result.push_back(c1);
}
else {
// the loop part
new_vertex = c1[index];
new_centerline.push_back(new_vertex);
for (size_t p = 0; p < index + 1; p++) {
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
new_centerline.clear();

// the tail part
for (size_t p = index; p < num; p++) {
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
}
}
else {	// there is one potential problem that two positions have to be stitched
// test geometry end vertex
stim::vec3<T> v1 = c1[num - 2] - c1[num - 1];
for (size_t p = num - 2; p > 0; p--) {		// 90ยฐ standard
stim::vec3<T> v2 = c1[p - 1] - c1[num - 1];
float cosine = v2.dot(v1);
if (cosine < 0) {
id = p;
threshold = v2.len();
break;
}
}
if (id != 100000) {						// find a downsample position
T* c;
c = (T*)malloc(sizeof(T) * (id + 1) * 3);
for (size_t p = 0; p < id + 1; p++) {
for (size_t d = 0; d < 3; d++) {
c[p * 3 + d] = c1[p][d];
}
}
stim::kdtree<T, 3> kdt;
kdt.create(c, id + 1, 5);				// create tree

T* query = (T*)malloc(sizeof(T) * 1 * 3);
for (size_t d = 0; d < 3; d++)
query[d] = c1[num - 1][d];
size_t index;
T dist;

kdt.search(query, 1, &index, &dist);

free(query);
free(c);

if (dist > threshold) {
result.push_back(c1);
}
else {
// the tail part
for (size_t p = 0; p < index + 1; p++) {
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
new_centerline.clear();

// the loop part
for (size_t p = index; p < num; p++) {
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
new_vertex = c1[index];
new_centerline.push_back(new_vertex);
result.push_back(new_centerline);
}
}
else {	// no stitch position
result.push_back(c1);
}
}
}
}

// two centerlines
else {
// find stitch position based on nearest neighbors
size_t num1 = c1.size();
T* c = (T*)malloc(sizeof(T) * num1 * 3);		// c1 as reference point
for (size_t p = 0; p < num1; p++)				// centerline to array
for (size_t d = 0; d < 3; d++)				// because right now my kdtree code is a relatively close code, it has its own structure
c[p * 3 + d] = c1[p][d];				// I will merge it into stimlib totally in the near future

stim::kdtree<T, 3> kdt;							// kdtree object
kdt.create(c, num1, 5);							// create tree

size_t num2 = c2.size();
T* query = (T*)malloc(sizeof(T) * num2 * 3);	// c2 as query point
for (size_t p = 0; p < num2; p++) {
for (size_t d = 0; d < 3; d++) {
query[p * 3 + d] = c2[p][d];
}
}
std::vector<size_t> index(num2);
std::vector<T> dist(num2);

kdt.search(query, num2, &index[0], &dist[0]);	// find the nearest neighbors in c1 for c2

// clear up
free(query);
free(c);

// find the average vertex distance of one centerline
T sigma1 = 0;
T sigma2 = 0;
for (size_t p = 0; p < num1 - 1; p++)
sigma1 += (c1[p] - c1[p + 1]).len();
for (size_t p = 0; p < num2 - 1; p++)
sigma2 += (c2[p] - c2[p + 1]).len();
sigma1 /= (num1 - 1);
sigma2 /= (num2 - 1);
float threshold = 4 * (sigma1 + sigma2) / 2;			// better way to do this?

T min_d = *std::min_element(dist.begin(), dist.end());	// find the minimum distance between c1 and c2

if (min_d > threshold) {								// if the minimum distance is too large
result.push_back(c1);
result.push_back(c2);

#ifdef DEBUG
std::cout << "The distance between these two centerlines is too large" << std::endl;
#endif
}
else {
//	auto smallest = std::min_element(dist.begin(), dist.end());
unsigned int smallest = std::min_element(dist.begin(), dist.end());
//	auto i = std::distance(dist.begin(), smallest);		// find the index of min-distance in distance list
unsigned int i = std::distance(dist.begin(), smallest);		// find the index of min-distance in distance list

// stitch position in c1 and c2
int id1 = index[i];
int id2 = i;

// actually there are two cases
// first one inacceptable
// second one acceptable
if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) {		// only stitch one end vertex to another centerline
result.push_back(c1);
result.push_back(c2);
}
else {
if (id1 == 0 || id1 == num1 - 1) {			// if the stitch vertex is the first or last vertex of c1
// for c2, consider two cases(one degenerate case)
if (id2 == 0 || id2 == num2 - 1) {		// case 1, if stitch position is also on the end of c2
// we have to decide which centerline get a new vertex, based on direction
// for c1, computer the direction change angle
stim::vec3<T> v1, v2;
float alpha1, alpha2;				// direction change angle
if (id1 == 0)
v1 = (c1[1] - c1[0]).norm();
else
v1 = (c1[num1 - 2] - c1[num1 - 1]).norm();
v2 = (c2[id2] - c1[id1]).norm();
alpha1 = v1.dot(v2);
if (id2 == 0)
v1 = (c2[1] - c2[0]).norm();
else
v1 = (c2[num2 - 2] - c2[num2 - 1]).norm();
v2 = (c1[id1] - c2[id2]).norm();
alpha2 = v1.dot(v2);
if (abs(alpha1) > abs(alpha2)) {					// add the vertex to c1 in order to get smooth connection
// push back c1
if (id1 == 0) {									// keep geometry information
new_vertex = c2[id2];
new_centerline.push_back(new_vertex);
for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
}
else {
for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
new_vertex = c2[id2];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
new_centerline.clear();

// push back c2
for (size_t p = 0; p < num2; p++) {
new_vertex = c2[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
}
else {												// add the vertex to c2 in order to get smooth connection
// push back c1
for (size_t p = 0; p < num1; p++) {
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
new_centerline.clear();

// push back c2
if (id2 == 0) {									// keep geometry information
new_vertex = c1[id1];
new_centerline.push_back(new_vertex);
for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
new_vertex = c2[p];
new_centerline.push_back(new_vertex);
}
}
else {
for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
new_vertex = c2[p];
new_centerline.push_back(new_vertex);
}
new_vertex = c1[id1];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
}
}
else {												// case 2, the stitch position is on c2
// push back c1
if (id1 == 0) {									// keep geometry information
new_vertex = c2[id2];
new_centerline.push_back(new_vertex);
for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
}
else {
for (size_t p = 0; p < num1; p++) {			// geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
new_vertex = c2[id2];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
new_centerline.clear();

// push back c2
for (size_t p = 0; p < id2 + 1; p++) {			// first part
new_vertex = c2[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
new_centerline.clear();

for (size_t p = id2; p < num2; p++) {			// second part
new_vertex = c2[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
}
}
else {							// if the stitch vertex is the first or last vertex of c2
// push back c2
if (id2 == 0) {										// keep geometry information
new_vertex = c1[id1];
new_centerline.push_back(new_vertex);
for (size_t p = 0; p < num2; p++) {				// stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2
new_vertex = c2[p];
new_centerline.push_back(new_vertex);
}
}
else {
for (size_t p = 0; p < num2; p++) {				// geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1
new_vertex = c2[p];
new_centerline.push_back(new_vertex);
}
new_vertex = c1[id1];
new_centerline.push_back(new_vertex);
result.push_back(new_centerline);
new_centerline.clear();

// push back c1
for (size_t p = 0; p < id1 + 1; p++) {			// first part
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
new_centerline.clear();

for (size_t p = id1; p < num1; p++) {			// second part
new_vertex = c1[p];
new_centerline.push_back(new_vertex);
}
result.push_back(new_centerline);
}
}
}
}
}
return result;
}

/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
std::vector< stim::centerline<T> > split(unsigned int idx){

std::vector< stim::centerline<T> > fl;				//create an array to store up to two fibers
size_t N = size();

//if the index is an end point, only the existing fiber is returned
if(idx == 0 || idx == N-1){
fl.resize(1);							//set the size of the fiber to 1
fl[0] = *this;							//copy the current fiber
}

//if the index is not an end point
else{

unsigned int N1 = idx + 1;					//calculate the size of both fibers
unsigned int N2 = N - idx;

fl.resize(2);								//set the array size to 2

fl[0] = stim::centerline<T>(N1);			//set the size of each fiber
fl[1] = stim::centerline<T>(N2);

//copy both halves of the fiber
unsigned int i;

//first half
for(i = 0; i < N1; i++)					//for each centerline point
fl[0][i] = std::vector< stim::vec3<T> >::at(i);
fl[0].init();							//initialize the length vector

//second half
for(i = 0; i < N2; i++)
fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
fl[1].init();							//initialize the length vector
}

return fl;										//return the array

}

/// Outputs the fiber as a string
std::string str(){
std::stringstream ss;
size_t N = std::vector< stim::vec3<T> >::size();
ss << "---------[" << N << "]---------" << std::endl;
for (size_t i = 0; i < N; i++)
ss << std::vector< stim::vec3<T> >::at(i) << std::endl;
ss << "--------------------" << std::endl;

return ss.str();
}

/// Back method returns the last point in the fiber
stim::vec3<T> back(){
return std::vector< stim::vec3<T> >::back();
}

////resample a fiber in the network
stim::centerline<T> resample(T spacing)
{
//std::cout<<"fiber::resample()"<<std::endl;

stim::vec3<T> v;    //v-direction vector of the segment
stim::vec3<T> p;      //- intermediate point to be added
stim::vec3<T> p1;   // p1 - starting point of an segment on the fiber,
stim::vec3<T> p2;   // p2 - ending point,
//double sum=0;  //distance summation

size_t N = size();

centerline<T> new_c; // initialize list of new resampled points on the fiber
// for each point on the centerline (skip if it is the last point on centerline)
for(unsigned int f=0; f< N-1; f++)
{
p1 = at(f);
p2 = at(f+1);
v = p2 - p1;

T lengthSegment = v.len();			//find Length of the segment as distance between the starting and ending points of the segment

if(lengthSegment >= spacing){ // if length of the segment is greater than standard deviation resample

// repeat resampling until accumulated stepsize is equsl to length of the segment
for(T step=0.0; step<lengthSegment; step+=spacing){
// calculate the resampled point by travelling step size in the direction of normalized gradient vector
p = p1 + v * (step / lengthSegment);

// add this resampled points to the new fiber list
new_c.push_back(p);
}
}
else       // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
new_c.push_back(at(f));
}
new_c.push_back(at(N-1));   //add the last point on the fiber to the new fiber list
//centerline newFiber(newPointList);
return new_c;
}

};

}	//end namespace stim

#endif
``````