From 72175dcf5c397dbdc95061cfc0451141a8617943 Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Thu, 9 Mar 2017 16:43:46 -0600 Subject: [PATCH] add stitch function in centerline for pavel --- stim/biomodels/centerline.h | 357 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 357 insertions(+), 0 deletions(-) diff --git a/stim/biomodels/centerline.h b/stim/biomodels/centerline.h index 46d9369..a464a71 100644 --- a/stim/biomodels/centerline.h +++ b/stim/biomodels/centerline.h @@ -3,6 +3,7 @@ #include #include +#include namespace stim{ @@ -134,6 +135,362 @@ public: return L.back(); } + /// stitch two centerlines + ///@param c1, c2: two centerlines + ///@param sigma: sample rate + static std::vector< stim::centerline > stitch(stim::centerline c1, stim::centerline c2 = stim::centerline()) { + + std::vector< stim::centerline > result; + stim::centerline new_centerline; + stim::vec3 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 v1 = c1[1] - c1[0]; // vector from c1[0] to c1[1] + for (size_t p = 2; p < num; p++) { // 90° standard??? + stim::vec3 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 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 v1 = c1[num - 2] - c1[num - 1]; + for (size_t p = num - 2; p > 0; p--) { // 90° standard + stim::vec3 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 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 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 index(num2); + std::vector 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()); + auto 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 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 > split(unsigned int idx){ -- libgit2 0.21.4