codebase / stimlib

Commit 72175dcf5c397dbdc95061cfc0451141a8617943

Authored by Jiaming Guo
1 parent cea40664

add stitch function in centerline for pavel

Showing 1 changed file with 357 additions and 0 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/centerline.h
  Show/Hide comments   View file @72175dc
... ... @@ -3,6 +3,7 @@
3 3  
4 4 #include <vector>
5 5 #include <stim/math/vec3.h>
  6 +#include <stim/structures/kdtree.cuh>
6 7  
7 8 namespace stim{
8 9  
... ... @@ -134,6 +135,362 @@ public:
134 135 return L.back();
135 136 }
136 137  
  138 + /// stitch two centerlines
  139 + ///@param c1, c2: two centerlines
  140 + ///@param sigma: sample rate
  141 + static std::vector< stim::centerline<T> > stitch(stim::centerline<T> c1, stim::centerline<T> c2 = stim::centerline<T>()) {
  142 +
  143 + std::vector< stim::centerline<T> > result;
  144 + stim::centerline<T> new_centerline;
  145 + stim::vec3<T> new_vertex;
  146 +
  147 + // if only one centerline, stitch itself!
  148 + if (c2.size() == 0) {
  149 + size_t num = c1.size();
  150 + size_t id = 100000; // store the downsample start position
  151 + T threshold;
  152 + if (num < 4) { // if the number of vertex is less than 4, do nothing
  153 + result.push_back(c1);
  154 + return result;
  155 + }
  156 + else {
  157 + // test geometry start vertex
  158 + stim::vec3<T> v1 = c1[1] - c1[0]; // vector from c1[0] to c1[1]
  159 + for (size_t p = 2; p < num; p++) { // 90° standard???
  160 + stim::vec3<T> v2 = c1[p] - c1[0];
  161 + float cosine = v2.dot(v1);
  162 + if (cosine < 0) {
  163 + id = p;
  164 + threshold = v2.len();
  165 + break;
  166 + }
  167 + }
  168 + if (id != 100000) { // find a downsample position on the centerline
  169 + T* c;
  170 + c = (T*)malloc(sizeof(T) * (num - id) * 3);
  171 + for (size_t p = id; p < num; p++) {
  172 + for (size_t d = 0; d < 3; d++) {
  173 + c[p * 3 + d] = c1[p][d];
  174 + }
  175 + }
  176 + stim::kdtree<T, 3> kdt;
  177 + kdt.create(c, num - id, 5); // create tree
  178 +
  179 + T* query = (T*)malloc(sizeof(T) * 3);
  180 + for (size_t d = 0; d < 3; d++)
  181 + query[d] = c1[0][d];
  182 + size_t index;
  183 + T dist;
  184 +
  185 + kdt.search(query, 1, &index, &dist);
  186 +
  187 + free(query);
  188 + free(c);
  189 +
  190 + if (dist > threshold) {
  191 + result.push_back(c1);
  192 + }
  193 + else {
  194 + // the loop part
  195 + new_vertex = c1[index];
  196 + new_centerline.push_back(new_vertex);
  197 + for (size_t p = 0; p < index + 1; p++) {
  198 + new_vertex = c1[p];
  199 + new_centerline.push_back(new_vertex);
  200 + }
  201 + result.push_back(new_centerline);
  202 + new_centerline.clear();
  203 +
  204 + // the tail part
  205 + for (size_t p = index; p < num; p++) {
  206 + new_vertex = c1[p];
  207 + new_centerline.push_back(new_vertex);
  208 + }
  209 + result.push_back(new_centerline);
  210 + }
  211 + }
  212 + else { // there is one potential problem that two positions have to be stitched
  213 + // test geometry end vertex
  214 + stim::vec3<T> v1 = c1[num - 2] - c1[num - 1];
  215 + for (size_t p = num - 2; p > 0; p--) { // 90° standard
  216 + stim::vec3<T> v2 = c1[p - 1] - c1[num - 1];
  217 + float cosine = v2.dot(v1);
  218 + if (cosine < 0) {
  219 + id = p;
  220 + threshold = v2.len();
  221 + break;
  222 + }
  223 + }
  224 + if (id != 100000) { // find a downsample position
  225 + T* c;
  226 + c = (T*)malloc(sizeof(T) * (id + 1) * 3);
  227 + for (size_t p = 0; p < id + 1; p++) {
  228 + for (size_t d = 0; d < 3; d++) {
  229 + c[p * 3 + d] = c1[p][d];
  230 + }
  231 + }
  232 + stim::kdtree<T, 3> kdt;
  233 + kdt.create(c, id + 1, 5); // create tree
  234 +
  235 + T* query = (T*)malloc(sizeof(T) * 1 * 3);
  236 + for (size_t d = 0; d < 3; d++)
  237 + query[d] = c1[num - 1][d];
  238 + size_t index;
  239 + T dist;
  240 +
  241 + kdt.search(query, 1, &index, &dist);
  242 +
  243 + free(query);
  244 + free(c);
  245 +
  246 + if (dist > threshold) {
  247 + result.push_back(c1);
  248 + }
  249 + else {
  250 + // the tail part
  251 + for (size_t p = 0; p < index + 1; p++) {
  252 + new_vertex = c1[p];
  253 + new_centerline.push_back(new_vertex);
  254 + }
  255 + result.push_back(new_centerline);
  256 + new_centerline.clear();
  257 +
  258 + // the loop part
  259 + for (size_t p = index; p < num; p++) {
  260 + new_vertex = c1[p];
  261 + new_centerline.push_back(new_vertex);
  262 + }
  263 + new_vertex = c1[index];
  264 + new_centerline.push_back(new_vertex);
  265 + result.push_back(new_centerline);
  266 + }
  267 + }
  268 + else { // no stitch position
  269 + result.push_back(c1);
  270 + }
  271 + }
  272 + }
  273 + }
  274 +
  275 +
  276 + // two centerlines
  277 + else {
  278 + // find stitch position based on nearest neighbors
  279 + size_t num1 = c1.size();
  280 + T* c = (T*)malloc(sizeof(T) * num1 * 3); // c1 as reference point
  281 + for (size_t p = 0; p < num1; p++) // centerline to array
  282 + for (size_t d = 0; d < 3; d++) // because right now my kdtree code is a relatively close code, it has its own structure
  283 + c[p * 3 + d] = c1[p][d]; // I will merge it into stimlib totally in the near future
  284 +
  285 + stim::kdtree<T, 3> kdt; // kdtree object
  286 + kdt.create(c, num1, 5); // create tree
  287 +
  288 + size_t num2 = c2.size();
  289 + T* query = (T*)malloc(sizeof(T) * num2 * 3); // c2 as query point
  290 + for (size_t p = 0; p < num2; p++) {
  291 + for (size_t d = 0; d < 3; d++) {
  292 + query[p * 3 + d] = c2[p][d];
  293 + }
  294 + }
  295 + std::vector<size_t> index(num2);
  296 + std::vector<T> dist(num2);
  297 +
  298 + kdt.search(query, num2, &index[0], &dist[0]); // find the nearest neighbors in c1 for c2
  299 +
  300 + // clear up
  301 + free(query);
  302 + free(c);
  303 +
  304 + // find the average vertex distance of one centerline
  305 + T sigma1 = 0;
  306 + T sigma2 = 0;
  307 + for (size_t p = 0; p < num1 - 1; p++)
  308 + sigma1 += (c1[p] - c1[p + 1]).len();
  309 + for (size_t p = 0; p < num2 - 1; p++)
  310 + sigma2 += (c2[p] - c2[p + 1]).len();
  311 + sigma1 /= (num1 - 1);
  312 + sigma2 /= (num2 - 1);
  313 + float threshold = 4 * (sigma1 + sigma2) / 2; // better way to do this?
  314 +
  315 + T min_d = *std::min_element(dist.begin(), dist.end()); // find the minimum distance between c1 and c2
  316 +
  317 + if (min_d > threshold) { // if the minimum distance is too large
  318 + result.push_back(c1);
  319 + result.push_back(c2);
  320 +
  321 +#ifdef DEBUG
  322 + std::cout << "The distance between these two centerlines is too large" << std::endl;
  323 +#endif
  324 + }
  325 + else {
  326 + auto smallest = std::min_element(dist.begin(), dist.end());
  327 + auto i = std::distance(dist.begin(), smallest); // find the index of min-distance in distance list
  328 +
  329 + // stitch position in c1 and c2
  330 + int id1 = index[i];
  331 + int id2 = i;
  332 +
  333 + // actually there are two cases
  334 + // first one inacceptable
  335 + // second one acceptable
  336 + if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) { // only stitch one end vertex to another centerline
  337 + result.push_back(c1);
  338 + result.push_back(c2);
  339 + }
  340 + else {
  341 + if (id1 == 0 || id1 == num1 - 1) { // if the stitch vertex is the first or last vertex of c1
  342 + // for c2, consider two cases(one degenerate case)
  343 + if (id2 == 0 || id2 == num2 - 1) { // case 1, if stitch position is also on the end of c2
  344 + // we have to decide which centerline get a new vertex, based on direction
  345 + // for c1, computer the direction change angle
  346 + stim::vec3<T> v1, v2;
  347 + float alpha1, alpha2; // direction change angle
  348 + if (id1 == 0)
  349 + v1 = (c1[1] - c1[0]).norm();
  350 + else
  351 + v1 = (c1[num1 - 2] - c1[num1 - 1]).norm();
  352 + v2 = (c2[id2] - c1[id1]).norm();
  353 + alpha1 = v1.dot(v2);
  354 + if (id2 == 0)
  355 + v1 = (c2[1] - c2[0]).norm();
  356 + else
  357 + v1 = (c2[num2 - 2] - c2[num2 - 1]).norm();
  358 + v2 = (c1[id1] - c2[id2]).norm();
  359 + alpha2 = v1.dot(v2);
  360 + if (abs(alpha1) > abs(alpha2)) { // add the vertex to c1 in order to get smooth connection
  361 + // push back c1
  362 + if (id1 == 0) { // keep geometry information
  363 + new_vertex = c2[id2];
  364 + new_centerline.push_back(new_vertex);
  365 + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
  366 + new_vertex = c1[p];
  367 + new_centerline.push_back(new_vertex);
  368 + }
  369 + }
  370 + else {
  371 + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
  372 + new_vertex = c1[p];
  373 + new_centerline.push_back(new_vertex);
  374 + }
  375 + new_vertex = c2[id2];
  376 + new_centerline.push_back(new_vertex);
  377 + }
  378 + result.push_back(new_centerline);
  379 + new_centerline.clear();
  380 +
  381 + // push back c2
  382 + for (size_t p = 0; p < num2; p++) {
  383 + new_vertex = c2[p];
  384 + new_centerline.push_back(new_vertex);
  385 + }
  386 + result.push_back(new_centerline);
  387 + }
  388 + else { // add the vertex to c2 in order to get smooth connection
  389 + // push back c1
  390 + for (size_t p = 0; p < num1; p++) {
  391 + new_vertex = c1[p];
  392 + new_centerline.push_back(new_vertex);
  393 + }
  394 + result.push_back(new_centerline);
  395 + new_centerline.clear();
  396 +
  397 + // push back c2
  398 + if (id2 == 0) { // keep geometry information
  399 + new_vertex = c1[id1];
  400 + new_centerline.push_back(new_vertex);
  401 + for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
  402 + new_vertex = c2[p];
  403 + new_centerline.push_back(new_vertex);
  404 + }
  405 + }
  406 + else {
  407 + for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
  408 + new_vertex = c2[p];
  409 + new_centerline.push_back(new_vertex);
  410 + }
  411 + new_vertex = c1[id1];
  412 + new_centerline.push_back(new_vertex);
  413 + }
  414 + result.push_back(new_centerline);
  415 + }
  416 + }
  417 + else { // case 2, the stitch position is on c2
  418 + // push back c1
  419 + if (id1 == 0) { // keep geometry information
  420 + new_vertex = c2[id2];
  421 + new_centerline.push_back(new_vertex);
  422 + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
  423 + new_vertex = c1[p];
  424 + new_centerline.push_back(new_vertex);
  425 + }
  426 + }
  427 + else {
  428 + for (size_t p = 0; p < num1; p++) { // geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2
  429 + new_vertex = c1[p];
  430 + new_centerline.push_back(new_vertex);
  431 + }
  432 + new_vertex = c2[id2];
  433 + new_centerline.push_back(new_vertex);
  434 + }
  435 + result.push_back(new_centerline);
  436 + new_centerline.clear();
  437 +
  438 + // push back c2
  439 + for (size_t p = 0; p < id2 + 1; p++) { // first part
  440 + new_vertex = c2[p];
  441 + new_centerline.push_back(new_vertex);
  442 + }
  443 + result.push_back(new_centerline);
  444 + new_centerline.clear();
  445 +
  446 + for (size_t p = id2; p < num2; p++) { // second part
  447 + new_vertex = c2[p];
  448 + new_centerline.push_back(new_vertex);
  449 + }
  450 + result.push_back(new_centerline);
  451 + }
  452 + }
  453 + else { // if the stitch vertex is the first or last vertex of c2
  454 + // push back c2
  455 + if (id2 == 0) { // keep geometry information
  456 + new_vertex = c1[id1];
  457 + new_centerline.push_back(new_vertex);
  458 + for (size_t p = 0; p < num2; p++) { // stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2
  459 + new_vertex = c2[p];
  460 + new_centerline.push_back(new_vertex);
  461 + }
  462 + }
  463 + else {
  464 + for (size_t p = 0; p < num2; p++) { // geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1
  465 + new_vertex = c2[p];
  466 + new_centerline.push_back(new_vertex);
  467 + }
  468 + new_vertex = c1[id1];
  469 + new_centerline.push_back(new_vertex);
  470 + result.push_back(new_centerline);
  471 + new_centerline.clear();
  472 +
  473 + // push back c1
  474 + for (size_t p = 0; p < id1 + 1; p++) { // first part
  475 + new_vertex = c1[p];
  476 + new_centerline.push_back(new_vertex);
  477 + }
  478 + result.push_back(new_centerline);
  479 + new_centerline.clear();
  480 +
  481 + for (size_t p = id1; p < num1; p++) { // second part
  482 + new_vertex = c1[p];
  483 + new_centerline.push_back(new_vertex);
  484 + }
  485 + result.push_back(new_centerline);
  486 + }
  487 + }
  488 + }
  489 + }
  490 + }
  491 + return result;
  492 + }
  493 +
137 494 /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
138 495 std::vector< stim::centerline<T> > split(unsigned int idx){
139 496  
... ...