LEFT | RIGHT |
| 1 /* |
| 2 * ***** BEGIN GPL LICENSE BLOCK ***** |
| 3 * |
| 4 * This program is free software; you can redistribute it and/or |
| 5 * modify it under the terms of the GNU General Public License |
| 6 * as published by the Free Software Foundation; either version 2 |
| 7 * of the License, or (at your option) any later version. |
| 8 * |
| 9 * This program is distributed in the hope that it will be useful, |
| 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 12 * GNU General Public License for more details. |
| 13 * |
| 14 * You should have received a copy of the GNU General Public License |
| 15 * along with this program; if not, write to the Free Software Foundation, |
| 16 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
| 17 * |
| 18 * Contributor(s): Tao Ju |
| 19 * |
| 20 * ***** END GPL LICENSE BLOCK ***** |
| 21 */ |
| 22 |
1 #include "octree.h" | 23 #include "octree.h" |
2 | 24 #include <Eigen/Dense> |
| 25 #include <limits> |
3 #include <time.h> | 26 #include <time.h> |
4 | 27 |
5 /** | 28 /** |
6 * Implementations of Octree member functions. | 29 * Implementations of Octree member functions. |
7 * | 30 * |
8 * @author Tao Ju | 31 * @author Tao Ju |
9 */ | 32 */ |
10 | 33 |
11 /* set to non-zero value to enable debugging output */ | 34 /* set to non-zero value to enable debugging output */ |
12 #define DC_DEBUG 0 | 35 #define DC_DEBUG 0 |
(...skipping 1389 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1402 newnode = patchSplit( newnode, st, len, ylists[3], 2, zlists[6], zlists[
7] ) ; | 1425 newnode = patchSplit( newnode, st, len, ylists[3], 2, zlists[6], zlists[
7] ) ; |
1403 ········ | 1426 ········ |
1404 // Recur | 1427 // Recur |
1405 len >>= 1 ; | 1428 len >>= 1 ; |
1406 int count = 0 ; | 1429 int count = 0 ; |
1407 for ( int i = 0 ; i < 8 ; i ++ ) | 1430 for ( int i = 0 ; i < 8 ; i ++ ) |
1408 { | 1431 { |
1409 if ( zlists[i] != NULL ) | 1432 if ( zlists[i] != NULL ) |
1410 { | 1433 { |
1411 int nori[3] = {· | 1434 int nori[3] = {· |
1412 » » » » st[0] + len * vertMap[i][0] ,· | 1435 » » » » st[0] + len * vertmap[i][0] ,· |
1413 » » » » st[1] + len * vertMap[i][1] ,· | 1436 » » » » st[1] + len * vertmap[i][1] ,· |
1414 » » » » st[2] + len * vertMap[i][2] } ; | 1437 » » » » st[2] + len * vertmap[i][2] } ; |
1415 patch( getChild( newnode , count ), nori, len, zlists[i]
) ; | 1438 patch( getChild( newnode , count ), nori, len, zlists[i]
) ; |
1416 } | 1439 } |
1417 | 1440 |
1418 if ( hasChild( newnode, i ) ) | 1441 if ( hasChild( newnode, i ) ) |
1419 { | 1442 { |
1420 count ++ ; | 1443 count ++ ; |
1421 } | 1444 } |
1422 } | 1445 } |
1423 #ifdef IN_DEBUG_MODE | 1446 #ifdef IN_DEBUG_MODE |
1424 dc_printf("Return from PATCH\n") ; | 1447 dc_printf("Return from PATCH\n") ; |
(...skipping 561 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
1986 float off[3] ; | 2009 float off[3] ; |
1987 int num = 0, num2 = 0 ; | 2010 int num = 0, num2 = 0 ; |
1988 | 2011 |
1989 UCHAR* leafnode = locateLeaf( leaf->pos ) ; | 2012 UCHAR* leafnode = locateLeaf( leaf->pos ) ; |
1990 for ( int i = 0 ; i < 4 ; i ++ ) | 2013 for ( int i = 0 ; i < 4 ; i ++ ) |
1991 { | 2014 { |
1992 int edgeind = faceMap[ dir * 2 ][ i ] ; | 2015 int edgeind = faceMap[ dir * 2 ][ i ] ; |
1993 int nst[3] ; | 2016 int nst[3] ; |
1994 for ( int j = 0 ; j < 3 ; j ++ ) | 2017 for ( int j = 0 ; j < 3 ; j ++ ) |
1995 { | 2018 { |
1996 » » » nst[j] = leaf->pos[j] + mindimen * vertMap[ edgevmap[ ed
geind][ 0 ] ][ j ] ; | 2019 » » » nst[j] = leaf->pos[j] + mindimen * vertmap[ edgemap[ edg
eind][ 0 ] ][ j ] ; |
1997 } | 2020 } |
1998 | 2021 |
1999 if ( getEdgeIntersectionByIndex( nst, edgeind / 4, off, 1 ) ) | 2022 if ( getEdgeIntersectionByIndex( nst, edgeind / 4, off, 1 ) ) |
2000 { | 2023 { |
2001 avg[0] += off[0] ; | 2024 avg[0] += off[0] ; |
2002 avg[1] += off[1] ; | 2025 avg[1] += off[1] ; |
2003 avg[2] += off[2] ; | 2026 avg[2] += off[2] ; |
2004 num ++ ; | 2027 num ++ ; |
2005 } | 2028 } |
2006 if ( getEdgeParity( leafnode, edgeind ) ) | 2029 if ( getEdgeParity( leafnode, edgeind ) ) |
(...skipping 196 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2203 for ( int i = 0 ; i < size ; i ++ ) | 2226 for ( int i = 0 ; i < size ; i ++ ) |
2204 { | 2227 { |
2205 table[i] = 0 ; | 2228 table[i] = 0 ; |
2206 } | 2229 } |
2207 for ( int i = 0 ; i < 256 ; i ++ ) | 2230 for ( int i = 0 ; i < 256 ; i ++ ) |
2208 { | 2231 { |
2209 int ind = 0 ; | 2232 int ind = 0 ; |
2210 for ( int j = 11 ; j >= 0 ; j -- ) | 2233 for ( int j = 11 ; j >= 0 ; j -- ) |
2211 { | 2234 { |
2212 ind <<= 1 ; | 2235 ind <<= 1 ; |
2213 » » » if ( ( ( i >> edgevmap[j][0] ) & 1 ) ^ ( ( i >> edgevmap
[j][1] ) & 1 ) ) | 2236 » » » if ( ( ( i >> edgemap[j][0] ) & 1 ) ^ ( ( i >> edgemap[j
][1] ) & 1 ) ) |
2214 { | 2237 { |
2215 ind |= 1 ; | 2238 ind |= 1 ; |
2216 } | 2239 } |
2217 } | 2240 } |
2218 | 2241 |
2219 table[ ind ] = i ; | 2242 table[ ind ] = i ; |
2220 } | 2243 } |
2221 | 2244 |
2222 // Next, traverse the grid | 2245 // Next, traverse the grid |
2223 int sg = 1 ; | 2246 int sg = 1 ; |
(...skipping 111 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2335 | 2358 |
2336 if ( par == 1 && inp == 0 ) | 2359 if ( par == 1 && inp == 0 ) |
2337 { | 2360 { |
2338 // Intersection edge, hasn't been processed | 2361 // Intersection edge, hasn't been processed |
2339 // Let's start filling | 2362 // Let's start filling |
2340 GridQueue* queue = new GridQueue() ; | 2363 GridQueue* queue = new GridQueue() ; |
2341 int total = 1 ; | 2364 int total = 1 ; |
2342 | 2365 |
2343 // Set to in process | 2366 // Set to in process |
2344 int mst[3] ; | 2367 int mst[3] ; |
2345 » » » » mst[0] = st[0] + vertmap[edgevmap[i][0]][0] * le
n ; | 2368 » » » » mst[0] = st[0] + vertmap[edgemap[i][0]][0] * len
; |
2346 » » » » mst[1] = st[1] + vertmap[edgevmap[i][0]][1] * le
n ; | 2369 » » » » mst[1] = st[1] + vertmap[edgemap[i][0]][1] * len
; |
2347 » » » » mst[2] = st[2] + vertmap[edgevmap[i][0]][2] * le
n; | 2370 » » » » mst[2] = st[2] + vertmap[edgemap[i][0]][2] * len
; |
2348 int mdir = i / 4 ; | 2371 int mdir = i / 4 ; |
2349 setInProcessAll( mst, mdir ) ; | 2372 setInProcessAll( mst, mdir ) ; |
2350 | 2373 |
2351 // Put this edge into queue | 2374 // Put this edge into queue |
2352 queue->pushQueue( mst, mdir ) ; | 2375 queue->pushQueue( mst, mdir ) ; |
2353 | 2376 |
2354 // Queue processing | 2377 // Queue processing |
2355 int nst[3], dir ; | 2378 int nst[3], dir ; |
2356 while ( queue->popQueue( nst, dir ) == 1 ) | 2379 while ( queue->popQueue( nst, dir ) == 1 ) |
2357 { | 2380 { |
(...skipping 59 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2417 } | 2440 } |
2418 } | 2441 } |
2419 ················································ | 2442 ················································ |
2420 if ( eind == 3 || eind == -1 ) | 2443 if ( eind == 3 || eind == -1 ) |
2421 { | 2444 { |
2422 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); | 2445 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); |
2423 } | 2446 } |
2424 else· | 2447 else· |
2425 { | 2448 { |
2426 int est[3] ; | 2449 int est[3] ; |
2427 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgevmap[edge][0]][0] * len ; | 2450 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgemap[edge][0]][0] * len ; |
2428 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgevmap[edge][0]][1] * len ; | 2451 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgemap[edge][0]][1] * len ; |
2429 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgevmap[edge][0]][2] * len ; | 2452 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgemap[edge][0]][2] * len ; |
2430 int edir = edge / 4 ; | 2453 int edir = edge / 4 ; |
2431 ························································ | 2454 ························································ |
2432 if ( isInProcess( cs[cin
d], edge ) == 0 ) | 2455 if ( isInProcess( cs[cin
d], edge ) == 0 ) |
2433 { | 2456 { |
2434 setInProcessAll(
est, edir ) ; | 2457 setInProcessAll(
est, edir ) ; |
2435 queue->pushQueue
( est, edir ) ; | 2458 queue->pushQueue
( est, edir ) ; |
2436 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; | 2459 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; |
2437 total ++ ; | 2460 total ++ ; |
2438 } | 2461 } |
2439 else | 2462 else |
(...skipping 89 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2529 } | 2552 } |
2530 } | 2553 } |
2531 ················································ | 2554 ················································ |
2532 if ( eind == 3 || eind == -1 ) | 2555 if ( eind == 3 || eind == -1 ) |
2533 { | 2556 { |
2534 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); | 2557 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); |
2535 } | 2558 } |
2536 else· | 2559 else· |
2537 { | 2560 { |
2538 int est[3] ; | 2561 int est[3] ; |
2539 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgevmap[edge][0]][0] * len ; | 2562 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgemap[edge][0]][0] * len ; |
2540 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgevmap[edge][0]][1] * len ; | 2563 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgemap[edge][0]][1] * len ; |
2541 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgevmap[edge][0]][2] * len ; | 2564 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgemap[edge][0]][2] * len ; |
2542 int edir = edge / 4 ; | 2565 int edir = edge / 4 ; |
2543 ························································ | 2566 ························································ |
2544 if ( getEdgeParity( cs[c
ind], edge ) == 1 ) | 2567 if ( getEdgeParity( cs[c
ind], edge ) == 1 ) |
2545 { | 2568 { |
2546 flipParityAll( e
st, edir ) ; | 2569 flipParityAll( e
st, edir ) ; |
2547 queue->pushQueue
( est, edir ) ; | 2570 queue->pushQueue
( est, edir ) ; |
2548 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; | 2571 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; |
2549 total ++ ; | 2572 total ++ ; |
2550 } | 2573 } |
2551 else | 2574 else |
(...skipping 49 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2601 | 2624 |
2602 if ( par == 1 && inp == 0 ) | 2625 if ( par == 1 && inp == 0 ) |
2603 { | 2626 { |
2604 // Intersection edge, hasn't been processed | 2627 // Intersection edge, hasn't been processed |
2605 // Let's start filling | 2628 // Let's start filling |
2606 GridQueue* queue = new GridQueue() ; | 2629 GridQueue* queue = new GridQueue() ; |
2607 int total = 1 ; | 2630 int total = 1 ; |
2608 | 2631 |
2609 // Set to in process | 2632 // Set to in process |
2610 int mst[3] ; | 2633 int mst[3] ; |
2611 » » » » mst[0] = st[0] + vertmap[edgevmap[i][0]][0] * le
n ; | 2634 » » » » mst[0] = st[0] + vertmap[edgemap[i][0]][0] * len
; |
2612 » » » » mst[1] = st[1] + vertmap[edgevmap[i][0]][1] * le
n ; | 2635 » » » » mst[1] = st[1] + vertmap[edgemap[i][0]][1] * len
; |
2613 » » » » mst[2] = st[2] + vertmap[edgevmap[i][0]][2] * le
n; | 2636 » » » » mst[2] = st[2] + vertmap[edgemap[i][0]][2] * len
; |
2614 int mdir = i / 4 ; | 2637 int mdir = i / 4 ; |
2615 setInProcessAll( mst, mdir ) ; | 2638 setInProcessAll( mst, mdir ) ; |
2616 | 2639 |
2617 // Put this edge into queue | 2640 // Put this edge into queue |
2618 queue->pushQueue( mst, mdir ) ; | 2641 queue->pushQueue( mst, mdir ) ; |
2619 | 2642 |
2620 // Queue processing | 2643 // Queue processing |
2621 int nst[3], dir ; | 2644 int nst[3], dir ; |
2622 while ( queue->popQueue( nst, dir ) == 1 ) | 2645 while ( queue->popQueue( nst, dir ) == 1 ) |
2623 { | 2646 { |
(...skipping 59 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2683 } | 2706 } |
2684 } | 2707 } |
2685 ················································ | 2708 ················································ |
2686 if ( eind == 3 || eind == -1 ) | 2709 if ( eind == 3 || eind == -1 ) |
2687 { | 2710 { |
2688 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); | 2711 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); |
2689 } | 2712 } |
2690 else· | 2713 else· |
2691 { | 2714 { |
2692 int est[3] ; | 2715 int est[3] ; |
2693 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgevmap[edge][0]][0] * len ; | 2716 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgemap[edge][0]][0] * len ; |
2694 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgevmap[edge][0]][1] * len ; | 2717 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgemap[edge][0]][1] * len ; |
2695 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgevmap[edge][0]][2] * len ; | 2718 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgemap[edge][0]][2] * len ; |
2696 int edir = edge / 4 ; | 2719 int edir = edge / 4 ; |
2697 ························································ | 2720 ························································ |
2698 if ( isInProcess( cs[cin
d], edge ) == 0 ) | 2721 if ( isInProcess( cs[cin
d], edge ) == 0 ) |
2699 { | 2722 { |
2700 setInProcessAll(
est, edir ) ; | 2723 setInProcessAll(
est, edir ) ; |
2701 queue->pushQueue
( est, edir ) ; | 2724 queue->pushQueue
( est, edir ) ; |
2702 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; | 2725 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; |
2703 total ++ ; | 2726 total ++ ; |
2704 } | 2727 } |
2705 else | 2728 else |
(...skipping 100 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
2806 } | 2829 } |
2807 } | 2830 } |
2808 ················································ | 2831 ················································ |
2809 if ( eind == 3 || eind == -1 ) | 2832 if ( eind == 3 || eind == -1 ) |
2810 { | 2833 { |
2811 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); | 2834 dc_printf("Wrong! this i
s not a consistent sign. %d\n", eind ); |
2812 } | 2835 } |
2813 else· | 2836 else· |
2814 { | 2837 { |
2815 int est[3] ; | 2838 int est[3] ; |
2816 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgevmap[edge][0]][0] * len ; | 2839 » » » » » » » est[0] = cst[cind][0] +
vertmap[edgemap[edge][0]][0] * len ; |
2817 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgevmap[edge][0]][1] * len ; | 2840 » » » » » » » est[1] = cst[cind][1] +
vertmap[edgemap[edge][0]][1] * len ; |
2818 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgevmap[edge][0]][2] * len ; | 2841 » » » » » » » est[2] = cst[cind][2] +
vertmap[edgemap[edge][0]][2] * len ; |
2819 int edir = edge / 4 ; | 2842 int edir = edge / 4 ; |
2820 ························································ | 2843 ························································ |
2821 if ( getEdgeParity( cs[c
ind], edge ) == 1 ) | 2844 if ( getEdgeParity( cs[c
ind], edge ) == 1 ) |
2822 { | 2845 { |
2823 flipParityAll( e
st, edir ) ; | 2846 flipParityAll( e
st, edir ) ; |
2824 queue->pushQueue
( est, edir ) ; | 2847 queue->pushQueue
( est, edir ) ; |
2825 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; | 2848 // dc_printf("Pu
shed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; |
2826 total ++ ; | 2849 total ++ ; |
2827 } | 2850 } |
2828 else | 2851 else |
(...skipping 330 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3159 int parity[12] ; | 3182 int parity[12] ; |
3160 fillEdgeOffsetsNormals( node, st, len, pts, norms, parity ) ; | 3183 fillEdgeOffsetsNormals( node, st, len, pts, norms, parity ) ; |
3161 | 3184 |
3162 numtype zero = 0, one = 1 ; | 3185 numtype zero = 0, one = 1 ; |
3163 for ( int i = 0 ; i < 12 ; i ++ ) | 3186 for ( int i = 0 ; i < 12 ; i ++ ) |
3164 { | 3187 { |
3165 int par = getEdgeParity( node, i ) ; | 3188 int par = getEdgeParity( node, i ) ; |
3166 // Let's check first | 3189 // Let's check first |
3167 if ( par ) | 3190 if ( par ) |
3168 { | 3191 { |
3169 » » » » if ( sgn[ edgevmap[i][0] ] == sgn[ edgevmap[i][1
] ] ) | 3192 » » » » if ( sgn[ edgemap[i][0] ] == sgn[ edgemap[i][1]
] ) |
3170 { | 3193 { |
3171 » » » » » dc_printf("Wrong! Parity: %d Sign: %d %d
\n", parity[i], sgn[ edgevmap[i][0] ], sgn[ edgevmap[i][1] ]); | 3194 » » » » » dc_printf("Wrong! Parity: %d Sign: %d %d
\n", parity[i], sgn[ edgemap[i][0] ], sgn[ edgemap[i][1] ]); |
3172 exit(0) ; | 3195 exit(0) ; |
3173 } | 3196 } |
3174 if ( parity[ i ] == 0 ) | 3197 if ( parity[ i ] == 0 ) |
3175 { | 3198 { |
3176 dc_printf("Wrong! No intersection found.
\n"); | 3199 dc_printf("Wrong! No intersection found.
\n"); |
3177 exit(0) ; | 3200 exit(0) ; |
3178 } | 3201 } |
3179 fwrite( &one, sizeof ( numtype ) , 1, fout ) ; | 3202 fwrite( &one, sizeof ( numtype ) , 1, fout ) ; |
3180 fwrite( &(pts[i]), sizeof( float ), 1, fout ) ; | 3203 fwrite( &(pts[i]), sizeof( float ), 1, fout ) ; |
3181 fwrite( norms[i], sizeof( float ), 3, fout ) ; | 3204 fwrite( norms[i], sizeof( float ), 3, fout ) ; |
3182 | 3205 |
3183 } | 3206 } |
3184 else | 3207 else |
3185 { | 3208 { |
3186 » » » » if ( sgn[ edgevmap[i][0] ] != sgn[ edgevmap[i][1
] ] ) | 3209 » » » » if ( sgn[ edgemap[i][0] ] != sgn[ edgemap[i][1]
] ) |
3187 { | 3210 { |
3188 » » » » » dc_printf("Wrong! Parity: %d Sign: %d %d
\n", parity[i], sgn[ edgevmap[i][0] ], sgn[ edgevmap[i][1] ]); | 3211 » » » » » dc_printf("Wrong! Parity: %d Sign: %d %d
\n", parity[i], sgn[ edgemap[i][0] ], sgn[ edgemap[i][1] ]); |
3189 exit(0) ; | 3212 exit(0) ; |
3190 } | 3213 } |
3191 fwrite ( &zero, sizeof ( numtype ) , 1, fout ); | 3214 fwrite ( &zero, sizeof ( numtype ) , 1, fout ); |
3192 } | 3215 } |
3193 } | 3216 } |
3194 } | 3217 } |
3195 } | 3218 } |
3196 #endif | 3219 #endif |
3197 | 3220 |
3198 | 3221 |
(...skipping 227 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3426 for ( int i = 0 ; i < 3 ; i ++ ) | 3449 for ( int i = 0 ; i < 3 ; i ++ ) |
3427 { | 3450 { |
3428 if ( getFaceEdgeNum( node, i * 2 ) ) | 3451 if ( getFaceEdgeNum( node, i * 2 ) ) |
3429 { | 3452 { |
3430 nface ++ ; | 3453 nface ++ ; |
3431 } | 3454 } |
3432 } | 3455 } |
3433 } | 3456 } |
3434 } | 3457 } |
3435 | 3458 |
3436 void Octree::computeMinimizer( UCHAR* leaf, int st[3], int len, float rvalue[3]
) | 3459 /* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */ |
3437 { | 3460 template<typename _Matrix_Type_> |
3438 » int i ; | 3461 void pseudoInverse(const _Matrix_Type_ &a, |
3439 »······· | 3462 » » » » _Matrix_Type_ &result, |
3440 » // First, gather all edge intersections | 3463 » » » » double epsilon = std::numeric_limits<typename
_Matrix_Type_::Scalar>::epsilon()) |
3441 » float pts[12][3], norms[12][3] ; | 3464 { |
3442 » // fillEdgeIntersections( leaf, st, len, pts, norms ) ; | 3465 » Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU
| |
3443 » int parity[12] ; | 3466 » » » » » » » » » »
» » » » Eigen::ComputeFullV); |
3444 » fillEdgeIntersections( leaf, st, len, pts, norms, parity ) ; | 3467 |
3445 | 3468 » typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(), |
3446 » // Next, construct QEF and minimizer | 3469 » » » » » » » » » »
» » » » » » a.rows()) * |
3447 » float mp[3] = {0, 0, 0 } ; | 3470 » » svd.singularValues().array().abs().maxCoeff(); |
| 3471 |
| 3472 » result = svd.matrixV() * |
| 3473 » » _Matrix_Type_((svd.singularValues().array().abs() > |
| 3474 » » » » » tolerance).select(svd.singularValues(
). |
| 3475 » » » » » » » » » »
array().inverse(), 0)).asDiagonal() * |
| 3476 » » svd.matrixU().adjoint(); |
| 3477 } |
| 3478 |
| 3479 void solve_least_squares(const float halfA[], const float b[], |
| 3480 » » » » » » const float midpoint[], float r
value[]) |
| 3481 { |
| 3482 » /* calculate pseudo-inverse */ |
| 3483 » Eigen::MatrixXf A(3, 3), pinv(3, 3); |
| 3484 » A << halfA[0], halfA[1], halfA[2], |
| 3485 » » halfA[1], halfA[3], halfA[4], |
| 3486 » » halfA[2], halfA[4], halfA[5]; |
| 3487 » pseudoInverse(A, pinv); |
| 3488 |
| 3489 » Eigen::Vector3f b2(b), mp(midpoint), result; |
| 3490 » b2 = b2 + A * -mp; |
| 3491 » result = pinv * b2 + mp; |
| 3492 |
| 3493 » for(int i = 0; i < 3; i++) |
| 3494 » » rvalue[i] = result(i); |
| 3495 } |
| 3496 |
| 3497 void minimize(float rvalue[3], float mp[3], const float pts[12][3], |
| 3498 » » » const float norms[12][3], const int parity[12]) |
| 3499 { |
3448 float ata[6] = { 0, 0, 0, 0, 0, 0 }; | 3500 float ata[6] = { 0, 0, 0, 0, 0, 0 }; |
3449 float atb[3] = { 0, 0, 0 } ; | 3501 float atb[3] = { 0, 0, 0 } ; |
3450 float btb = 0 ; | |
3451 int ec = 0 ; | 3502 int ec = 0 ; |
3452 ········ | 3503 ········ |
3453 » for ( i = 0 ; i < 12 ; i ++ ) | 3504 » for ( int i = 0 ; i < 12 ; i ++ ) |
3454 { | 3505 { |
3455 // if ( getEdgeParity( leaf, i) ) | 3506 // if ( getEdgeParity( leaf, i) ) |
3456 if ( parity[ i ] ) | 3507 if ( parity[ i ] ) |
3457 { | 3508 { |
3458 » » » float* norm = norms[i] ; | 3509 » » » const float* norm = norms[i] ; |
3459 » » » float* p = pts[i] ; | 3510 » » » const float* p = pts[i] ; |
3460 | |
3461 » » » // dc_printf("Norm: %f, %f, %f Pts: %f, %f, %f\n", norm[
0], norm[1], norm[2], p[0], p[1], p[2] ) ; | |
3462 » » »······· | |
3463 » » » /* | |
3464 » » » if ( p[0] < 0 || p[1] < 0 || p[2] < 0 || | |
3465 » » » p[0] > dimen || p[1] > dimen || p[2] > dimen ) | |
3466 » » » { | |
3467 » » » » dc_printf("Bad@!!! At: (%d %d %d) \n", p[0], p[1
], p[2] ) ; | |
3468 » » » } | |
3469 » » » */ | |
3470 | 3511 |
3471 // QEF | 3512 // QEF |
3472 ata[ 0 ] += (float) ( norm[ 0 ] * norm[ 0 ] ); | 3513 ata[ 0 ] += (float) ( norm[ 0 ] * norm[ 0 ] ); |
3473 ata[ 1 ] += (float) ( norm[ 0 ] * norm[ 1 ] ); | 3514 ata[ 1 ] += (float) ( norm[ 0 ] * norm[ 1 ] ); |
3474 ata[ 2 ] += (float) ( norm[ 0 ] * norm[ 2 ] ); | 3515 ata[ 2 ] += (float) ( norm[ 0 ] * norm[ 2 ] ); |
3475 ata[ 3 ] += (float) ( norm[ 1 ] * norm[ 1 ] ); | 3516 ata[ 3 ] += (float) ( norm[ 1 ] * norm[ 1 ] ); |
3476 ata[ 4 ] += (float) ( norm[ 1 ] * norm[ 2 ] ); | 3517 ata[ 4 ] += (float) ( norm[ 1 ] * norm[ 2 ] ); |
3477 ata[ 5 ] += (float) ( norm[ 2 ] * norm[ 2 ] ); | 3518 ata[ 5 ] += (float) ( norm[ 2 ] * norm[ 2 ] ); |
3478 ························ | 3519 ························ |
3479 double pn = p[0] * norm[0] + p[1] * norm[1] + p[2] * nor
m[2] ; | 3520 double pn = p[0] * norm[0] + p[1] * norm[1] + p[2] * nor
m[2] ; |
3480 ························ | 3521 ························ |
3481 atb[ 0 ] += (float) ( norm[ 0 ] * pn ) ; | 3522 atb[ 0 ] += (float) ( norm[ 0 ] * pn ) ; |
3482 atb[ 1 ] += (float) ( norm[ 1 ] * pn ) ; | 3523 atb[ 1 ] += (float) ( norm[ 1 ] * pn ) ; |
3483 atb[ 2 ] += (float) ( norm[ 2 ] * pn ) ; | 3524 atb[ 2 ] += (float) ( norm[ 2 ] * pn ) ; |
3484 | 3525 |
3485 btb += (float) pn * (float) pn ; | |
3486 ························ | |
3487 // Minimizer | 3526 // Minimizer |
3488 mp[0] += p[0] ; | 3527 mp[0] += p[0] ; |
3489 mp[1] += p[1] ; | 3528 mp[1] += p[1] ; |
3490 mp[2] += p[2] ; | 3529 mp[2] += p[2] ; |
3491 ························ | 3530 ························ |
3492 ec ++ ; | 3531 ec ++ ; |
3493 } | 3532 } |
3494 } | 3533 } |
3495 | 3534 |
3496 if ( ec == 0 ) | 3535 if ( ec == 0 ) |
3497 { | 3536 { |
3498 return ; | 3537 return ; |
3499 } | 3538 } |
3500 mp[0] /= ec ; | 3539 mp[0] /= ec ; |
3501 mp[1] /= ec ; | 3540 mp[1] /= ec ; |
3502 mp[2] /= ec ; | 3541 mp[2] /= ec ; |
3503 ········ | 3542 ········ |
3504 » // Solve | 3543 » // Solve least squares |
3505 » BoundingBoxf * box = new BoundingBoxf ; | 3544 » solve_least_squares(ata, atb, mp, rvalue); |
3506 » box->begin.x = (float) st[0] ; | 3545 } |
3507 » box->begin.y = (float) st[1] ; | 3546 |
3508 » box->begin.z = (float) st[2] ; | 3547 void Octree::computeMinimizer( UCHAR* leaf, int st[3], int len, float rvalue[3]
) |
3509 » box->end.x = (float) st[0] + len ; | 3548 { |
3510 » box->end.y = (float) st[1] + len ; | 3549 » // First, gather all edge intersections |
3511 » box->end.z = (float) st[2] + len ; | 3550 » float pts[12][3], norms[12][3] ; |
| 3551 » // fillEdgeIntersections( leaf, st, len, pts, norms ) ; |
| 3552 » int parity[12] ; |
| 3553 » fillEdgeIntersections( leaf, st, len, pts, norms, parity ) ; |
| 3554 |
| 3555 » // Next, construct QEF and minimizer |
| 3556 » float mp[3] = {0, 0, 0}; |
| 3557 » minimize(rvalue, mp, pts, norms, parity); |
3512 ········ | 3558 ········ |
3513 float mat[10] ; | |
3514 calcPoint( ata, atb, btb, mp, rvalue, box, mat ) ; | |
3515 // dc_printf("%f ", error) ; | |
3516 | |
3517 /* Restraining the location of the minimizer */ | 3559 /* Restraining the location of the minimizer */ |
3518 float nh1 = hermite_num * len ; | 3560 float nh1 = hermite_num * len ; |
3519 float nh2 = ( 1 + hermite_num ) * len ; | 3561 float nh2 = ( 1 + hermite_num ) * len ; |
3520 if((mode == DUALCON_MASS_POINT || mode == DUALCON_CENTROID) || | 3562 if((mode == DUALCON_MASS_POINT || mode == DUALCON_CENTROID) || |
3521 ( rvalue[0] < st[0] - nh1 || rvalue[1] < st[1] - nh1 || rvalue[2
] < st[2] - nh1 || | 3563 ( rvalue[0] < st[0] - nh1 || rvalue[1] < st[1] - nh1 || rvalue[2
] < st[2] - nh1 || |
3522 rvalue[0] > st[0] + nh2 || rvalue[1] > st[1] + nh2 || rvalue[2
] > st[2] + nh2 )) | 3564 rvalue[0] > st[0] + nh2 || rvalue[1] > st[1] + nh2 || rvalue[2
] > st[2] + nh2 )) |
3523 { | 3565 { |
3524 if(mode == DUALCON_CENTROID) { | 3566 if(mode == DUALCON_CENTROID) { |
3525 // Use centroids | 3567 // Use centroids |
3526 » » » rvalue[0] = (float) st[0] + len / 2 ; // mp[0] ; | 3568 » » » rvalue[0] = (float) st[0] + len / 2 ; |
3527 » » » rvalue[1] = (float) st[1] + len / 2 ; // mp[1] ; | 3569 » » » rvalue[1] = (float) st[1] + len / 2 ; |
3528 » » » rvalue[2] = (float) st[2] + len / 2 ; // mp[2] ; | 3570 » » » rvalue[2] = (float) st[2] + len / 2 ; |
3529 } | 3571 } |
3530 else { | 3572 else { |
3531 // Use mass point instead | 3573 // Use mass point instead |
3532 rvalue[0] = mp[0] ; | 3574 rvalue[0] = mp[0] ; |
3533 rvalue[1] = mp[1] ; | 3575 rvalue[1] = mp[1] ; |
3534 rvalue[2] = mp[2] ; | 3576 rvalue[2] = mp[2] ; |
3535 } | 3577 } |
3536 } | 3578 } |
3537 } | 3579 } |
3538 | 3580 |
(...skipping 68 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
3607 void Octree::processEdgeWrite( UCHAR* node[4], int depth[4], int maxdep, int dir
) | 3649 void Octree::processEdgeWrite( UCHAR* node[4], int depth[4], int maxdep, int dir
) |
3608 { | 3650 { |
3609 //int color = 0 ; | 3651 //int color = 0 ; |
3610 | 3652 |
3611 int i = 3 ; | 3653 int i = 3 ; |
3612 { | 3654 { |
3613 if ( getEdgeParity( node[i], processEdgeMask[dir][i] ) ) | 3655 if ( getEdgeParity( node[i], processEdgeMask[dir][i] ) ) |
3614 { | 3656 { |
3615 int flip = 0 ; | 3657 int flip = 0 ; |
3616 int edgeind = processEdgeMask[dir][i] ; | 3658 int edgeind = processEdgeMask[dir][i] ; |
3617 » » » if ( getSign( node[i], edgevmap[ edgeind ][ 1 ] ) > 0 ) | 3659 » » » if ( getSign( node[i], edgemap[ edgeind ][ 1 ] ) > 0 ) |
3618 { | 3660 { |
3619 flip = 1 ; | 3661 flip = 1 ; |
3620 } | 3662 } |
3621 ························ | 3663 ························ |
3622 int num = 0 ; | 3664 int num = 0 ; |
3623 { | 3665 { |
3624 int ind[8]; | 3666 int ind[8]; |
3625 if(use_manifold) | 3667 if(use_manifold) |
3626 { | 3668 { |
3627 /* Deprecated | 3669 /* Deprecated |
(...skipping 531 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
4159 le[j] = isLeaf( node, c[j] ) ; | 4201 le[j] = isLeaf( node, c[j] ) ; |
4160 ne[j] = chd[ c[j] ] ; | 4202 ne[j] = chd[ c[j] ] ; |
4161 } | 4203 } |
4162 | 4204 |
4163 edgeProcParity( ne, le, de, depth - 1, cellProcEdgeMask[
i ][ 4 ] ) ; | 4205 edgeProcParity( ne, le, de, depth - 1, cellProcEdgeMask[
i ][ 4 ] ) ; |
4164 } | 4206 } |
4165 } | 4207 } |
4166 ········ | 4208 ········ |
4167 } | 4209 } |
4168 | 4210 |
| 4211 /* definitions for global arrays */ |
| 4212 const int edgemask[3] = {5, 3, 6}; |
| 4213 |
| 4214 const int faceMap[6][4] = { |
| 4215 {4, 8, 5, 9}, |
| 4216 {6, 10, 7, 11}, |
| 4217 {0, 8, 1, 10}, |
| 4218 {2, 9, 3, 11}, |
| 4219 {0, 4, 2, 6}, |
| 4220 {1, 5, 3, 7} |
| 4221 }; |
| 4222 |
| 4223 const int cellProcFaceMask[12][3] = { |
| 4224 {0, 4, 0}, |
| 4225 {1, 5, 0}, |
| 4226 {2, 6, 0}, |
| 4227 {3, 7, 0}, |
| 4228 {0, 2, 1}, |
| 4229 {4, 6, 1}, |
| 4230 {1, 3, 1}, |
| 4231 {5, 7, 1}, |
| 4232 {0, 1, 2}, |
| 4233 {2, 3, 2}, |
| 4234 {4, 5, 2}, |
| 4235 {6, 7, 2} |
| 4236 }; |
| 4237 |
| 4238 const int cellProcEdgeMask[6][5] = { |
| 4239 {0, 1, 2, 3, 0}, |
| 4240 {4, 5, 6, 7, 0}, |
| 4241 {0, 4, 1, 5, 1}, |
| 4242 {2, 6, 3, 7, 1}, |
| 4243 {0, 2, 4, 6, 2}, |
| 4244 {1, 3, 5, 7, 2} |
| 4245 }; |
| 4246 |
| 4247 const int faceProcFaceMask[3][4][3] = { |
| 4248 {{4, 0, 0}, |
| 4249 {5, 1, 0}, |
| 4250 {6, 2, 0}, |
| 4251 {7, 3, 0}}, |
| 4252 {{2, 0, 1}, |
| 4253 {6, 4, 1}, |
| 4254 {3, 1, 1}, |
| 4255 {7, 5, 1}}, |
| 4256 {{1, 0, 2}, |
| 4257 {3, 2, 2}, |
| 4258 {5, 4, 2}, |
| 4259 {7, 6, 2}} |
| 4260 }; |
| 4261 |
| 4262 const int faceProcEdgeMask[3][4][6] = { |
| 4263 {{1, 4, 0, 5, 1, 1}, |
| 4264 {1, 6, 2, 7, 3, 1}, |
| 4265 {0, 4, 6, 0, 2, 2}, |
| 4266 {0, 5, 7, 1, 3, 2}}, |
| 4267 {{0, 2, 3, 0, 1, 0}, |
| 4268 {0, 6, 7, 4, 5, 0}, |
| 4269 {1, 2, 0, 6, 4, 2}, |
| 4270 {1, 3, 1, 7, 5, 2}}, |
| 4271 {{1, 1, 0, 3, 2, 0}, |
| 4272 {1, 5, 4, 7, 6, 0}, |
| 4273 {0, 1, 5, 0, 4, 1}, |
| 4274 {0, 3, 7, 2, 6, 1}} |
| 4275 }; |
| 4276 |
| 4277 const int edgeProcEdgeMask[3][2][5] = { |
| 4278 {{3, 2, 1, 0, 0}, |
| 4279 {7, 6, 5, 4, 0}}, |
| 4280 {{5, 1, 4, 0, 1}, |
| 4281 {7, 3, 6, 2, 1}}, |
| 4282 {{6, 4, 2, 0, 2}, |
| 4283 {7, 5, 3, 1, 2}}, |
| 4284 }; |
| 4285 |
| 4286 const int processEdgeMask[3][4] = { |
| 4287 {3, 2, 1, 0}, |
| 4288 {7, 5, 6, 4}, |
| 4289 {11, 10, 9, 8} |
| 4290 }; |
| 4291 |
| 4292 const int dirCell[3][4][3] = { |
| 4293 {{0, -1, -1}, |
| 4294 {0, -1, 0}, |
| 4295 {0, 0, -1}, |
| 4296 {0, 0, 0}}, |
| 4297 {{-1, 0, -1}, |
| 4298 {-1, 0, 0}, |
| 4299 {0, 0, -1}, |
| 4300 {0, 0, 0}}, |
| 4301 {{-1, -1, 0}, |
| 4302 {-1, 0, 0}, |
| 4303 {0, -1, 0}, |
| 4304 {0, 0, 0}} |
| 4305 }; |
| 4306 |
| 4307 const int dirEdge[3][4] = { |
| 4308 {3, 2, 1, 0}, |
| 4309 {7, 6, 5, 4}, |
| 4310 {11, 10, 9, 8} |
| 4311 }; |
LEFT | RIGHT |