// // Copyright (c) 1993-2022 Robert McNeel & Associates. All rights reserved. // OpenNURBS, Rhinoceros, and Rhino3D are registered trademarks of Robert // McNeel & Associates. // // THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY. // ALL IMPLIED WARRANTIES OF FITNESS FOR ANY PARTICULAR PURPOSE AND OF // MERCHANTABILITY ARE HEREBY DISCLAIMED. // // For complete openNURBS copyright information see . // //////////////////////////////////////////////////////////////// #include "opennurbs.h" #if !defined(ON_COMPILING_OPENNURBS) // This check is included in all opennurbs source .c and .cpp files to insure // ON_COMPILING_OPENNURBS is defined when opennurbs source is compiled. // When opennurbs source is being compiled, ON_COMPILING_OPENNURBS is defined // and the opennurbs .h files alter what is declared and how it is declared. #error ON_COMPILING_OPENNURBS must be defined when compiling opennurbs #endif #include "opennurbs_subd_data.h" ///////////////////////////////////////////////////////// // // Built-in quad subdivision // static double eigenvalue_cos(unsigned int k, unsigned int N) { double a = ((double)k) / ((double)N); if (2 * k > N) { ON_SubDIncrementErrorCount(); ON_ERROR("bogus k"); } double cos_pia = cos(a*ON_PI); // 0 <= k <= N/2 // cos(0) = 1 // cos(pi/6) = sqrt(3)/2 // cos(pi/4) = 1/sqrt(2) // cos(pi/3) = 1/2 // cos(pi/2) = 0 if (0 == k) cos_pia = 1.0; else if (6 * k == N) cos_pia = (0.5*sqrt(3)); else if (4 * k == N) cos_pia = sqrt(0.5); else if (3 * k == N) cos_pia = 0.5; else if (2 * k == N) cos_pia = 0.0; return cos_pia; } static void GetQuadSubdivisionMatrix_eigenvalue_pair( unsigned int k, unsigned int N, double* eigenvalues ) { // 0 <= k <= N, N >= 3 // a = 2*pi*k/N // lambdas = (5 + cos(a) +/- cos(a/2)*sqrt(18 + 2*cos(a)))/16 // Use cos(a) = cos(2pi - a) and cos(pi/2 + a) = cos(pi/2 - a) // to limit all cosine evaluations to 0 <= a <= pi/2. This // insures that eigenvalues that are theoretically equal // get the same value when stored in a double. // c1 = cos(2*pi*k/N) = sign1*cos(pi*k1/N) // c2 = cos(pi*k/N) = sign2*cos(pi*k2/N); double sign1, sign2; unsigned int k1 = 2 * k; if (k1 > N) k1 = 2 * N - k1; // cos(pi*k1/N) = cos(2pi*k/N) if (2 * k1 > N) { sign1 = -1.0; k1 = N - k1; } else { sign1 = 1.0; } // cos(pi*k2/N) = cos(pi*k/N); unsigned int k2 = k; if (2 * k2 > N) { sign2 = -1.0; k2 = N - k2; } else { sign2 = 1.0; } double c1 = sign1*eigenvalue_cos(k1, N); double c2 = sign2*eigenvalue_cos(k2, N); if (fabs(c1 - cos(2.0*(ON_PI*k / ((double)N)))) > 1e-6) { ON_ERROR("bogus c1"); } if (fabs(c2 - cos(ON_PI*k / ((double)N))) > 1e-6) { ON_ERROR("bogus c2"); } double x = 5.0 + c1; double y = 3.0*((4 * k2 == N) ? (((c2 < 0.0) ? -1.0 : 1.0)*sqrt(1.0 + c1 / 9.0)) : c2*sqrt(2.0 + c1 / 4.5)); eigenvalues[0] = 0.0625*(x - y); eigenvalues[1] = 0.0625*(x + y); double lambda0 = (5.0 + c1 - c2*sqrt(18.0 + 2.0 * c1)) / 16.0; double lambda1 = (5.0 + c1 + c2*sqrt(18.0 + 2.0 * c1)) / 16.0; if (fabs(lambda0 - eigenvalues[0]) > 1e-6) { ON_ERROR("bogus lambda0"); } if (fabs(lambda1 - eigenvalues[1]) > 1e-6) { ON_ERROR("bogus lambda1"); } } static unsigned int GetQuadCreaseEigenvalues( unsigned int N, unsigned int R, double* eigenvalues ) { static const double e2[4] = {1.0, 0.5, 0.25, 0.25}; static const double e3[6] = {1.0, 0.5, 0.5, 0.25, 0.25, 0.125}; static const double e4[8] = {1.0, 0.5, 0.5, 0.3475970508005518921819191, 0.25, 0.25, 0.08990294919944810781808094, 0.0625}; static const double e5[10] = {1.0, 0.5, 0.5, 0.4027739218546748178592416, 0.2985112089851678742713826, 0.25, 0.25, 0.06132374806987680407830079, 0.04544938284868830104054729, 0.03661165235168155944989445}; static const double e6[12] = {1.0, 0.5, 0.5, 0.4338607343486338835457207, 0.3500160863488253415435189, 0.2788059129733551205593399, 0.25, 0.25, 0.0428127141360395954023001, 0.0341025407605693744181211, 0.02751214135449768844149263, 0.02387287570313157198721332}; static const double e7[14] = {1.0, 0.5, 0.5, 0.4524883960194854484584697, 0.3868484484678842764124981, 0.3198112959966680476694149, 0.269127266963170621672735, 0.25, 0.25, 0.03111320661766488579086877, 0.02618235305722229063965432, 0.02164520058600606189657112, 0.01850525303440488985059954, 0.0167468245269451691545346}; static const double e8[16] = {1.0, 0.5, 0.5, 0.4643651755137316635248271, 0.4125707561012937086818205, 0.3544821369293089604984253, 0.3012775591854784677914885, 0.2636617247210804425225423, 0.25, 0.25, 0.02347494981551238513891938, 0.02054399860657505933386027, 0.01746052935054695740643632, 0.01500214366764081029526673, 0.01332883274300524623177295, 0.01237889151219760922048721}; static const double e9[18] = {1.0, 0.5, 0.5, 0.4723463441863384015192462, 0.4308226511075866703470591, 0.3815615404738973077099547, 0.3318586781599466572198922, 0.2891905790294515501929928, 0.2602656637115420192311205, 0.25, 0.25, 0.01827951159672544722081085, 0.01645119019440832022485581, 0.01433600966659543228200454, 0.01246857639828100325799948, 0.01104289481022786208695261, 0.01007212033415834999881356, 0.009515058436089155483977101}; static const double e10[20] = {1.0, 0.5, 0.5, 0.4779461323788179580284805, 0.4440893381762303346169582, 0.4024164222230339077518639, 0.3578375190704969729339551, 0.3156351744644154185583903, 0.2809030992590270743668487, 0.2580062768878958596628428, 0.25, 0.25, 0.0146089903173884963168662, 0.01341819015462357521932491, 0.0119416703391074854280824, 0.01053330352465963744605301, 0.009366444788855289841073367, 0.008487506627292569369514451, 0.00788626781457720035829123, 0.00753842240176145199323634}; static const double e11[22] = {1.0, 0.5, 0.5, 0.4820168104737422559487429, 0.453970779385004836126948, 0.4185545871310304361025496, 0.3791682984415772365648197, 0.3396012178854845423282521, 0.3036900507512289303420477, 0.2749846048497542572933016, 0.256424494661140823867932, 0.25, 0.25, 0.01192931172817658658840325, 0.0111241417795889216648019, 0.01007266363842353548275418, 0.009007528743858636629851362, 0.0080675724846343704060705, 0.007308408092049598881127255, 0.00673824807776591199003043, 0.006346184749337779034933979, 0.006117935463105803485445083}; static const double e12[24] = {1.0, 0.5, 0.5, 0.4850641610369383823975897, 0.4614985474347087619667882, 0.4311821255710073423107508, 0.3965315266016385362937548, 0.3603110020217894833299555, 0.3254201240117244451121387, 0.2946777837281643548152523, 0.270614738818649360997862, 0.2552725419984611848187261, 0.25, 0.25, 0.009917592896102293972385845, 0.009355326173828645421819572, 0.008591381125050664705040022, 0.007779755959415404224109989, 0.007026399790425526641978355, 0.006384584778897758789127306, 0.00587150765060411340781648, 0.005485800904827523567735387, 0.005219287163334916237544707, 0.005063378298187826263703993}; static const double e13[26] = {1.0, 0.5, 0.5, 0.4874023599648010128335877, 0.4673500587464579739048183, 0.4411915336267357444040795, 0.4106853584910705583143491, 0.3778828393545584232161755, 0.3449931832638130534338588, 0.3142414753967076886176823, 0.2877286532750699869095651, 0.2672980989570970429233383, 0.2544068297436071576795717, 0.25, 0.25, 0.008370985397992234789274648, 0.007967268997581054543760475, 0.007401542504344500852893572, 0.006777068031025239694881906, 0.006172979526104779585092999, 0.005635704073174505096388671, 0.005185565574477465291827455, 0.00482700980099718390848469, 0.004556832329593394957851468, 0.004369358935986746324441877, 0.0042592717138664641562821}; static const double e14[28] = {1.0, 0.5, 0.5, 0.4892343474885726073240006, 0.4719804834258717874763822, 0.4492298614073128369092644, 0.42228946574606179458229, 0.3927069701810019940266986, 0.362183044731843177171771, 0.3324764840394177067868378, 0.3053079120274184419879482, 0.2822665844803848275730051, 0.2647220004036752431982968, 0.2537393296624997399764198, 0.25, 0.25, 0.00715748880273074983971188, 0.00686054203316051302491066, 0.006434117641714528308420467, 0.005948540274674075908830307, 0.005462450723752333177711865, 0.005014415879728434370349952, 0.004624660494400380852689325, 0.004300690777742151874669176, 0.004042777577568631615465928, 0.003847905739002843063701085, 0.003712201361565622874300741, 0.003632272821743496605377215}; static const double e15[30] = {1.0, 0.5, 0.5, 0.4906956940465353267709958, 0.4757028727019528632086783, 0.455766067824132773843117, 0.4318751128950106685349497, 0.4052160522897933526263431, 0.3771127514616678985622512, 0.3489624347972032671084194, 0.322168625035827813472657, 0.2980747069657881557286016, 0.2779003551379126765050381, 0.2626814932568613771801873, 0.2532135590826550284758983, 0.25, 0.25, 0.00618847384916111889230231, 0.005965420209880330285766977, 0.005638731508127695651873455, 0.005257089756414251201239889, 0.004863929529021519713339058, 0.004490470412801530600934818, 0.004155270492876199683215923, 0.003867086409290046155236813, 0.003628376449228194769988451, 0.003438179362753015717975701, 0.00329408456109496112534435, 0.003193436395811162253984047, 0.00313401097727204912273354}; static const double e16[32] = {1.0, 0.5, 0.5, 0.4918796373911240560348713, 0.4787373551247809344061374, 0.4611427069644219962279705, 0.4398580642032740903608299, 0.4158055761130113066382872, 0.390027416834009389044947, 0.3636412108488202632024345, 0.3377927419491292896466998, 0.3136080620499505782036071, 0.292146813032174258103034, 0.274357897741047252333337, 0.2610378215649526154634098, 0.2527918855047203998444193, 0.25, 0.25, 0.005402764220102485932368158, 0.005232096046270863116682577, 0.004978077778632910171734617, 0.004674960989517055686540183, 0.004355037766598012314251226, 0.00404323357055087285837187, 0.003755831059271643390444508, 0.003501740890995885397890572, 0.003284648000405711892357816, 0.003105035613274500157028433, 0.002961718646983871018171953, 0.002852868988636084124507635, 0.002776644630749646420752988, 0.002731549908274295258929157}; static const double e17[34] = {1.0, 0.5, 0.5, 0.4928519465484388397628774, 0.4812419214171066989423527, 0.4656128029920844325827301, 0.4465606941843988484692585, 0.424812165480439064558305, 0.4011966249932869370793232, 0.3766149590426153432153104, 0.3520057408143481537617202, 0.3283103707443307886219883, 0.3064384408280979168009157, 0.2872343561627713321664122, 0.2714458174511992239759956, 0.2596943795921930616164232, 0.2524484295514478534381793, 0.25, 0.25, 0.004757090297340728139502408, 0.004624358743088481586008327, 0.004424160910175009107610274, 0.004180976088102615001936721, 0.003918959943644192824684854, 0.00365788010922537763040869, 0.003411648832828200475698599, 0.003188720856577044503144021, 0.002993345157921484120166864, 0.002826943464389544626206934, 0.002689264842243817342049787, 0.002579224555426395685829897, 0.00249546001990384341095027, 0.002436674914664392471599904, 0.00240183994959619385922722}; static const double e18[36] = {1.0, 0.5, 0.5, 0.4936600572298639301751713, 0.4833321071081027545791269, 0.4693655637489887341991605, 0.4522328837092091861166156, 0.4325135619763201089802969, 0.4108745595057475305920414, 0.3880478731496155426360755, 0.3648060650769540190542029, 0.3419366335423175665278088, 0.320216108508285892492672, 0.3003846672037622620284949, 0.2831218712479197365270371, 0.2690238699263027046806524, 0.2585822504755439455599121, 0.2521649077398439413425099, 0.25, 0.25, 0.004220179878692891801759286, 0.004115445927936134297961187, 0.003955713186520081732092417, 0.003758739178522429411908854, 0.003542728327855244753694386, 0.003323322098671379512283609, 0.003112217777696630959319223, 0.002917115069157786970353485, 0.002742396862321762198393639, 0.00259003933228738277985616, 0.002460457574771729873773056, 0.002353170917222417960220722, 0.002267275982618210720919093, 0.002201759937123569867153837, 0.002155696524705100825980994, 0.002128362539512277714756394}; static const double e19[38] = {1.0, 0.5, 0.5, 0.494338860883546070532964, 0.4850938522059087731405798, 0.4725442488309184925419464, 0.457069112822645281805701, 0.4391358254141653594145081, 0.4192859786862247557426383, 0.3981190803632466700566333, 0.376274597652746018210246, 0.3544129175605437659005613, 0.333195821913103737936705, 0.3132670521652698526535814, 0.2952334593627969724433445, 0.2796471003777019238459064, 0.2669884954445580774594627, 0.2576512096332567606306801, 0.2519280988436448268724412, 0.25, 0.25, 0.003768993776777150865030121, 0.003685274515452676520870507, 0.003556390829335076853386117, 0.003395405979373806912108798, 0.003216151173333596924639297, 0.003031009581678132504732807, 0.002849721918135655591096795, 0.002679121978037925651288328, 0.002523464094201966948068253, 0.002385003592067608708145546, 0.002264600976431821046188419, 0.002162236332782625743806115, 0.002077400135120119142943681, 0.002009368305901747016666927, 0.0019573850140940428631998, 0.001920778461640462632113888, 0.001899030873473992579157122}; static const double e20[40] = {1.0, 0.5, 0.5, 0.4949144741229649487947415, 0.4865920796256587990369664, 0.4752586491350964339113918, 0.4612216953194131189045514, 0.4448620526473305969505388, 0.4266235443444738895244926, 0.4070009705959143560203706, 0.3865267637057714686536609, 0.3657566942813565517893462, 0.34525503501672240530079, 0.3255795888907262134327553, 0.3072669591378968612633435, 0.2908183749503575990130514, 0.2766862982234843174818616, 0.2652619581329413087246701, 0.2568639586300465015069328, 0.251728246273788614145149, 0.25, 0.25, 0.003386264950190495753910062, 0.003318560306693570214981728, 0.003213497115566963950845717, 0.003080812251285889215138478, 0.002931102745744173241052439, 0.002774195246119185997058834, 0.002618157176972010977012182, 0.002468953239997107683296335, 0.002330561683921314234983495, 0.002205328627589478520754814, 0.002094389446004944194288476, 0.001998057886466699264485282, 0.001916141267292168987804277, 0.001848175283123922040090202, 0.001793588688772172590353523, 0.001751813424471741486263179, 0.001722355238933792682089404, 0.001704837074659703299686351}; static const size_t e_sizeof[] = { sizeof(e2), sizeof(e3), sizeof(e4), sizeof(e5), sizeof(e6), sizeof(e7), sizeof(e8), sizeof(e9), sizeof(e10), sizeof(e11), sizeof(e12), sizeof(e13), sizeof(e14), sizeof(e15), sizeof(e16), sizeof(e17), sizeof(e18), sizeof(e19), sizeof(e20) }; static const double* e[] = { e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e20 }; static size_t e_count = sizeof(e_sizeof) / sizeof(e_sizeof[0]); if (e_count != sizeof(e) / sizeof(e[0])) return ON_SUBD_RETURN_ERROR(0); if (N < 2) return ON_SUBD_RETURN_ERROR(0); unsigned int edex = N - 2; if (edex >= e_count) return ON_SUBD_RETURN_ERROR(0); if (R != e_sizeof[edex] / sizeof(double)) return ON_SUBD_RETURN_ERROR(0); if (nullptr != eigenvalues) { const double* src = e[edex]; const double* src1 = src + R; double* dst = eigenvalues; while (src < src1) { *dst++ = *src++; } } return R; } unsigned int ON_SubDSectorType::GetAllEigenvalues( double* eigenvalues, size_t eigenvalues_capacity ) { if ( 0 == eigenvalues_capacity) eigenvalues = nullptr; else { if (nullptr == eigenvalues) return ON_SUBD_RETURN_ERROR(0); for (size_t i = 0; i < eigenvalues_capacity; i++) eigenvalues[i] = ON_UNSET_VALUE; } if (false == IsValid()) return ON_SUBD_RETURN_ERROR(0); const unsigned int R = PointRingCount(); if ( 0 == R || (nullptr != eigenvalues && eigenvalues_capacity < R)) return ON_SUBD_RETURN_ERROR(0); const ON_SubDVertexTag vertex_tag = VertexTag(); const unsigned int N = EdgeCount(); if (false == ON_SubD::IsValidSectorEdgeCount(vertex_tag, N)) return ON_SUBD_RETURN_ERROR(0); if (ON_SubDVertexTag::Smooth == vertex_tag) { if (nullptr == eigenvalues) { // caller is asking if eigenvalues are available return R; } eigenvalues[0] = 1.0; double a = (3 * N - 7); double b = (N >= 6) ? ((5 * N - 30)*N + 49) : ((5 * N*N + 49) - 30 * N); double c = sqrt(b); double d = 0.125; unsigned int D = N; while (D > 0 && 0 == (D % 2)) { D /= 2; d *= 0.5; } d /= ((double)D); eigenvalues[1] = d*(a + c); eigenvalues[2] = d*(a - c); for (unsigned int k = 1; k < N; k++) { double e2[2]; GetQuadSubdivisionMatrix_eigenvalue_pair(k, N, e2); eigenvalues[2 * k + 1] = e2[0]; eigenvalues[2 * k + 2] = e2[1]; } // return sorted in decreasing order ON_SortDoubleArrayDecreasing(eigenvalues + 1, R - 1); } else if (ON_SubDVertexTag::Crease == vertex_tag) { if (N <= 20) { if (nullptr == eigenvalues) { // caller is asking if eigenvalues are available return R; } if (R != GetQuadCreaseEigenvalues(N, R, eigenvalues) || !(1.0 == eigenvalues[0])) return ON_SUBD_RETURN_ERROR(0); } } if (nullptr == eigenvalues) { // caller is asking if eigenvalues are available return 0; // not an error } if (!(1.0 == eigenvalues[0])) return ON_SUBD_RETURN_ERROR(0); return R; } bool ON_SubDMatrix::EvaluateCosAndSin( unsigned int j, unsigned int n, double* cos_theta, double* sin_theta ) { double cos_angle; double sin_angle; if (n < 1) { if (cos_theta) *cos_theta = ON_DBL_QNAN; if (sin_theta) *sin_theta = ON_DBL_QNAN; return ON_SUBD_RETURN_ERROR(false); } #if defined(ON_DEBUG) #define ON_DEBUG_TEST_SUBD_EVSINCOS #endif #if defined(ON_DEBUG_TEST_SUBD_EVSINCOS) const double theta_dbg = (j*ON_PI)/n; const double cos_dbg = cos(theta_dbg); const double sin_dbg = sin(theta_dbg); #endif double cos_sign = 1.0; while ( j > n) { // cos(a) = -cos(a-pi) // sin(a) = -sin(a-pi) cos_sign *= -1.0; j -= n; } double sin_sign = cos_sign; if ( 0 == j ) { // cos(0) = 1 // sin(0) = 0 cos_angle = cos_sign; sin_angle = 0.0; } else if ( n == j ) { // cos(pi) = -1 // sin(pi) = 0 cos_angle = -cos_sign; sin_angle = 0.0; } else if ( 2*j == n) { // cos(pi/2) = 0 // sin(pi/2) = 1 cos_angle = 0.0; sin_angle = sin_sign; } else { if ( 2*j > n ) { // cos(a) = cos(pi-a) // sin(a) = sin(pi-a) j = n-j; cos_sign *= -1.0; } if (6 * j == n) { // cos(a) = cos(pi/6) // sin(a) = sin(pi/6) cos_angle = cos_sign*0.5*sqrt(3.0); sin_angle = sin_sign*0.5; } else if (4 * j == n) { // cos(a) = cos(pi/4) // sin(a) = sin(pi/4) const double sqrt_half = sqrt(0.5); cos_angle = cos_sign*sqrt_half; sin_angle = sin_sign*sqrt_half; } else if (3 * j == n) { // cos(a) = cos(pi/3) // sin(a) = sin(pi/3) cos_angle = cos_sign*0.5; sin_angle = sin_sign*0.5*sqrt(3.0); } else { const double theta = (j*ON_PI) / n; cos_angle = cos_sign*cos(theta); sin_angle = sin_sign*sin(theta); } } #if defined(ON_DEBUG_TEST_SUBD_EVSINCOS) double e = fabs(cos_dbg - cos_angle) + fabs(sin_dbg - sin_angle); if ( e > 1e-12 ) { ON_SubDIncrementErrorCount(); ON_ERROR("Bug in SinCosEv precision improvment code."); } #endif if (cos_theta) *cos_theta = cos_angle; if (sin_theta) *sin_theta = sin_angle; return true; } unsigned int ON_SubDSectorType::GetSubdivisionMatrix( double** S, size_t S_capacity ) const { if ( S_capacity < 3 || nullptr == S ) return ON_SUBD_RETURN_ERROR(0); const unsigned int R = PointRingCount(); if (R < 3) return ON_SUBD_RETURN_ERROR(0); if ( S_capacity < R ) return ON_SUBD_RETURN_ERROR(0); if (nullptr == S) return ON_SUBD_RETURN_ERROR(0); for (unsigned int i = 0; i < R; i++) { if (nullptr == S[i]) return ON_SUBD_RETURN_ERROR(0); } const unsigned int N = EdgeCount(); if (this->IsSmoothSector() || this->IsDartSector()) { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Smooth and Dart // See comments in ON_SubDSectorType::GetSurfaceEvaluationCoefficients() // for more details on how this case is handled. // Valence 2 has R = 5, N = 5 if (N < 2 || R < 5) return ON_SUBD_RETURN_ERROR(0); /* Parameters: N - [in] vertex valence (>= 3) S - [out] (2N+1)x(2N+1) matrix (4N*N + 4N + 1) doubles S[0,...,2N] = first row S[2N+1,...,4N+1] = second row S[4N+2,...,6N+2] = third row ... S[4N*N + 2N,...,4N*N + 4N] = 2N+1st (last) row */ const double a = 1.0 - 7.0 / (4.0*N); const double b = 3.0 / (2.0*N*N); const double c = 1.0 / (4.0*N*N); const double dart_factor = this->IsDartSector() ? 0.25*cos(2.0*ON_PI/((double)N)) : 0.0; const double d_center = 3.0 / 8.0 + dart_factor; const double d_ring = 3.0 / 8.0 - dart_factor; const double e = 1.0 / 16.0; const double f = 1.0 / 4.0; double* x; double* y; double* x1; if (5 == R) { if (this->IsDartSector()) { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Dart // See comments in ON_SubDSectorType::GetSurfaceEvaluationCoefficients() // for more details on how this case is handled. // // Dart Valence 2 subdivision matrix special case. S[0][0] = 1.0 / 8.0; S[0][1] = 3.0 / 8.0; S[0][2] = 1.0 / 16.0; S[0][3] = 3.0 / 8.0; S[0][4] = 1.0 / 16.0; // creased edge subdivision S[1][0] = 1.0 / 2.0; S[1][1] = 1.0 / 2.0; S[1][2] = 0.0; S[1][3] = 0.0; S[1][4] = 0.0; S[2][0] = 1.0 / 4.0; S[2][1] = 1.0 / 4.0; S[2][2] = 1.0 / 4.0; S[2][3] = 1.0 / 4.0; S[2][4] = 0.0; // smooth edge subdivision // k = 2, theta = 2pi/2 = pi, cos(theta) = cos(pi) = -1 S[3][0] = 1.0 / 8.0; // = 3/4 + cos(theta)/4 S[3][1] = 1.0 / 8.0; // deviates from general case because P0 is a corner of both quads attached to edge(C,P1) S[3][2] = 1.0 / 16.0; S[3][3] = 5.0 / 8.0; // = 3/4 - cos(theta)/4 S[3][4] = 1.0 / 16.0; S[4][0] = 1.0 / 4.0; S[4][1] = 1.0 / 4.0; S[4][2] = 0.0; S[4][3] = 1.0 / 4.0; S[4][4] = 1.0 / 4.0; } else { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Smooth // See comments in ON_SubDSectorType::GetSurfaceEvaluationCoefficients() // for more details on how this case is handled. // // Smooth Valence 2 subdivision matrix special case. S[0][0] = 1.0 / 8.0; S[0][1] = 3.0 / 8.0; S[0][2] = 1.0 / 16.0; S[0][3] = 3.0 / 8.0; S[0][4] = 1.0 / 16.0; S[1][0] = 3.0 / 8.0; S[1][1] = 3.0 / 8.0; S[1][2] = 1.0 / 16.0; S[1][3] = 1.0 / 8.0; S[1][4] = 1.0 / 16.0; S[2][0] = 1.0 / 4.0; S[2][1] = 1.0 / 4.0; S[2][2] = 1.0 / 4.0; S[2][3] = 1.0 / 4.0; S[2][4] = 0.0; S[3][0] = 3.0 / 8.0; S[3][1] = 1.0 / 8.0; S[3][2] = 1.0 / 16.0; S[3][3] = 3.0 / 8.0; S[3][4] = 1.0 / 16.0; S[4][0] = 1.0 / 4.0; S[4][1] = 1.0 / 4.0; S[4][2] = 0.0; S[4][3] = 1.0 / 4.0; S[4][4] = 1.0 / 4.0; } } else { // first row x = S[0]; *x++ = a; x1 = x + R; while (x < x1) { *x++ = b; *x++ = c; } // second row x = S[1]; x1 = x + R - 2; *x++ = d_center; const double* E = x; *x++ = d_ring; *x++ = e; *x++ = e; while (x < x1) *x++ = 0.0; *x++ = e; *x++ = e; // third row x = S[2]; x1 = x + R; *x++ = f; const double* F = x; *x++ = f; *x++ = f; *x++ = f; while (x < x1) *x++ = 0.0; unsigned int k = 0; const unsigned int kmax = R - 1; for (unsigned int i = 3; i < R; i += 2) { x = S[i]; y = S[i + 1]; k = (k + (kmax - 2)) % kmax; *x++ = d_center; *y++ = f; x1 = x + kmax; while (x < x1) { *x++ = E[k]; *y++ = F[k]; k = (k + 1) % kmax; } } if (this->IsDartSector()) { // fix creased edge row x = S[1]; x1 = x + R; *x++ = 0.5; // center vertex *x++ = 0.5; // crease edge vertex while (x < x1) *x++ = 0.0; } } return R; } if ( this->IsCreaseSector() || this->IsCornerSector() ) { // The Catmull-Clark subdivision matrix is singular for a crease // vertex with 2 crease edges and a single face. This case gets // special handling in all ON_SubD code. // S[4][4] = {{3./4, 1./8, 0, 1./8}, {1./2, 1./2, 0, 0}, {1./4, 1./4, 1./4, 1./4}, {1./2, 0, 0, 1./2}}; // eigenvalues[4] = {1, 1./2, 1./4, 1./4}; // eigenvectors[3][4] = {{1, 1, 1, 1}, {0, -1, 0, 1}, {0, 0, 1, 0}}; // (eigenvalue 1/4 has algebraic multiplicity = 2 and geometric multiplicity = 1) if (N < 2 || R < 4) return ON_SUBD_RETURN_ERROR(0); const double w = this->SectorCoefficient(); double* x; double* y; double* x1; if (this->IsCornerSector()) { // first row (1, 0, ..., 0) x = S[0]; x1 = x + R; *x++ =1.0; while (x < x1) *x++ = 0.0; } else { // first row (3/4, 1/8, 0, ..., 0, 1/8) x = S[0]; *x++ = 0.75; *x++ = 0.125; x1 = x + R - 3; while (x < x1) *x++ = 0.0; *x = 0.125; } // second row (1/2, 1/2, 0, ..., 0) x = S[1]; *x++ = 0.5; *x++ = 0.5; x1 = x + R - 2; while (x < x1) *x++ = 0.0; // third row (1/4, 1/4, 1/4, 1/4, 0, ..., 0) x = S[2]; *x++ = 0.25; *x++ = 0.25; *x++ = 0.25; *x++ = 0.25; x1 = x + R-4; while (x < x1) *x++ = 0.0; // last row (1/2, 0, ..., 0, 1/2) x = S[R - 1]; *x++ = 0.5; x1 = x + R - 2; while (x < x1) *x++ = 0.0; *x = 0.5; if (R > 4) { // penultimate row (1/4, 0, ..., 0, 1/4, 1/4, 1/4 ) x = S[R - 2]; *x++ = 0.25; x1 = x + R - 4; while (x < x1) *x++ = 0.0; *x++ = 0.25; *x++ = 0.25; *x = 0.25; // fourth row (3w/4, 1/16,1/16,3(1-w)/4,1/16,1/16, 0, ...., 0) x = S[3]; *x++ = 0.75*w; *x++ = 0.0625; *x++ = 0.0625; *x++ = 0.75*(1.0 - w); *x++ = 0.0625; *x++ = 0.0625; x1 = x + R - 6; while (x < x1) *x++ = 0.0; if (R > 6) { unsigned int k = 0; const unsigned int kmax = R - 1; const double* F = S[2]; const double* E = S[3]; const double F0 = *F++; const double E0 = *E++; for (unsigned int i = 4; i < R - 2; i += 2) { x = S[i]; y = S[i + 1]; *x++ = F0; *y++ = E0; k = (k + (kmax - 2)) % kmax; x1 = x + kmax; while (x < x1) { *x++ = F[k]; *y++ = E[k]; k = (k + 1) % kmax; } } } } return R; } return ON_SUBD_RETURN_ERROR(0); } unsigned int ON_SubDSectorType::GetSubdivisionMatrix( double* S, size_t S_capacity ) const { if ( S_capacity < 9 || nullptr == S ) return ON_SUBD_RETURN_ERROR(0); const unsigned int R = PointRingCount(); if (R < 3) return ON_SUBD_RETURN_ERROR(0); if ( S_capacity < R*R ) return ON_SUBD_RETURN_ERROR(0); double** Srows = new(std::nothrow) double* [R]; if ( nullptr == Srows) return ON_SUBD_RETURN_ERROR(0); Srows[0] = S; for (unsigned int i = 1; i < R; i++) { Srows[i] = Srows[i-1] + R; } const unsigned int rc = GetSubdivisionMatrix(Srows, R); delete[] Srows; return rc; } double ON_SubDSectorType::SubdominantEigenvalue() const { const double rc_error = ON_UNSET_VALUE; const unsigned int R = PointRingCount(); if (R < 3) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int F = FaceCount(); if (F < 1) return ON_SUBD_RETURN_ERROR(rc_error); switch (VertexTag()) { case ON_SubDVertexTag::Unset: break; case ON_SubDVertexTag::Smooth: if (1 == (R % 2)) { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Smooth // See comments in ON_SubDSectorType::GetSurfaceEvaluationCoefficients() // for more details on how this case is handled. // // R = 5 is the valence 2 case. // In this case the evaluation matrix is 5x5 and the eigenvalues are // In this case the evaluation matrix is 5x5 and the eigenvalues are // 1, 1/4, 1/4, -1/4, 1/8. // The multiplicity of 1/4 is two, but -1/4 is also a subdominant eigenvalue. // This function returns 1/4 when R = 5 and tag = Smooth. double cos_2pi_over_F, cos_pi_over_F, lambda; ON_SubDMatrix::EvaluateCosAndSin(2, F, &cos_2pi_over_F, &lambda); ON_SubDMatrix::EvaluateCosAndSin(1, F, &cos_pi_over_F, &lambda); lambda = (5.0 + cos_2pi_over_F + 3.0*cos_pi_over_F*sqrt(2.0*(1.0 + cos_2pi_over_F/9.0)))/16.0; return lambda; } break; case ON_SubDVertexTag::Dart: // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Dart // See comments in ON_SubDSectorType::GetSurfaceEvaluationCoefficients() // for more details on how this case is handled. // // R = 5 is the valence 2 case. // In this case the evaluation matrix is 5x5 and the eigenvalues are // 1, 1/2, 1/4, sqrt(2)/8, -sqrt(2)/8. // All eigenvalues have multiplicity 1, including the subdominant 1/2. // This function returns 1/2 when R = 5 and tag = Dart. if (1 == (R % 2)) return 0.5; break; case ON_SubDVertexTag::Corner: case ON_SubDVertexTag::Crease: if (0 == (R % 2)) return 0.5; break; } return ON_SUBD_RETURN_ERROR(rc_error); } unsigned int ON_SubDSectorType::SubdominantEigenvalueMulitiplicity() const { if (false == IsValid() ) return 0; // Catmull-Clark quad subdivision special cases if (ON_SubDVertexTag::Crease == m_vertex_tag) { if (0 == m_sector_face_count) { // "wire" edges (no faces) // In this case the evaluation matrix is 3x3 and the eigenvalues are // 1, 1/2, 1/4. // All eigenvalues have multiplicity 1, including the subdominant 1/2. return 1; // for the 1/2 } if (1 == m_sector_face_count) { // single face // In this case the evaluation matrix is 4x4 and the eigenvalues are // 1, 1/2, 1/4, 1/4. // The Jordan block for 1/4, 1/4 had a 1 on the super diagonal. return 1; // for the 1/2 } } else if (2 == m_sector_face_count) { // valance 2 special cases if (ON_SubDVertexTag::Smooth == m_vertex_tag) { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Smooth // See comments in ON_SubDSectorType::GetSurfaceEvaluationCoefficients() // for more details on how this case is handled. // // In this case the evaluation matrix is 5x5 and the eigenvalues are // 1, 1/4, 1/4, -1/4, 1/8. // The multiplicity of 1/4 is two, but -1/4 is also a subdominant eigenvalue. // For valence >= 3 the Smooth eigenvalues are 1, 1/4, 1/4, smaller values ... return 2; // for the 1/4, 1/4 pair } if (ON_SubDVertexTag::Dart == m_vertex_tag) { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Dart // See comments in ON_SubDSectorType::GetSurfaceEvaluationCoefficients() // for more details on how this case is handled. // // In this case the evaluation matrix is 5x5 and the eigenvalues are // 1, 1/2, 1/4, sqrt(2)/8, -sqrt(2)/8. // All eigenvalues have multiplicity 1, including the subdominant 1/2. return 1; // for the 1/2 } } // For a Corner vertex with valence = 0 (wire edges), the matrix is 3x3 and the eigenvalues are 1, 1/2, 1/2. // For a Corner vertex with valence = 1 (single face), the matrix is 4x4 and the eigenvalues are 1, 1/2, 1/2, 1/4. // For a Smooth vertex with valence >= 3, the eigenvalues are 1, 1/4, 1/4, smaller values ... // For a Dart vertex with valence >= 3, the eigenvalues are 1, 1/2, 1/2, smaller values ... return 2; } double ON_SubDSectorType::GetSubdominantEigenvectors( double* E1, size_t E1_capacity, double* E2, size_t E2_capacity ) const { const double rc_error = ON_UNSET_VALUE; const double lambda = SubdominantEigenvalue(); if (!(lambda > 0.0 && lambda < 1.0)) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int R = PointRingCount(); if (R < 3) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int F = FaceCount(); if (F < 1) return ON_SUBD_RETURN_ERROR(rc_error); if ( 0 == E1_capacity) E1 = nullptr; else if (E1_capacity < R || nullptr == E1) { return ON_SUBD_RETURN_ERROR(rc_error); } if ( 0 == E2_capacity) E2 = nullptr; else if (E2_capacity < R || nullptr == E2) { return ON_SUBD_RETURN_ERROR(rc_error); } if (nullptr == E1 || nullptr == E2) { // If one of LT0ev or E2 is null, then both must be null. if (nullptr != E1 || nullptr != E2) { return ON_SUBD_RETURN_ERROR(rc_error); } } double y, cos0, cos1, sin0, sin1; switch (VertexTag()) { case ON_SubDVertexTag::Unset: break; case ON_SubDVertexTag::Smooth: if (1 == (R % 2)) { if (nullptr != E1) { // cos0 = cos(2pi * 0/F) // sin0 = sin(2pi * 0/F) // cos1 = cos(2pi * 1/F) // sin1 = sin(2pi * 1/F) cos0 = 1.0; sin0 = 0.0; ON_SubDMatrix::EvaluateCosAndSin(2, F, &cos1, &sin1); //y = 0.5*(3.0 * sqrt((1.0 + cos1 / 9.0) / (1.0 + cos1)) - 1.0); y = (0.25 == lambda) ? -1.0 : (2.0*(4.0*lambda - 1.0) / (1.0 + cos1) - 1.0); E1[0] = 0.0; E2[0] = 0.0; E1[1] = sin0; E2[1] = cos0; unsigned int i = 2; for (;;) { E1[i] = y*(sin0 + sin1); E2[i] = y*(cos0 + cos1); i++; if (i==R) break; E1[i] = sin1; E2[i] = cos1; i++; cos0 = cos1; sin0 = sin1; ON_SubDMatrix::EvaluateCosAndSin(i, F, &cos1, &sin1); // E1[] and E2[] values are symmetric and we could stop halfway and copy // current loop debugged and tested Feb 10, 2015 } } return lambda; } break; case ON_SubDVertexTag::Dart: if (1 == (R % 2)) { if (nullptr != E1) { // cos0 = cos(2pi * 0/F) // sin0 = sin(2pi * 0/F) // cos1 = cos(2pi * 1/F) // sin1 = sin(2pi * 1/F) cos0 = 1.0; sin0 = 0.0; ON_SubDMatrix::EvaluateCosAndSin(2, F, &cos1, &sin1); E1[0] = 0.0; E2[0] = 0.0; E1[1] = sin0; E2[1] = cos0; unsigned int i = 2; for (;;) { E1[i] = sin0 + sin1; E2[i] = cos0 + cos1; i++; if (i==R) break; E1[i] = sin1; E2[i] = cos1; i++; cos0 = cos1; sin0 = sin1; ON_SubDMatrix::EvaluateCosAndSin(i, F, &cos1, &sin1); // E1[] and E2[] values are symmetric and we could stop halfway and copy // current loop debugged and tested Feb 10, 2015 } } return lambda; } break; case ON_SubDVertexTag::Corner: if (0 == (R % 2)) { const unsigned int sector_angle_index = CornerSectorAngleIndex(); const unsigned int M = ON_SubDSectorType::MaximumCornerAngleIndex; //const unsigned int I = ((2*sector_angle_index <= M) ? sector_angle_index : (M-sector_angle_index)); const unsigned int I = 2*((2*sector_angle_index <= M) ? sector_angle_index : (M-sector_angle_index)); // F faces, F-1 interior smooth edges and 2 boundary edges // theta = (i/M*2pi)/F cos1 = 1.0; // = cos(0) sin1 = 0.0; // = sin(0) E1[0] = 0.0; E2[0] = 0.0; unsigned int i = 1; for (;;) { cos0 = cos1; sin0 = sin1; ON_SubDMatrix::EvaluateCosAndSin((i/2 + 1)*I,F*M,&cos1,&sin1); E1[i] = sin0; E2[i] = cos0; i++; E1[i] = sin0+sin1; E2[i] = cos0+cos1; i++; if (R - 1 == i) { E1[i] = sin1; // = sin(pi) E2[i] = cos1; // = cos(pi) break; } } return lambda; } break; case ON_SubDVertexTag::Crease: if (0 == (R % 2)) { if (1 == F) { // one face and 2 boundary edges // E1 = eigenvector with eigenvalue = 1/2 E1[0] = 0.0; // center point coefficients E1[1] = 1.0; // initial boundary edge point coefficient E1[2] = 0.0; // face point coefficient E1[3] = -1.0; // final boundary edge point coefficient // (0,0,1,0) = eigenvector with eigenvalue = 1/4. // E2 is not an eigenvector // S * E2 = (1/4, 1/4, 2/5, 1/4) E2[0] = 1.0; // center point coefficients E2[1] = -2.0; // boundary edge point coefficient E2[2] = -5.0; // face point coefficient E2[3] = -2.0; // boundary edge point coefficient } else { // F faces, F-1 interior smooth edges and 2 boundary edges // theta = pi/F cos1 = 1.0; // = cos(0) sin1 = 0.0; // = sin(0) E1[0] = 0.0; E2[0] = 0.0; unsigned int i = 1; for (;;) { cos0 = cos1; sin0 = sin1; ON_SubDMatrix::EvaluateCosAndSin(i/2 + 1,F,&cos1,&sin1); E1[i] = sin0; E2[i] = cos0; i++; E1[i] = sin0+sin1; E2[i] = cos0+cos1; i++; if (R - 1 == i) { E1[i] = sin1; // = sin(pi) E2[i] = cos1; // = cos(pi) break; } } } return lambda; } break; } return ON_SUBD_RETURN_ERROR(rc_error); } static double ON_SubDSectorType_LimitSurfaceNormalSign( unsigned int R, double sector_angle, const double* L1, const double* L2 ) { // Evaluate a point ring is a counterclockwise circle in the x-y plane with center vertex at (0,0). const double delta_a = sector_angle / (R-1); ON_2dVector T[2] = { ON_2dVector::ZeroVector, ON_2dVector::ZeroVector }; double a = 0.0; for (unsigned int i = 1; i < R; i++, a += delta_a) { double ring_x = cos(a); double ring_y = sin(a); T[0].x += L1[i] * ring_x; T[0].y += L1[i] * ring_y; T[1].x += L2[i] * ring_x; T[1].y += L2[i] * ring_y; } T[0].Unitize(); T[1].Unitize(); // z-component of the cross product should be positive. double z = T[0].x * T[1].y - T[0].y * T[1].x; return z; } double ON_SubDSectorType::SurfaceNormalSign() const { const double rc_error = ON_UNSET_VALUE; if (false == IsValid()) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int R = PointRingCount(); const ON_SubDVertexTag vertex_tag = VertexTag(); ON_SimpleArray buffer; double* LP = buffer.Reserve(3*R); if (nullptr == LP) return ON_SUBD_RETURN_ERROR(rc_error); double* L1 = LP + R; double* L2 = L1 + R; if (R != GetSurfaceEvaluationCoefficients(LP, R, L1, R, L2, R)) return ON_SUBD_RETURN_ERROR(rc_error); double sector_angle; switch (vertex_tag) { case ON_SubDVertexTag::Smooth: case ON_SubDVertexTag::Dart: sector_angle = 2.0*ON_PI; break; case ON_SubDVertexTag::Crease: sector_angle = 0.5*ON_PI; break; case ON_SubDVertexTag::Corner: sector_angle = CornerSectorAngleRadians(); break; default: return ON_SUBD_RETURN_ERROR(rc_error); break; } double z = ON_SubDSectorType_LimitSurfaceNormalSign(R,sector_angle,L1,L2); return z; } bool ON_SubDSectorType::SurfaceEvaluationCoefficientsAvailable() const { if (IsValid()) { // Available as of March 23, 2015 //if (ON_SubDVertexTag::Dart == m_vertex_tag && FaceCount() > 5) //{ // // temporary limit // return false; //} return true; } return false; } unsigned int ON_SubDSectorType::GetSurfaceEvaluationCoefficients( double* LP, size_t LP_capacity, double* L1, size_t L1_capacity, double* L2, size_t L2_capacity ) const { const double lambda = SubdominantEigenvalue(); if (!(lambda > 0.0 && lambda < 1.0)) return ON_SUBD_RETURN_ERROR(0); const unsigned int R = PointRingCount(); if (R < 3) return ON_SUBD_RETURN_ERROR(0); const unsigned int F = FaceCount(); if (F < 1) return ON_SUBD_RETURN_ERROR(0); const double dF = F; if ( 0 == LP_capacity) LP = nullptr; else if (LP_capacity < R || nullptr == LP) { return ON_SUBD_RETURN_ERROR(0); } if ( 0 == L1_capacity) L1 = nullptr; else if (L1_capacity < R || nullptr == L1) { return ON_SUBD_RETURN_ERROR(0); } if ( 0 == L2_capacity) L2 = nullptr; else if (L2_capacity < R || nullptr == L2) { return ON_SUBD_RETURN_ERROR(0); } if (nullptr == L1 || nullptr == L2) { // If one of L1 or L2 is null, then both must be null. if (nullptr != L1 || nullptr != L2) { return ON_SUBD_RETURN_ERROR(0); } } double p, q, y, z, cos0, cos1, sin0, sin1; bool b180degreeCorner = false; switch (VertexTag()) { case ON_SubDVertexTag::Unset: break; case ON_SubDVertexTag::Smooth: if (R >= 5 && 1 == (R % 2)) { if (nullptr != LP) { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Smooth // // General case code works for R = 5 case: LP[5] = {2/7, 2/7, 1/14, 2/7, 1/14}; LP[0] = dF / (5.0 + dF); // center point coefficient y = dF * (dF + 5.0); p = 4.0 / y; // edge point coefficient q = 1.0 / y; // face point coefficient for (unsigned int i = 1; i < R; i++) { LP[i] = p; i++; LP[i] = q; } } if (nullptr != L1) { if (5 == R) { // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Smooth // // Evaluation matrix has eigenvalues 1, 1/4, 1/4, -1/4, 1/8, // so we have a well defined limit point and 3 subdominant eigenvalues. // The 1 eigenspace is 1 dimensional and spanned by (1,1,1,1,1). // The 1/4 eigenspace is 2 dimensional and spanned by (0,1,0,-1,0) and (0,0,1,0,-1). // The -1/4 eigenspace is 1 dimensional and spanned by (2,-1,0,-1,0). // The 1/8 is 1 dimensional and spanned by (1,1,-6,1,-6). // The general case code below will divide by zero when y is evaluated because cos1 = -1. // Turns out there is 1 well defined tangent vector (Q1-Q0) (quad points) // This is a "pseudo" limit tangent L1[0] = 0.0; L1[1] = 1.0; L1[2] = 0.0; L1[3] = -1.0; L1[4] = 0.0; // This is a bonafide limit tangent L2[0] = 0.0; L2[1] = 0.0; L2[2] = 1.0; L2[3] = 0.0; L2[4] = -1.0; } else { // cos0 = cos(2pi * 0/F) // sin0 = sin(2pi * 0/F) // cos1 = cos(2pi * 1/F) // sin1 = sin(2pi * 1/F) cos0 = 1.0; sin0 = 0.0; ON_SubDMatrix::EvaluateCosAndSin(2, F, &cos1, &sin1); y = 2.0* (lambda - 0.125*(3.0 + cos1)) / (1.0 + cos1); L1[0] = 0.0; L1[1] = cos0; L2[0] = 0.0; L2[1] = sin0; unsigned int i = 2; for (;;) { L1[i] = y * (cos0 + cos1); L2[i] = y * (sin0 + sin1); i++; if (i == R) break; L1[i] = cos1; L2[i] = sin1; i++; cos0 = cos1; sin0 = sin1; ON_SubDMatrix::EvaluateCosAndSin(i, F, &cos1, &sin1); // E0[] and E1[] values are symmetric and we could stop halfway and copy // current loop debugged and tested Feb 10, 2015 } } } return R; } break; case ON_SubDVertexTag::Dart: if (1 == (R % 2)) { ON_Matrix Sbuffer; const double*const* S = nullptr; double* E1 = nullptr; double* E2 = nullptr; if (F > 5) { if (false == Sbuffer.Create(R+2,R)) return ON_SUBD_RETURN_ERROR(0); if (R != ON_SubDSectorType::GetSubdivisionMatrix(Sbuffer.m, R)) return ON_SUBD_RETURN_ERROR(0); S = Sbuffer.m; E1 = Sbuffer.m[R]; E2 = Sbuffer.m[R+1]; } const double* termination_tolerances = nullptr; double* eigenvectors[2] = { E1, E2 }; double eigenprecision[2] = { 0 }; double eigenpivots[3] = { 0 }; if (nullptr != LP) { switch (F) { case 2: // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Dart // // Evaluation matrix is 5x5 with has eigenvalues 1, 1/2, 1/4, sqrt(2)/8, - sqrt(2)/8. LP[0] = 8.0 / 31.0; LP[1] = 10.0 / 31.0; LP[2] = 3.0 / 62.0; LP[3] = 10.0 / 31.0; LP[4] = 3.0 / 62.0; break; case 3: LP[0] = 243/19.0; LP[1] = 130/19.0; LP[2] = 1.0; LP[3] = 120/19.0; LP[4] = 29/19.0; //LP[5] = 120/19.0; //LP[6] = 1.0; break; case 4: LP[0] = 752/29.0; LP[1] = 210/29.0; LP[2] = 1.0; LP[3] = 160/29.0; LP[4] = 44/29.0; LP[5] = 180/29.0; //LP[6] = 44/29.0; //LP[7] = 160/29.0; //LP[8] = 1.0; break; case 5: LP[0] = 43.222764462468341719; LP[1] = 7.4578211569974673376; LP[2] = 1.0; LP[3] = 5.0843576860050653249; LP[4] = 1.4771459186191773181; LP[5] = 5.7257510234301278178; LP[6] = 1.5305953634045991926; //LP[7] = 5.7257510234301278178; //LP[8] = 1.4771459186191773181; //LP[9] = 5.0843576860050653249; //LP[10] = 1.0000000000000000000; break; default: { // TODO // Finish general case formula for dominant eigenvector of Transpose(S) // Until I derive the formula, I can get the dominant eigenvector as the kernel of the matrix (S - I). // This calculation has plenty of precision, but it's a less elegant way to get the solution. if (nullptr == S) return ON_SUBD_RETURN_ERROR(0); const double lambda_local = 1.0; const unsigned int lambda_multiplicity = 1; if (lambda_multiplicity != ON_GetEigenvectors(R, S, true, lambda_local, lambda_multiplicity, termination_tolerances, eigenvectors, eigenprecision, eigenpivots)) return ON_SUBD_RETURN_ERROR(0); if (E1[0] < 0.0) { for ( unsigned int i = 0; i < R; i++ ) E1[i] = -E1[i]; } for (unsigned int i = 0; i < R; i++) { if (E1[i] < 0.0) return ON_SUBD_RETURN_ERROR(0); } LP[0] = E1[0]; LP[1] = E1[1]; unsigned int i0 = 2; unsigned int i1 = R-1; while (i0 <= i1) { LP[i0] = (E1[i0] == E1[i1]) ? E1[i0] : 0.5*(E1[i0] + E1[i1]); i0++; i1--; } } break; } if (F >= 3) { // The code below works for all F >= 3. const double* q0 = &LP[2]; double* q1 = &LP[R - 1]; while (q0 < q1) *q1-- = *q0++; double LPsum = 0.0; for (unsigned int i = 0; i < R; i++) LPsum += LP[i]; for (unsigned int i = 0; i < R; i++) LP[i] /= LPsum; } } if (nullptr != L1) { switch (F) { case 2: // April 2019 Dale Lear - Catmull Clark Valence 2 Evaluation: Dart // // Evaluation matrix is 5x5 with has eigenvalues 1, 1/2, 1/4, sqrt(2)/8, - sqrt(2)/8. // The multiplicity of the subdominant eigen value is 1/2 and there is a single well defined // tangent in the limit surface. The other is a "fake" tangent that is close enough. // This is pseudo tangent associated with the eigenvalue 1/4 // (1/4)^n goes to zero much faster than (1/2)^n, but ... // If the 5 points were translated so the limit point was zero // and then a rotation about the origin were applied so bonafide limit tangent // was parallel to the x-axis (L2*(C,P0,Q0,P1,P1) = 0), then, as n goes to infinity, // Q0 - Q1 would be the farthest point from the x-axis. // the 5 L1[0] = 0.0; L1[1] = 0.0; L1[2] = -1.0; L1[3] = 0.0; L1[4] = 1.0; // This is a bonafide limit tangent associated with the eigenvalue 1/2. // I L2[0] = 4.0; // 2^2 L2[1] = 6.0; L2[2] = -1.0; L2[3] = -8.0; L2[4] = -1.0; break; case 3: L1[0] = -9/4.0; // -L2[0]/4 L1[1] = -3.0; L1[2] = 1/2.0; L1[3] = 3.0; L1[4] = 3/4.0; L1[5] = 1; L1[6] = 0; L2[0] = 9.0; // 3^2 L2[1] = 12.0; L2[2] = -3; L2[3] = -16; // -16 // -3, -16 = L2[3] = -16, -3 L2[4] = -3; // -3 L2[5] = 0; L2[6] = 1; break; case 4: L1[0] = -4.0; // -L2[0]/4 L1[1] = -7.0; L1[2] = 1/2.0; L1[3] = 3.0; L1[4] = 3/2.0; L1[5] = 4.0; L1[6] = 1.0; L1[7] = 1; L1[8] = 0; L2[0] = 16.0; // 4^2 L2[1] = 28.0; L2[2] = -3; L2[3] = -16; L2[4] = -7.0; // -3, -16, L2[4], -16, -3 L2[5] = -16; L2[6] = -3; L2[7] = 0; L2[8] = 1; break; case 5: L1[0] = -6.25; // -L2[0]/4 L1[1] = -14.208203932499369089; L1[2] = 0.50000000000000000000; L1[3] = 3.0000000000000000000; L1[4] = 1.9635254915624211362; L1[5] = 5.8541019662496845446; L1[6] = 2.3680339887498948482; L1[7] = 4.6180339887498948482; L1[8] = 1.1545084971874737121; L1[9] = 1; L1[10] = 0; L2[0] = 25.0; // 5^2 L2[1] = 56.832815729997476357; L2[2] = -3; L2[3] = -16; L2[4] = -9.4721359549995793928; L2[5] = -25.888543819998317571; // -3, -16, ..., L2[5], ..., -16, -3 L2[6] = -9.4721359549995793928; L2[7] = -16; // -16 L2[8] = -3; // -3 L2[9] = 0; L2[10] = 1; break; default: { // TODO // Finish general case formula for subdominant eigenvectors of Transpose(S) // Until I derive the formula, I can get the subdominant eigenvectors as the kernel of the matrix (S - 1/2*I). // This calculation has plenty of precision, but it's a less elegant way to get the solution. if (nullptr == S) return ON_SUBD_RETURN_ERROR(0); const double lambda_local = SubdominantEigenvalue(); if (!(lambda_local > 0.0 && lambda_local < 1.0)) return ON_SUBD_RETURN_ERROR(0); const unsigned int lambda_multiplicity = 2; if (lambda_multiplicity != ON_GetEigenvectors(R, S, true, lambda_local, lambda_multiplicity, termination_tolerances, eigenvectors, eigenprecision, eigenpivots)) return ON_SUBD_RETURN_ERROR(0); const double E[2][2] = { { E1[R - 2], E2[R - 2] }, { E1[R - 1], E2[R - 1] } }; const double det = E[0][0] * E[1][1] - E[0][1] * E[1][0]; if (!(0.0 != det)) return 0; const double EtoL[2][2] = { { E[1][1] / det, (0.0 == E[0][1]) ? 0.0 : (-E[0][1] / det) }, { (0.0 == E[1][0]) ? 0.0 : (-E[1][0] / det), E[0][0] / det } }; for (unsigned int i = 0; i < R; i++) { L1[i] = E1[i] * EtoL[0][0] + E2[i] * EtoL[1][0]; L2[i] = E1[i] * EtoL[0][1] + E2[i] * EtoL[1][1]; } const double int_tol = 1e-12; double x; for (unsigned int i = 0; i < R; i++) { x = floor(L1[i]+0.49); if (fabs(x - L1[i]) <= int_tol) L1[i] = x; x = floor(L2[i]+0.49); if (fabs(x - L2[i]) <= int_tol) L2[i] = x; } L1[0] = -0.25*F*F; L1[R - 2] = 1.0; L1[R - 1] = 0.0; L2[0] = F*F; L2[2] = -3.0; L2[3] = -16.0; L2[R-4] = -16.0; L2[R-3] = -3.0; L2[R-2] = 0.0; L2[R-1] = 1.0; // L2[2, ..., R-1] is symmetric about L2[2+R/2] unsigned int i0 = 4; unsigned int i1 = R - 5; while (i0 < i1) { if (L2[i0] != L2[i1]) { x = 0.5*(L2[i0] + L2[i1]); L2[i0] = x; L2[i1] = x; } i0++; i1--; } } break; } } return R; } break; case ON_SubDVertexTag::Corner: if (0 == (R % 2)) { const unsigned int angle_index = CornerSectorAngleIndex(); if ( 36 == angle_index && F > 1) b180degreeCorner = true; // use crease calculation else { // the corner angle alpha is not equal to pi (180 degrees) // so limit tangents will not be parallel // Lpev[] LP[0] = 1.0; // center point coefficient for (unsigned int i = 1; i < R; i++) LP[i++] = 0.0; // all face and edge coefficients = 0 L1[0] = -1.0; // center point coefficient L1[1] = 1.0; // initial boundary edge point coefficient L1[R - 1] = 0.0; // final boundary edge point coefficient // The sign of s is selected to keep the normal pointing "up". double s = (2 * angle_index <= ON_SubDSectorType::MaximumCornerAngleIndex) ? 1.0 : -1.0; L2[0] = -s; // center point coefficient L2[1] = 0.0; // initial boundary edge point coefficient L2[R - 1] = s; // final boundary edge point coefficient for (unsigned int i = 2; i + 1 < R; i++) { L1[i] = 0.0; // all face and interior edge coefficients = 0 L2[i] = 0.0; // all face and interior edge coefficients = 0 } if (36 == angle_index && 1 == F && 4 == R) { // This case has no good solution. L2[0] = -1.0; L2[1] = 0.0; L2[2] = 1.0; L2[3] = 0.0; } // TODO deal with 180 degree "corner case". // L0 and L1 are the same as above. L2 will be similar to the crease case return R; } } if ( false == b180degreeCorner) break; case ON_SubDVertexTag::Crease: // NOTE: In the case when there are 2 crease edes and a single face, // The Catmull-Clark subdivision matrix is singular. if (0 == (R % 2)) { // Lpev[] if (b180degreeCorner) { LP[0] = 1.0; // center point coefficient for (unsigned int i = 1; i < R; i++) LP[i++] = 0.0; // all face and edge coefficients = 0 } else { LP[0] = 2.0 / 3.0; // center point coefficient LP[1] = LP[R - 1] = 1.0 / 6.0; // boundary edge point coefficients for (unsigned int i = 2; i + 1 < R; i++) LP[i++] = 0.0; // all face and interior edge coefficients = 0 } // Using L1[] computes the tangent to the crease curve. This // curve is a uniform cubic spline with control points // (..., final crease edge point, center point, initial crease edge point, ...) L1[ 0] = 0.0; // center point coefficient L1[ 1] = 1.0; // initial boundary edge point coefficient L1[R-1] = -1.0; // final boundary edge point coefficient for (unsigned int i = 2; i+1 < R; i++) L1[i] = 0.0; // all face and interior edge coefficients = 0 // Using L2[] computes a tangent that points from the limit point // into the subdivision surface. if (1 == F) { // This is the special where IsCatmullClarkCreaseOneFaceCase() returns true. // Catmull-Clark subdivision, // center vertex is a crease vertex, // one face and two crease edges. L2[0] = -2.0; // center point coefficients L2[1] = 1.0; // boundary edge point coefficient L2[2] = 0.0; // face point coefficient L2[3] = 1.0; // boundary edge point coefficient } else { // F faces, F-1 interior smooth edges and 2 boundary edges // theta = pi/F cos0 = 1.0; // cos(0) sin0 = 0.0; // sin(0) ON_SubDMatrix::EvaluateCosAndSin(1,F,&cos1,&sin1); // cos1 = cos(pi/F), sin1 = sin(pi/F) z = (1.0 + cos1) / sin1; L2[0] = -sin1; // center point coefficients L2[1] = -0.25*z*(1.0 + 2.0*cos1); // initial boundary edge point coefficient L2[2] = 0.25*(sin0 + sin1); // first interior face point coefficient L2[R-1] = L2[1]; // final boundary edge point coefficient unsigned int i = 3; for (;;) { cos0 = cos1; sin0 = sin1; ON_SubDMatrix::EvaluateCosAndSin(i/2 + 1, F, &cos1, &sin1); L2[i] = sin0; // interior edge point coefficient i++; L2[i] = 0.25*(sin0 + sin1); // interior face point coefficient i++; if ( i >= R-1 ) break; } } return R; } break; } return ON_SUBD_RETURN_ERROR(0); } bool ON_SubDMatrix::IsValid() const { return (m_sector_type.IsValid() && nullptr != m_S && m_R == m_sector_type.PointRingCount()); } bool ON_SubDMatrix::IsValidPointRing( const double* point_ring, size_t point_ring_count, size_t point_ring_stride ) const { return ( point_ring_count >= 4 && point_ring_stride >= 3 && nullptr != point_ring && m_R == point_ring_count); } bool ON_SubDMatrix::EvaluateSubdivisionPoint( unsigned int element_index, const double* point_ring, size_t point_ring_count, size_t point_ring_stride, double subd_point[3] ) const { if ( nullptr == m_S || element_index >= m_R ) return ON_SUBD_RETURN_ERROR(false); if (false == IsValidPointRing(point_ring,point_ring_count,point_ring_stride)) return ON_SUBD_RETURN_ERROR(false); subd_point[0] = 0.0; subd_point[1] = 0.0; subd_point[2] = 0.0; double c; const double* s = m_S[element_index]; const double* s1 = s + m_R; while (s < s1) { c = *s++; subd_point[0] += c*point_ring[0]; subd_point[1] += c*point_ring[1]; subd_point[2] += c*point_ring[2]; point_ring += point_ring_stride; } return true; } bool ON_SubDMatrix::EvaluateSurfacePoint( const double* point_ring, size_t point_ring_count, size_t point_ring_stride, bool bUndefinedNormalIsPossible, ON_SubDSectorSurfacePoint& limit_point ) const { limit_point.m_next_sector_limit_point = nullptr; limit_point.m_sector_face = nullptr; return EvaluateSurfacePoint( point_ring, point_ring_count, point_ring_stride, bUndefinedNormalIsPossible, limit_point.m_limitP, limit_point.m_limitT1, limit_point.m_limitT2, limit_point.m_limitN ); } static bool Internal_GetAlterateTangent( const ON_SubDMatrix& subd_matrix, unsigned int Ldex, size_t point_ring_count, size_t point_ring_stride, const double* point_ring, const double L[3][3], double* alternate_tangent ) { if (point_ring_count < 4 || point_ring_stride < 3 || nullptr == point_ring) return false; if (2 == Ldex) { if ( 4 == subd_matrix.m_R && ON_SubDVertexTag::Crease == subd_matrix.m_sector_type.VertexTag() ) { // valence 2 crease case when crease edges are collinear // F = face point, C = crease vertex point. // Default tangents: // L[1] = ON_3dPoint(point_ring + point_ring_stride) - ON_3dPoint(point_ring + 3 * point_ring_stride); const ON_3dPoint C(point_ring); const ON_3dPoint Q(point_ring + 2 * point_ring_stride); ON_3dVector V(Q - C); if (V.IsNotZero()) { alternate_tangent[0] = V.x; alternate_tangent[1] = V.y; alternate_tangent[2] = V.z; return true; } } } return false; } static bool Internal_GetAlterateNormal( const ON_SubDMatrix& subd_matrix, size_t point_ring_count, size_t point_ring_stride, const double* point_ring, double L[3][3], bool bHaveAlternateL[3], double alternate_normal[3] ) { if (point_ring_count < 4 || point_ring_stride < 3 || nullptr == point_ring) return false; if ( 4 == subd_matrix.m_R && ON_SubDVertexTag::Crease == subd_matrix.m_sector_type.VertexTag() ) { // valence 2 crease case when crease edges are collinear // F = face point, C = crease vertex point. ON_3dVector N(ON_3dVector::ZeroVector); if (false == bHaveAlternateL[2]) { bHaveAlternateL[2] = Internal_GetAlterateTangent(subd_matrix, 2, point_ring_count, point_ring_stride, point_ring, L, L[2]); if (bHaveAlternateL[2]) { ON_3dVector T1(L[1]); ON_3dVector T2(L[2]); N = ON_CrossProduct(T1, T2).UnitVector(); } } if (N.IsZero()) { N = ON_CrossProduct(L[1], L[2]).UnitVector(); if (N.IsZero()) { const ON_3dPoint C(point_ring); const ON_3dPoint P0(point_ring + 1 * point_ring_stride); // end of crease leaving C const ON_3dPoint Q(point_ring + 2 * point_ring_stride); const ON_3dPoint P1(point_ring + 1 * point_ring_stride); // end of crease entering C const ON_3dVector A = (P0 - P1).UnitVector(); const ON_3dVector B = (Q - C).UnitVector(); N = ON_CrossProduct(A, B).UnitVector(); } } if (N.IsNotZero()) { alternate_normal[0] = N.x; alternate_normal[1] = N.y; alternate_normal[2] = N.z; return true; } } return false; } bool ON_SubDMatrix::EvaluateSurfacePoint( const double* point_ring, size_t point_ring_count, size_t point_ring_stride, bool bUndefinedNormalIsPossible, double limit_point[3], double limit_tangent1[3], double limit_tangent2[3], double limit_normal[3] ) const { if (nullptr != limit_point) { limit_point[0] = ON_DBL_QNAN; limit_point[1] = ON_DBL_QNAN; limit_point[2] = ON_DBL_QNAN; } if (nullptr != limit_normal) { limit_normal[0] = ON_DBL_QNAN; limit_normal[1] = ON_DBL_QNAN; limit_normal[2] = ON_DBL_QNAN; } if (nullptr != limit_tangent1) { limit_tangent1[0] = ON_DBL_QNAN; limit_tangent1[1] = ON_DBL_QNAN; limit_tangent1[2] = ON_DBL_QNAN; } if (nullptr != limit_tangent2) { limit_tangent2[0] = ON_DBL_QNAN; limit_tangent2[1] = ON_DBL_QNAN; limit_tangent2[2] = ON_DBL_QNAN; } if (nullptr == m_LP || nullptr == m_L1 || nullptr == m_L2 ) return ON_SUBD_RETURN_ERROR(false); if (false == IsValidPointRing(point_ring,point_ring_count,point_ring_stride)) return ON_SUBD_RETURN_ERROR(false); double x, y, z, c, L[3][3] = { {0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0} }; const double* s[3] = { m_LP, m_L1, m_L2 }; const double* s1 = s[0]+m_R; const double* P = point_ring; while (s[0] < s1) { x = P[0]; y = P[1]; z = P[2]; P += point_ring_stride; c = *s[0]++; L[0][0] += c*x; L[0][1] += c*y; L[0][2] += c*z; c = *s[1]++; L[1][0] += c*x; L[1][1] += c*y; L[1][2] += c*z; c = *s[2]++; L[2][0] += c*x; L[2][1] += c*y; L[2][2] += c*z; } if (nullptr != limit_point) { limit_point[0] = L[0][0]; limit_point[1] = L[0][1]; limit_point[2] = L[0][2]; } if (nullptr == limit_tangent1 && nullptr == limit_tangent2 && nullptr == limit_normal) return true; bool bHaveAlternateL[3] = {}; for (unsigned int Ldex = 1; Ldex < 3; ++Ldex) { if (0.0 == L[Ldex][0] && 0.0 == L[Ldex][1] && 0.0 == L[Ldex][2]) { Internal_GetAlterateTangent(*this, Ldex, point_ring_count, point_ring_stride, point_ring, L, L[Ldex]); bHaveAlternateL[Ldex] = true; } } ON_3dVector T1(L[1]); T1.Unitize(); ON_3dVector T2(L[2]); T2.Unitize(); ON_3dVector N = ON_CrossProduct(T1, T2); N.Unitize(); if (N.IsZero()) { Internal_GetAlterateNormal( *this, point_ring_count, point_ring_stride, point_ring, L, bHaveAlternateL, &N.x ); } if (nullptr != limit_tangent1) { limit_tangent1[0] = T1.x; limit_tangent1[1] = T1.y; limit_tangent1[2] = T1.z; if (0.0 == limit_tangent1[0] && 0.0 == limit_tangent1[1] && 0.0 == limit_tangent1[2]) { if (false == bUndefinedNormalIsPossible) { ON_ERROR("limit_tangent1[0] = zero vector"); bUndefinedNormalIsPossible = true; // one error is sufficient } } } if (nullptr != limit_tangent2) { limit_tangent2[0] = T2.x; limit_tangent2[1] = T2.y; limit_tangent2[2] = T2.z; if (0.0 == limit_tangent2[0] && 0.0 == limit_tangent2[1] && 0.0 == limit_tangent2[2]) { if ( false == bUndefinedNormalIsPossible ) { ON_ERROR("limit_tangent2[0] = zero vector"); bUndefinedNormalIsPossible = true; // one error is sufficient } } } if (nullptr != limit_normal) { limit_normal[0] = N.x; limit_normal[1] = N.y; limit_normal[2] = N.z; if (0.0 == limit_normal[0] && 0.0 == limit_normal[1] && 0.0 == limit_normal[2]) { if ( false == bUndefinedNormalIsPossible ) { ON_ERROR("limit_normal[0] = zero vector"); } } } return true; } class ON_SubDMatrixHashElement { public: ON_SubDMatrix m_sm; class ON_SubDMatrixHashElement* m_next = nullptr; }; static ON_SubDMatrixHashElement* FindMatrixHashElement( const ON_SubDSectorType& sector_type, ON_SubDMatrixHashElement* element_list ) { for (/*empty init*/; nullptr != element_list; element_list = element_list->m_next) { if ( 0 == ON_SubDSectorType::Compare(§or_type,&element_list->m_sm.m_sector_type) ) return element_list; } return nullptr; } static int CompareMatrixHashElement( const ON_SubDMatrixHashElement* a, const ON_SubDMatrixHashElement* b ) { if ( a == b ) return 0; // put nulls at end of list if (nullptr == a) return 1; if (nullptr == b) return -1; if (0 == a->m_sm.m_R || 0 == b->m_sm.m_R) { // at least one matrix is not valid - put invalid matrices at the end if ( a->m_sm.m_R > 0 ) return -1; if ( b->m_sm.m_R > 0 ) return 1; } return ON_SubDSectorType::Compare(&a->m_sm.m_sector_type,&b->m_sm.m_sector_type); } const ON_SubDMatrix& ON_SubDMatrix::FromCache( ON_SubDSectorType sector_type ) { // The ON_SubDMatrix cache is a global resource for the entire application. // This is a thread safe function that manages the ON_SubDMatrix cache. static ON_SleepLock lock; // The sector_type.Hash() is used to find elements in the cache. static ON_SubDMatrixHashElement* hash_table[256] = { 0 }; const unsigned int hash = sector_type.SectorTypeHash(); if (0 == hash) return ON_SUBD_RETURN_ERROR(ON_SubDMatrix::Empty); const unsigned int hash_capacity = (unsigned int)sizeof(hash_table)/sizeof(hash_table[0]); const unsigned int hash_index = hash % hash_capacity; ON_SubDMatrixHashElement* hash_element = FindMatrixHashElement(sector_type,hash_table[hash_index]); if (nullptr != hash_element) return hash_element->m_sm; if ( false == sector_type.IsValid()) return ON_SUBD_RETURN_ERROR(ON_SubDMatrix::Empty); // We need to add this case to the cache. // Lock the cache bool bReturnLock = lock.GetLock(0,30*ON_SleepLock::OneSecond); if (false == bReturnLock) bReturnLock = lock.GetLock(0, ON_SleepLock::OneSecond); // see if another thread made the sector matrix while we waited hash_element = FindMatrixHashElement(sector_type,hash_table[hash_index]); if (nullptr == hash_element) { // These matrices are created one time and then used many times. // The memory used to store the matrices is app workspace memory and // is not leaked. ON_MemoryAllocationTracking disable_tracking(false); // calculate the matrix and add it to the cache unsigned int R = 0; ON_SubDMatrixHashElement* new_hash_element = new(std::nothrow) ON_SubDMatrixHashElement(); if (nullptr != new_hash_element) { R = new_hash_element->m_sm.SetFromSectorType(sector_type); if (0 == R) { delete new_hash_element; new_hash_element = new(std::nothrow) ON_SubDMatrixHashElement(); if (nullptr != new_hash_element) new_hash_element->m_sm.m_sector_type = sector_type; } if (nullptr != new_hash_element) { if (nullptr == hash_table[hash_index]) hash_table[hash_index] = new_hash_element; else { // insert elements used frequently first ON_SubDMatrixHashElement* prev_hash_element = nullptr; hash_element = hash_table[hash_index]; while (nullptr != hash_element) { if (CompareMatrixHashElement(new_hash_element, hash_element) < 0) { new_hash_element->m_next = hash_element; if (nullptr != prev_hash_element) prev_hash_element->m_next = new_hash_element; else hash_table[hash_index] = new_hash_element; break; } prev_hash_element = hash_element; hash_element = hash_element->m_next; if (nullptr == hash_element) { prev_hash_element->m_next = new_hash_element; break; } } } } hash_element = new_hash_element; } } // Unlock the cache if (bReturnLock) lock.ReturnLock(); if ( nullptr != hash_element ) return hash_element->m_sm; return ON_SUBD_RETURN_ERROR(ON_SubDMatrix::Empty); } static unsigned int SetTypeAndValenceFailure() { ON_SubDIncrementErrorCount(); return 0; } unsigned int ON_SubDMatrix::SetFromSectorType( ON_SubDSectorType sector_type ) { m_sector_type = ON_SubDSectorType::Empty; m_R = 0; m_S = nullptr; m_LP = nullptr; m_L1 = nullptr; m_L2 = nullptr; if (false == sector_type.IsValid() ) return SetTypeAndValenceFailure(); const unsigned int R = sector_type.PointRingCount(); if (R < 3) return SetTypeAndValenceFailure(); if (m__max_R < R) m__max_R = 0; double* LP = m__buffer.Reserve(3*R); if ( nullptr == LP ) return SetTypeAndValenceFailure(); double* L1 = LP + R; double* L2 = L1 + R; const bool bLimitEvaluationCoefficientsAvailable = sector_type.SurfaceEvaluationCoefficientsAvailable(); if (bLimitEvaluationCoefficientsAvailable) { if (R != sector_type.GetSurfaceEvaluationCoefficients(LP, R, L1, R, L2, R)) return SetTypeAndValenceFailure(); } else { m__buffer.Zero(); } if (m__max_R < R) { if (!m__S.Create(R, R)) return SetTypeAndValenceFailure(); m__max_R = R; } if (R != sector_type.GetSubdivisionMatrix(m__S.m, R)) return SetTypeAndValenceFailure(); m_sector_type = sector_type; m_R = R; m_S = m__S.m; m_LP = LP; m_L1 = L1; m_L2 = L2; if (false == bLimitEvaluationCoefficientsAvailable) { // This is a temporary lame, slow and approximate hack - good thing it's cached. // WHen I finished calculating the general case for left eigenvectors // of the dart subdivision matrix, this will be deleted. ON_Matrix Sinfinity[2]; Sinfinity[0].Multiply(m__S,m__S); for (unsigned int i = 1; i < 4; i++) { Sinfinity[(i % 2)].Multiply(m__S, Sinfinity[(i+1)%2]); } const double* Sinfinity0 = Sinfinity[1].m[0]; double LPsum = 0.0; for (unsigned int i = 0; i < R; i++) { LP[i] = Sinfinity0[i]; if ( !(LP[i] >= 0.0) ) LP[i] = 0.0; LPsum += LP[i]; } if (fabs(1.0 - LPsum) > 1e-15) { for (unsigned int i = 0; i < R; i++) LP[i] /= LPsum; } } return m_R; } static double TestMatrixReturnValue( double d, double rc ) { if (!(d >= 0.0 && rc >= 0.0)) ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (d > rc) { #if defined(ON_DEBUG) if (d > 0.0001) { // almost certainly a bug ON_SubDIncrementErrorCount(); } #endif rc = d; } return rc; } double ON_SubDMatrix::TestMatrix() const { if (nullptr == m_S || m_R < 3) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if ( false == m_sector_type.IsValid() ) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (m_R != m_sector_type.PointRingCount()) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); const double lambda = m_sector_type.SubdominantEigenvalue(); if (!(lambda > 0.0 && lambda < 1.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); const unsigned int lambda_multiplicity = m_sector_type.SubdominantEigenvalueMulitiplicity(); if ( lambda_multiplicity <= 0 ) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); double rc = 0.0; const bool bTestLimitEvaluation = m_sector_type.SurfaceEvaluationCoefficientsAvailable(); // Eigen information test ON_SimpleArray buffer; double* E1 = buffer.Reserve(5*m_R); double* E2 = E1 + m_R; double* LP = E2 + m_R; double* L1 = LP + m_R; double* L2 = L1 + m_R; double d = m_sector_type.GetSubdominantEigenvectors(E1, m_R, E2, m_R); if (!(d == lambda)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (bTestLimitEvaluation) { if (!m_sector_type.GetSurfaceEvaluationCoefficients(LP, m_R, L1, m_R, L2, m_R)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); } else { for (unsigned int i = 0; i < m_R; i++) { LP[i] = 0.0; L1[i] = 0.0; L2[i] = 0.0; } } // A smooth sector with 2 faces is degenerate and does not have nice eigenvalues and eigenvectors // that give a well defined surface normal. In this case we use a heuristic for the normal. // When bSmoothTwoFaceCase E1 = {0,0,0,0,0}, lengthE1 = 0, lengthE1 = 0. const bool bSmoothTwoFaceCase = (m_sector_type.IsSmoothSector() || m_sector_type.IsDartSector()) && 2u == m_sector_type.FaceCount() && 2u == m_sector_type.EdgeCount() ; const bool bDartTwoFaceCase = (m_sector_type.IsSmoothSector() || m_sector_type.IsDartSector()) && 2u == m_sector_type.FaceCount() && 2u == m_sector_type.EdgeCount() ; const bool bSmoothOrDartTwoFaceCase = bSmoothTwoFaceCase || bDartTwoFaceCase; double lengthE1 = 0.0; double lengthE2 = 0.0; double lengthL1 = 0.0; double lengthL2 = 0.0; for (unsigned int i = 0; i < m_R; i++) { lengthE1 += E1[i] * E1[i]; lengthE2 += E2[i] * E2[i]; lengthL1 += L1[i] * L1[i]; lengthL2 += L2[i] * L2[i]; } if (false == bSmoothOrDartTwoFaceCase && !(lengthE1 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (!(lengthE2 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (!(lengthL1 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (!(lengthL2 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); lengthE1 = sqrt(lengthE1); lengthE2 = sqrt(lengthE2); lengthL1 = sqrt(lengthL1); lengthL2 = sqrt(lengthL2); if (false == bSmoothOrDartTwoFaceCase && !(lengthE1 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (!(lengthE2 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (!(lengthL1 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); if (!(lengthL2 > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); double LPsum = 0.0; double L1sum = 0.0; double L2sum = 0.0; double E1oLP = 0.0; double E2oLP = 0.0; for (unsigned int i = 0; i < m_R; i++) { if (!(LP[i] >= 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); const double* Si = m_S[i]; double x1 = 0.0; double x2 = 0.0; for (unsigned int j = 0; j < m_R; j++) { // E1, E2 should be eigenvectors of S with eigenvalue lambda x1 += Si[j] * E1[j]; x2 += Si[j] * E2[j]; } d = bSmoothOrDartTwoFaceCase ? 0.0 : fabs((x1 - lambda*E1[i])/lengthE1); rc = TestMatrixReturnValue(d,rc); if (!(rc >= 0.0)) break; if (2 == lambda_multiplicity) { d = fabs((x2 - lambda*E2[i])/lengthE2); rc = TestMatrixReturnValue(d, rc); if (!(rc >= 0.0)) break; } if ( false == bTestLimitEvaluation ) continue; if (nullptr != m_LP) { if ( !(m_LP[i] == LP[i]) ) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); } double y1 = 0.0; double y2 = 0.0; for (unsigned int j = 0; j < m_R; j++) { // L1, L2 should be eigenvectors of Transpose(S) with eigenvalue lambda y1 += m_S[j][i] * L1[j]; y2 += m_S[j][i] * L2[j]; } d = (bDartTwoFaceCase && (2u == i || 4u == i)) ? 0.0 : fabs((y1 - lambda*L1[i])/lengthL1); rc = TestMatrixReturnValue(d,rc); if (!(rc >= 0.0)) break; if (2 == lambda_multiplicity) { d = fabs((y2 - lambda*L2[i])/lengthL2); rc = TestMatrixReturnValue(d, rc); if (!(rc >= 0.0)) break; } LPsum += LP[i]; L1sum += L1[i]; L2sum += L2[i]; E1oLP += E1[i] * LP[i]; E2oLP += E2[i] * LP[i]; } if (bTestLimitEvaluation) { while (rc >= 0.0) { // LP coefficients should sum to 1 rc = TestMatrixReturnValue(fabs(1.0 - LPsum), rc); if (!(rc >= 0.0)) break; // E1 and E2 should be orthogonal to LP which means // E0oLP and E1oLP should be zero // // Why? // Set T = Transpose(S). // T*LP = LP (LP is an eigenvector or T with eigenvalue = 1) // S*E = lambda*E (E is an eigenvector of S with eigenvalue lambda != 1) // LP o E1 // = Transpose(LP) * E // = Transpose( T*LP ) * E // = Transpose(LP) * Transpose(T) * E // = Transpose(LP) * S * E // = Transpose(LP) * (lambda E) // = lambda * (Transpose(LP) * E ) // = lambda * LPoE // If LPoE != 0, then lambda = 1, which is not the case. rc = bSmoothOrDartTwoFaceCase ? 0.0 : TestMatrixReturnValue(fabs(E1oLP) / lengthE1, rc); if (!(rc >= 0.0)) break; rc = TestMatrixReturnValue(fabs(E2oLP) / lengthE2, rc); if (!(rc >= 0.0)) break; // L1 and L2 should be orthogonal to (1,1,...,1) which means // L1 L2 coefficients should sum to zero. rc = TestMatrixReturnValue(fabs(L1sum) / lengthL1, rc); if (!(rc >= 0.0)) break; rc = TestMatrixReturnValue(fabs(L2sum) / lengthL2, rc); if (!(rc >= 0.0)) break; break; } // See if L1 and L2 can produce a reasonable normal in the simplest possible case double z = m_sector_type.SurfaceNormalSign(); if (!(z > 0.0)) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); } if (rc >= 0.0) return rc; return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); } double ON_SubDMatrix::TestEvaluation( ON_SubDSectorType sector_type, unsigned int minimum_sector_face_count, unsigned int maximum_sector_face_count, ON_TextLog* text_log ) { ON_SubDVertexTag vertex_tags[] = { ON_SubDVertexTag::Smooth ,ON_SubDVertexTag::Crease ,ON_SubDVertexTag::Corner ,ON_SubDVertexTag::Dart }; const char* vertex_tag_names[sizeof(vertex_tags) / sizeof(vertex_tags[0])] = { "smooth" ,"crease" ,"corner" ,"dart" }; unsigned int corner_sector_angle_index0 = ON_UNSET_UINT_INDEX-1; unsigned int corner_sector_angle_index1 = ON_UNSET_UINT_INDEX; const double corner_sector_angle_radians = (ON_SubDVertexTag::Corner == sector_type.VertexTag()) ? sector_type.CornerSectorAngleRadians() : ON_SubDSectorType::UnsetCornerSectorAngle; ON_SubDVertexTag vertex_tag = sector_type.VertexTag(); size_t subd_type_count = 1; size_t subd_type_index0 = 0; if ( 0 == subd_type_index0 && 2 == subd_type_count) subd_type_count = 1; // tri stuff no ready size_t vertex_tag_count = sizeof(vertex_tags) / sizeof(vertex_tags[0]); size_t vertex_tag_index0 = 0; if (ON_SubDVertexTag::Unset != vertex_tag) { for (size_t vertex_tag_index = vertex_tag_index0; vertex_tag_index < vertex_tag_count; vertex_tag_index++) { if (vertex_tag == vertex_tags[vertex_tag_index]) { vertex_tag_index0 = vertex_tag_index; vertex_tag_count = vertex_tag_index + 1; if (ON_SubDVertexTag::Corner == vertex_tag && ON_SubDSectorType::IsValidCornerSectorAngleRadians(corner_sector_angle_radians) ) { unsigned int angle_dex = sector_type.CornerSectorAngleIndex(); if (angle_dex <= ON_SubDSectorType::MaximumCornerAngleIndex) { corner_sector_angle_index0 = angle_dex; corner_sector_angle_index1 = angle_dex+1; } else { corner_sector_angle_index1 = ON_UNSET_UINT_INDEX; corner_sector_angle_index0 = corner_sector_angle_index1-1; } } break; } } } const unsigned int Fmax = ( 0 == maximum_sector_face_count) ? 20 : maximum_sector_face_count; unsigned int pass_count = 0; unsigned int fail_count = 0; double max_d = 0.0; const unsigned int maximum_fail_count = 10; for (size_t vertex_tag_index = vertex_tag_index0; vertex_tag_index < vertex_tag_count; vertex_tag_index++) { const ON_SubDVertexTag vertex_tag_for_scope = vertex_tags[vertex_tag_index]; const char* sVertexTagName = vertex_tag_names[vertex_tag_index]; unsigned int Fmin = ON_SubDSectorType::MinimumSectorFaceCount(vertex_tag_for_scope); if ( minimum_sector_face_count > Fmin ) Fmin = minimum_sector_face_count; unsigned int angle_i0 = corner_sector_angle_index0; unsigned int angle_i1 = corner_sector_angle_index1; if (ON_SubDVertexTag::Corner == vertex_tag_for_scope && ON_SubDSectorType::UnsetCornerSectorAngle == corner_sector_angle_radians) { angle_i0 = 2; angle_i1 = ON_SubDSectorType::MaximumCornerAngleIndex/2 - 1; } for (unsigned int F = Fmin; F <= Fmax; F++) { for (unsigned int corner_sector_angle_index = angle_i0; corner_sector_angle_index < angle_i1; corner_sector_angle_index++) { double angle_radians = corner_sector_angle_radians; if (ON_SubDVertexTag::Corner == vertex_tag_for_scope && ON_SubDSectorType::UnsetCornerSectorAngle == angle_radians) angle_radians = ON_SubDSectorType::AngleRadiansFromCornerAngleIndex(corner_sector_angle_index); ON_SubDSectorType test_sector_type = ON_SubDSectorType::Create( vertex_tag_for_scope, F, angle_radians); if (false == test_sector_type.SurfaceEvaluationCoefficientsAvailable()) continue; const unsigned int N = test_sector_type.EdgeCount(); const ON_SubDMatrix& SM = ON_SubDMatrix::FromCache(test_sector_type); double d = SM.TestEvaluation(); if (d >= 0.0) { pass_count++; if (d > max_d) max_d = d; } else { fail_count++; } if (nullptr != text_log) { ON_String test_description; if (ON_SubDVertexTag::Corner == vertex_tag_for_scope) test_description.Format("%s, %u faces, %u edges, angle = %u/%u 2pi", sVertexTagName, F, N, corner_sector_angle_index, ON_SubDSectorType::MaximumCornerAngleIndex); else test_description.Format("%s, %u faces, %u edges", sVertexTagName, F, N); if (d >= 0.0) text_log->Print("Test( %s) passed. Deviation = %g\n", (const char*)test_description, d); else text_log->Print("Test( %s ) failed\n", (const char*)test_description); if (SM.m_R > 0) { // Print evaluation coefficients. const ON_TextLogIndent indent1(*text_log); text_log->Print("Limit suface evaluation coefficients:\n"); const ON_TextLogIndent indent2(*text_log); for (unsigned Ldex = 0; Ldex < 3; ++Ldex) { const double* L; ON_String Lid; switch (Ldex) { case 0u: L = SM.m_LP; Lid = "point"; break; case 1u: L = SM.m_L1; Lid = "tangent[1]"; break; case 2u: L = SM.m_L2; Lid = "tangent[2]"; break; default: L = nullptr; break; } if (nullptr == L) continue; char sep[] = { ' ', '=', ' ', 0 }; unsigned termcount = 0; if ( 0.0 == L[0]) text_log->Print("%s", static_cast(Lid)); else { if (1.0 == L[0]) text_log->Print("%s%sV", static_cast(Lid), sep, L[0]); else if (-1.0 == L[0]) text_log->Print("%s%s-V", static_cast(Lid), sep, L[0]); else text_log->Print("%s%s%g*V", static_cast(Lid), sep, L[0]); ++termcount; } bool bPopIndent = false; for (unsigned r = 1u; r < SM.m_R; ++r) { double c = fabs(L[r]); if (0.0 == L[r]) continue; ++termcount; if (termcount >= 2) { c = fabs(c); sep[1] = (L[r] < 0.0) ? '-' : '+'; } if (8 == termcount && SM.m_R > 9 && false == bPopIndent) { text_log->PrintNewLine(); text_log->PushIndent(); bPopIndent = true; } if (1.0 == c) text_log->Print("%s", sep); else text_log->Print("%s%g*", sep, c); if (1u == r % 2u) text_log->Print("E%u", (r+1u)/2u); else text_log->Print("Q%u", r/2u); if (14 == termcount && r < SM.m_R - 2) { text_log->Print(" + ..."); r = SM.m_R - 2; } } text_log->PrintNewLine(); if (bPopIndent) text_log->PopIndent(); } } } if (ON_SubDVertexTag::Corner != vertex_tag_for_scope) break; if (fail_count >= maximum_fail_count) break; } if (fail_count >= maximum_fail_count) break; } if (fail_count >= maximum_fail_count) break; } if (text_log) { text_log->Print("%u tests. %d failed. Maximum deviation = %g\n",pass_count+fail_count, fail_count, max_d); if (fail_count >= maximum_fail_count) text_log->Print("Additional tests canceled because there were too many failures.\n"); } if (fail_count >= maximum_fail_count) return ON_SUBD_RETURN_ERROR(ON_UNSET_VALUE); return max_d; } static void TestPrecision( double d, double& precision ) { if (d > 1e-3) { ON_SubDIncrementErrorCount(); } if (d > precision) precision = d; } double ON_SubDMatrix::TestComponentRing( const ON_SubDComponentPtr* component_ring, size_t component_ring_count ) const { const double rc_error = ON_UNSET_VALUE; if (!m_sector_type.IsValid()) return ON_SUBD_RETURN_ERROR(rc_error); if (component_ring_count < 4 || component_ring_count != m_sector_type.ComponentRingCount()) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int F = m_sector_type.FaceCount(); const unsigned int N = m_sector_type.EdgeCount(); const unsigned int R = m_sector_type.PointRingCount(); if (m_R != R) return ON_SUBD_RETURN_ERROR(rc_error); const ON_SubDVertexTag vertex_tag = m_sector_type.VertexTag(); const unsigned int face_edge_count = m_sector_type.FacetEdgeCount(); const bool bSubdivideFaces = (R == F+N+1 && 4 == face_edge_count); const ON_SubDVertex* center_vertex = component_ring[0].Vertex(); if (nullptr == center_vertex) return ON_SUBD_RETURN_ERROR(rc_error); if (vertex_tag != center_vertex->m_vertex_tag || center_vertex->IsStandard()) return ON_UNSET_VALUE; // not an error - need a round of subdivision before matrix can be used ON_SimpleArray point_ring_array; if ( R != ON_SubD::GetSectorPointRing(false,component_ring_count,component_ring,point_ring_array) ) return ON_SUBD_RETURN_ERROR(rc_error); const ON_3dPoint* vertexR = point_ring_array.Array(); ON_3dPoint Q; if (false == center_vertex->EvaluateCatmullClarkSubdivisionPoint(&Q.x)) return ON_SUBD_RETURN_ERROR(rc_error); double Pdelta = 0.0; double Edelta = 0.0; double Fdelta = 0.0; ON_3dPoint P = ON_3dPoint::Origin; if (false == EvaluateSubdivisionPoint(0, &vertexR[0].x, R,3,P)) return ON_SUBD_RETURN_ERROR(rc_error); double d = P.DistanceTo(Q); TestPrecision(d, Pdelta); for (unsigned int i = 1; i < component_ring_count; i++) { const ON_SubDEdge* edge = component_ring[i].Edge(); if (nullptr == edge) return ON_SUBD_RETURN_ERROR(rc_error); if (false == edge->EvaluateCatmullClarkSubdivisionPoint(&Q.x)) return ON_SUBD_RETURN_ERROR(rc_error); P = ON_3dPoint::Origin; if (false == EvaluateSubdivisionPoint(i, &vertexR[0].x, R,3,P)) return ON_SUBD_RETURN_ERROR(rc_error); double dist_PQ = P.DistanceTo(Q); TestPrecision(dist_PQ, Edelta); if ( !bSubdivideFaces || i+1 >= component_ring_count ) continue; i++; const ON_SubDFace* face = component_ring[i].Face(); if (nullptr == face) return ON_SUBD_RETURN_ERROR(rc_error); if (face_edge_count != face->m_edge_count) return ON_SUBD_RETURN_ERROR(rc_error); if (false == face->EvaluateCatmullClarkSubdivisionPoint(&Q.x)) return ON_SUBD_RETURN_ERROR(rc_error); P = ON_3dPoint::Origin; if (false == EvaluateSubdivisionPoint(i, &vertexR[0].x, R,3,P)) return ON_SUBD_RETURN_ERROR(rc_error); dist_PQ = P.DistanceTo(Q); TestPrecision(dist_PQ, Fdelta); } double max_delta = Pdelta; if (Edelta > max_delta) max_delta = Edelta; if (Fdelta > max_delta) max_delta = Fdelta; return max_delta; } double ON_SubDMatrix::TestEvaluation( const unsigned int subd_recursion_count, ON_SubDSectorIterator sit, ON_SimpleArray& component_ring, ON_SimpleArray< ON_3dPoint >& subd_points, ON_SubDSectorSurfacePoint& limit_point ) { const double rc_error = ON_UNSET_VALUE; limit_point = ON_SubDSectorSurfacePoint::Unset; component_ring.SetCount(0); subd_points.SetCount(0); const ON_SubDSectorType sector_type = ON_SubDSectorType::Create(sit); const unsigned int component_ring_count = ON_SubD::GetSectorComponentRing(sit,component_ring); if ( component_ring_count < 4 || component_ring_count != sector_type.ComponentRingCount()) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int N = ON_SubD::ComponentRingEdgeCount(component_ring_count); const unsigned int F = ON_SubD::ComponentRingFaceCount(component_ring_count); if ( N != sector_type.EdgeCount()) return ON_SUBD_RETURN_ERROR(rc_error); if ( F != sector_type.FaceCount()) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int R = N + F; if (R != sector_type.PointRingCount()) return ON_SUBD_RETURN_ERROR(rc_error); const unsigned int point_ring_stride = 3; ON_3dPoint* point_ring = subd_points.Reserve(R*(subd_recursion_count+1)); const ON_3dPoint* point_ring0 = point_ring; const ON_3dPoint* point_ringEnd = point_ring + R*(subd_recursion_count+1); unsigned int rc = ON_SubD::GetSectorPointRing( true, // bSubdivideIfNeeded, component_ring.Array(), component_ring_count, &point_ring[0].x, R, 3 // point_ring_stride, ); if (rc != R) return ON_SUBD_RETURN_ERROR(rc_error); const ON_SubDVertex* center_vertex = sit.CenterVertex(); const ON_SubDFace* sector_face = sit.CurrentFace(); if (subd_recursion_count > 0) { ON_SubD_FixedSizeHeap fsh[2]; const ON_SubDComponentPtr* element = component_ring.Array(); ON_SimpleArray< ON_SubDComponentPtr > component_ring_buffer; const ON_SubDVertex* v0 = nullptr; for (unsigned int recursion_level = 0; recursion_level < subd_recursion_count; recursion_level++) { const ON_SubDVertex* v1 = ON_SubD::SubdivideSector( v0, element, component_ring_count, fsh[recursion_level%2]); if ( nullptr == v1 || N != v1->m_edge_count || F != v1->m_face_count) return ON_SUBD_RETURN_ERROR(rc_error); if (0 == recursion_level) { if (point_ring0[0].x == v1->m_P[0] && point_ring0[0].y == v1->m_P[1] && point_ring0[0].z == v1->m_P[2]) { // subdivision was required to get initial point ring. v0 = v1; v1 = ON_SubD::SubdivideSector( v0, element, component_ring_count, fsh[recursion_level % 2]); if (nullptr == v1 || N != v1->m_edge_count || F != v1->m_face_count) return ON_SUBD_RETURN_ERROR(rc_error); } } v0 = v1; subd_points.AppendNew() = v1->m_P; for (unsigned i = 0; i < N; i++) { const ON_SubDEdge* e = v1->Edge(i); if ( nullptr == e ) return ON_SUBD_RETURN_ERROR(rc_error); const ON_SubDVertex* v = e->m_vertex[1]; if ( nullptr == v ) return ON_SUBD_RETURN_ERROR(rc_error); subd_points.AppendNew() = v->m_P; if (i < F) { const ON_SubDFace* f = v1->Face(i); if ( nullptr == f ) return ON_SUBD_RETURN_ERROR(rc_error); v = f->Vertex(2); if ( nullptr == v ) return ON_SUBD_RETURN_ERROR(rc_error); subd_points.AppendNew() = v->m_P; } } } } const ON_SubDMatrix& SM = ON_SubDMatrix::FromCache(sector_type); if (false == SM.IsValid() || nullptr == SM.m_S || R != SM.m_R) return ON_SUBD_RETURN_ERROR(rc_error); double maxd = 0.0; double d; // calculate subdivision deviation (should be nearly zero) { const ON_3dPoint* point_ring1 = point_ring0; for (const ON_3dPoint* point_ring2 = point_ring1 + R; point_ring2 < point_ringEnd; point_ring1 += R, point_ring2 += R) { for (unsigned int i = 0; i < R; i++) { ON_3dPoint P(ON_3dPoint::Origin); if ( false == SM.EvaluateSubdivisionPoint(i, &point_ring1[0].x, R,point_ring_stride,P) ) return ON_SUBD_RETURN_ERROR(rc_error); d = P.DistanceTo(point_ring2[i]); if (d > maxd) maxd = d; } point_ring1 = point_ring2; } } // calculate limit point using ON_SubDVertex evaluation code ON_SubDSectorSurfacePoint center_vertex_limit_point = ON_SubDSectorSurfacePoint::Unset; center_vertex->GetSurfacePoint(sector_face,center_vertex_limit_point); // calculate limit point using matrix if ( false == SM.EvaluateSurfacePoint(&point_ring0[0].x, R, point_ring_stride, false, limit_point ) ) return ON_SUBD_RETURN_ERROR(rc_error); // calculate limit point deviation (should be nearly zero) const ON_3dPoint center_vertex_P(center_vertex_limit_point.m_limitP); const ON_3dPoint matrix_P(limit_point.m_limitP); d = center_vertex_P.DistanceTo(matrix_P); if (!ON_IsValid(d)) return ON_SUBD_RETURN_ERROR(rc_error); if (d > maxd) maxd = d; // calculate limit normal deviation (should be nearly zero) const ON_3dVector center_vertex_N(center_vertex_limit_point.m_limitN); const ON_3dVector matrix_N(limit_point.m_limitN); d = (matrix_N - center_vertex_N).Length(); if (!ON_IsValid(d)) return ON_SUBD_RETURN_ERROR(rc_error); if (d > maxd) maxd = d; return maxd; } //#define ON_SUBD_JORDON_DECOMP_CODE #if defined (ON_SUBD_JORDON_DECOMP_CODE) // This code is not required for ON_SubD evaulation. // It was of theoretical interest early in the project. class ON_NumberPolisher { public: ON_NumberPolisher(unsigned int valence, unsigned int R, const double* const * S) : m_N(valence) , m_R(R) , m_S(S) { unsigned int a = valence; while (a > 0 && 0 == a % 2) a /= 2; m_polish_base = 1024.0*a*a; const double polish[] = { 1.0, static_cast(a), static_cast(a*a), static_cast(2 * valence + 1), static_cast(4 * valence - 7), 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0, 25.0, 27.0 }; m_polish_count = sizeof(polish) / sizeof(polish[0]); m_up = std::unique_ptr< double >(new double[m_polish_count + 2 * valence + 1]); double* x = m_up.get(); m_polish = x; m_E = x + m_polish_count; for (size_t i = 0; i < m_polish_count; i++) { x[i] = polish[i]; } } bool PolishEigenvector( double lambda, double* eigenvector, double* eigenprecision ) const; bool PolishGeneralizedEigenvector( double lambda, const double* eigenvector, double* generalized_eigenvector, double* generalized_eigenprecision ) const; double PolishZeroOneValue( double v ) const { const double polish = 1.0; return ON_NumberPolisher::PolishValue(m_polish_01_base, m_polish_01_tol, &polish, 1, v); } double PolishRationalValue( double v ) const { return ON_NumberPolisher::PolishValue(m_polish_base, m_polish_tol, m_polish, m_polish_count, v); } private: const unsigned int m_N = 0; const unsigned int m_R = 0; const double* const * m_S = nullptr; bool m_bTransposeS = false; double m_polish_base = 1024.0; double m_polish_tol = 1.0e-8; const double* m_polish = nullptr; size_t m_polish_count; double m_polish_01_base = 16.0; double m_polish_01_tol = 1.0e-12; double m_polish_01[3]; double* m_E = nullptr; std::unique_ptr< double > m_up; static double PolishValue( const double polish_base, const double polish_tol, const double* polish, const size_t polish_count, double v); }; double ON_NumberPolisher::PolishValue( const double polish_base, const double polish_tol, const double* polish, const size_t polish_count, double v) { if (0.0 == v) return 0.0; if (floor(v) == v) return v; double tol = polish_base*polish_tol; if (tol > 0.015625) tol = 0.015625; for (size_t i = 0; i < polish_count; i++) { double t = polish[i]; if (!(t > 0.0)) continue; const double x = (t * fabs(v))*polish_base; double f = floor(x); if (f == x) return v; double e = x - f; if (e >= 0.5) { f += 1.0; e = f - x; } if (!(e >= 0.0)) continue; if (!(e <= tol)) continue; double polished_v = (f / t) / polish_base; if (0.0 != polished_v && v < 0.0) polished_v = -polished_v; if (v == polished_v) return v; if (fabs(v - polished_v) <= polish_tol) return polished_v; // polish modified v } return v; } bool ON_NumberPolisher::PolishEigenvector( double lambda, double* eigenvector, double* eigenprecision ) const { bool bPolishApplied = false; bool bRationalPolishApplied = false; for (unsigned int j = 0; j < m_R; j++) { double v = PolishZeroOneValue(eigenvector[j]); if (v != eigenvector[j]) { bPolishApplied = true; eigenvector[j] = v; } else { v = PolishRationalValue(eigenvector[j]); if (v != eigenvector[j]) bRationalPolishApplied = true; } m_E[j] = v; } double eprecision = ON_EigenvectorPrecision(m_R, m_S, m_bTransposeS, lambda, eigenvector); if (bRationalPolishApplied) { double Eprecision = ON_EigenvectorPrecision(m_R, m_S, m_bTransposeS, lambda, m_E); if (Eprecision <= eprecision) { for (unsigned int j = 0; j < m_R; j++) { eigenvector[j] = m_E[j]; } eprecision = Eprecision; } } if (eigenprecision) *eigenprecision = eprecision; return (bPolishApplied || bRationalPolishApplied); } bool ON_NumberPolisher::PolishGeneralizedEigenvector( double lambda, const double* xeigenvector, double* generalized_eigenvector, double* generalized_eigenprecision ) const { bool bPolishApplied = false; bool bRationalPolishApplied = false; for (unsigned int j = 0; j < m_R; j++) { double v = PolishZeroOneValue(generalized_eigenvector[j]); if (v != generalized_eigenvector[j]) { bPolishApplied = true; generalized_eigenvector[j] = v; } else { v = PolishRationalValue(generalized_eigenvector[j]); if (v != generalized_eigenvector[j]) bRationalPolishApplied = true; } m_E[j] = v; } double eprecision = ON_MatrixSolutionPrecision(m_R, m_S, m_bTransposeS, lambda, generalized_eigenvector, xeigenvector); if (bRationalPolishApplied) { double Eprecision = ON_MatrixSolutionPrecision(m_R, m_S, m_bTransposeS, lambda, m_E, xeigenvector); if (Eprecision <= eprecision) { for (unsigned int j = 0; j < m_R; j++) { generalized_eigenvector[j] = m_E[j]; } eprecision = Eprecision; } } if (generalized_eigenprecision) *generalized_eigenprecision = eprecision; return (bPolishApplied || bRationalPolishApplied); } static bool GetQuadDominantEigenvector( unsigned int R, const double* const* S, double** eigenvectors, double* eigenprecision, double* eigenpivots ) { // Any matrix whose rows sum to 1 has (1,1,...1) as an eigenvector // and the quad subd matrix is such a matrix. For this matrix, // 1 is the dominant eigenvalue. // The call to ON_GetEigenvectors is still made to insure that // the ON_GetEigenvectors() code is robust enough to calculate this // "easy" eigenvector. const double lambda = 1.0; const unsigned int lambda_muliplicity = 1; const bool bTransposeS = false; // Use easy case as a validation check unsigned int ec = ON_GetEigenvectors(R, S, bTransposeS, lambda, lambda_muliplicity, nullptr, eigenvectors, eigenprecision, eigenpivots); if (1 != ec) return ON_SUBD_RETURN_ERROR(false); double* e1 = eigenvectors[0]; double e1accuracy = 0.0; for (unsigned int i = 0; i < R; i++) { double x = fabs(1.0 - e1[i]); if (x > e1accuracy) e1accuracy = x; e1[i] = 1.0; } if (e1accuracy > 1e-4) { // ON_GetEigenvectors returned junk on the easiest case. // We cannot expect anything good to come out for the // other eigenvalues and S is probably corrupt. return ON_SUBD_RETURN_ERROR(false); } double e1precision = ON_EigenvectorPrecision(R,S,bTransposeS,lambda,e1); if (e1precision > 1e-4) { // S is junk. return ON_SUBD_RETURN_ERROR(false); } if (nullptr != eigenprecision) *eigenprecision = e1precision; return true; } static bool ON_SolveGetGeneralizedEigenvector( const unsigned int N, const double*const* S, bool bTransposeS, double lambda, const double* eigenvector, double* generalized_eigenvector, double* generalized_eigenprecision, double* generalized_eigenpivots ) { if (N < 1 || nullptr == S || !ON_IsValid(lambda) || nullptr == eigenvector) return false; unsigned int i, j, maxi, maxj, n0; double x; const double* src; const double* src1; double* dst; ON_Matrix _M(N, N); double** M = _M.m; ON_SimpleArray< double > _B(N); double* B = _B.Array(); // Xdex records column permutations ON_SimpleArray< unsigned int > _Xdex(N); unsigned int* Xdex = _Xdex.Array(); if (bTransposeS) { // M = Transpose(S - lambda*I) for (i = 0; i < N; i++) { Xdex[i] = i; B[i] = eigenvector[i]; dst = M[i]; for (unsigned int j = 0; j < N; j++) { *dst++ = S[j][i]; } M[i][i] -= lambda; } } else { // M = S - lambda*I for (i = 0; i < N; i++) { Xdex[i] = i; B[i] = eigenvector[i]; src = S[i]; src1 = src + N; dst = M[i]; while (src < src1) { *dst++ = *src++; } M[i][i] -= lambda; } } // reduce system of equations to upper triangular double gpivots[3] = { 0 }; for (n0 = 0; n0 < N-1; n0++) { // find pivot = biggest entry in sub-matrix maxi = n0; maxj = n0; x = 0.0; for (i = n0; i < N; i++) { src = M[i]; src1 = src + N; for ( src += n0; src < src1; src++ ) { if (fabs(*src) > x) { x = fabs(*src); maxi = i; maxj = (unsigned int)(src - M[i]); } } } if (!(x > 0.0)) break; if (0 == n0) { gpivots[0] = x; gpivots[1] = x; } else if (x < gpivots[1]) gpivots[1] = x; else if (x > gpivots[0]) gpivots[0] = x; // put pivot in M[n0][n0] if (n0 != maxi) { // swap rows n0 and maxi dst = M[n0]; M[n0] = M[maxi]; M[maxi] = dst; x = B[n0]; B[n0] = B[maxi]; B[maxi] = x; } if (n0 != maxj) { // swap columns n0 and maxj for (i = 0; i < N; i++) { dst = M[i]; x =dst[n0]; dst[n0] = dst[maxj]; dst[maxj] = x; } j = Xdex[n0]; Xdex[n0] = Xdex[maxj]; Xdex[maxj] = j; } src1 = M[n0] + N; for (i = n0 + 1; i < N; i++) { src = M[n0] + n0; dst = M[i] + n0; x = -(*dst++) / (*src++); if (0.0 != x) { B[i] += x*B[n0]; while ( src < src1 ) *dst++ += x*(*src++); } } } // For debugging, save B ON_SimpleArray< double > savedB; savedB.Append(N, B); // In an ideal world, M[N-1][N-1] and B[N-1] // are zero at this point. // In reality, it should be much smaller than gpivots[1]. x = fabs(M[N - 1][N - 1]); if (fabs(B[N - 1]) > x) x = fabs(B[N - 1]); gpivots[2] = x; if (nullptr != generalized_eigenpivots) { generalized_eigenpivots[0] = gpivots[0]; generalized_eigenpivots[1] = gpivots[1]; generalized_eigenpivots[2] = gpivots[2]; } i = N - 1; B[i] = 0.0; src1 = B + i; while (i > 0) { src = B + i; const double* Mi = M[i-1] + i; x = 0.0; while (src < src1) { x += (*src++)*(*Mi++); } i--; B[i] = (B[i] - x) / M[i][i]; } // solution with column permutation applied is in B[] for (i = 0; i < N; i++) generalized_eigenvector[Xdex[i]] = B[i]; if (nullptr != generalized_eigenprecision) *generalized_eigenprecision = ON_MatrixSolutionPrecision(N, S, bTransposeS, lambda, generalized_eigenvector, eigenvector); return true; } static unsigned int ON_GetGeneralizedEigenvectors( const ON_NumberPolisher& polish, const unsigned int N, const double*const* M, bool bTransposeM, double lambda, unsigned int lambda_multiplicity, double** eigenvectors, double* eigenprecision, double* eigenpivots ) { if (N < 1) return ON_SUBD_RETURN_ERROR(0); if (nullptr == M) return ON_SUBD_RETURN_ERROR(0); if (!(lambda == lambda && ON_IsValid(lambda))) return ON_SUBD_RETURN_ERROR(0); if (lambda_multiplicity < 1 || lambda_multiplicity > N) return ON_SUBD_RETURN_ERROR(0); // Get eigenvectors double epivots[3] = { 0 }; unsigned int ec = ON_GetEigenvectors(N, M, bTransposeM, lambda, lambda_multiplicity, nullptr, eigenvectors, eigenprecision, epivots); if (0 == ec || ec > lambda_multiplicity) return ON_SUBD_RETURN_ERROR(0); if (ec > 0 && ec <= lambda_multiplicity) { for (unsigned int n = 0; n < ec; n++) { double ep = 0.0; //polish.PolishEigenvector(lambda, eigenvectors[n], &ep); if (eigenprecision) eigenprecision[n] = ep; } if (2 == lambda_multiplicity && 1 == ec) { // This case happens for quad subdivision matrices at crease vertices. // Get the generalized eigenvector by solving (M - lambda*I)*g = e double gpivots[3] = { 0 }; double gep = 0.0; bool rc = ON_SolveGetGeneralizedEigenvector( N, M, bTransposeM, lambda, eigenvectors[0], eigenvectors[1], &gep, gpivots ); if (rc) { //polish.PolishGeneralizedEigenvector(lambda, eigenvectors[0], eigenvectors[1], &gep); ec++; } if (eigenprecision) eigenprecision[1] = gep; } } return ec; } /* Description: Calculate the Jordon decomposition for the subdivision matrix. Parameters: eigenvalues - [in] R eigenvalues of S sorted in decreasing order. 1.0 = eigenvalues[0]. super_diagonal - [out] If super_diagonal is not nullpter, then the R-1 values of the Jordon block super diagonal are returned here (zeros and ones). eigenvectors - [out] If eigenvectors is not nullptr, then R generalized eigenvectors are returned so that if lambda = eigenvalues[i] and V = eigenvectors[i], then lambda*V = S * V. eigenprecision - [out] If eigenprecision is not nullptr, then eigenprecision[i] reports on the accuracy of the corresponding eigenvector. eigenprecision[i] = maximum value fabs(lambda*E[j] - (S*E[i])[j]), 0 <= j < 2N+1. eigenpivots - [out] If eigenprecision is not nullptr, then eigenpivots[i] reports on the pivots found when calculating the eigenvectors. eigenpivots[i][0] = maximum nonzero pivot eigenpivots[i][1] = minimum nonzero pivot eigenpivots[i][2] = minimum "zero" pivot E - [out] If E is not nullptr, then E[][] is the NxN matrix whose columns are the eigenvectors of S. Einverse - [out] If Einverse is not nullptr, then Einvers[][] is the inverse of the NxN matrix whose columns are the eigenvectors of S. Returns: R > 0: The size of the subdivision matrix is RxR. 0: failure. Remarks: S = E * J * Inverse(E), where J is the Jordon matrix with diagonal = eigenvalues[], and super diagonal = super_diagonal[], and zeros elsewhere. */ unsigned int ON_SubDMatrix_GetJordonDecomposition( const ON_SubDMatrix& subd_matrix, const double* eigenvalues, double* super_diagonal, double** eigenvectors, double* eigenprecision, double** eigenpivots, double** E, double** Einverse ) { unsigned int R = subd_matrix.m_R; if (R < 3 ) return ON_SUBD_RETURN_ERROR(0); const bool bTransposeS = false; const double* const* S = subd_matrix.m_S; if (nullptr == S) return ON_SUBD_RETURN_ERROR(0); if (nullptr == eigenvalues) return ON_SUBD_RETURN_ERROR(0); if (nullptr == S[0]) return ON_SUBD_RETURN_ERROR(0); if (!(1.0 == eigenvalues[0])) return ON_SUBD_RETURN_ERROR(0); for (unsigned int i = 1; i < R; i++) { if (nullptr == S[i]) return ON_SUBD_RETURN_ERROR(0); if (!(eigenvalues[i] <= eigenvalues[i-1])) return ON_SUBD_RETURN_ERROR(0); } const bool bNeed_Einverse = (nullptr != Einverse); const bool bNeed_eigenvectors = (bNeed_Einverse || nullptr != E || nullptr != eigenprecision || nullptr != eigenpivots || nullptr != super_diagonal); ON_Matrix _eigenvectors; if (bNeed_eigenvectors) { if (nullptr == eigenvectors) { _eigenvectors.Create(R, R); eigenvectors = _eigenvectors.m; } } ON_SimpleArray< double > _super_diagonal; if (bNeed_eigenvectors && nullptr == super_diagonal) super_diagonal = _super_diagonal.Reserve(R - 1); //////////////////////////////////////////////////////////////////////// // // Get eigenvectors[] // if (nullptr != eigenvectors) { const ON_NumberPolisher polish(R/2, R, S); ON_SimpleArray _eprecision; if (nullptr == eigenprecision) eigenprecision = _eprecision.Reserve(R); if (false == GetQuadDominantEigenvector(R, S, eigenvectors, eigenprecision, eigenpivots ? eigenpivots[0] : nullptr)) return ON_SUBD_RETURN_ERROR(0); unsigned int eigenvector_count = 1; double lambda = 1.0; while (eigenvector_count < R) { if (!(eigenvalues[eigenvector_count] < lambda)) return ON_SUBD_RETURN_ERROR(0); lambda = eigenvalues[eigenvector_count]; unsigned int lambda_multiplicity = 1; while (eigenvector_count + lambda_multiplicity < R && lambda == eigenvalues[eigenvector_count + lambda_multiplicity]) lambda_multiplicity++; double epivots[3] = { 0 }; unsigned int ec = ON_GetGeneralizedEigenvectors(polish, R, S, bTransposeS, lambda, lambda_multiplicity, eigenvectors + eigenvector_count, eigenprecision + eigenvector_count, epivots); if (ec != lambda_multiplicity) { return ON_SUBD_RETURN_ERROR(0); } for (unsigned int n = eigenvector_count; n < eigenvector_count + lambda_multiplicity; n++) { if (nullptr != eigenpivots) { eigenpivots[n][0] = epivots[0]; eigenpivots[n][1] = epivots[1]; eigenpivots[n][2] = epivots[2]; } } eigenvector_count += lambda_multiplicity; } //////////////////////////////////////////////////////////////////////// // // Get E[] // if (nullptr != E) { for (unsigned int i = 0; i < R; i++) { const double* Ecol = eigenvectors[i]; for (unsigned int j = 0; j < R; j++) { E[j][i] = Ecol[j]; } } } if (bNeed_Einverse) { ON_Matrix _Einverse; if (eigenvectors == _eigenvectors.m && nullptr == Einverse) { // use eigenvectors memory _Einverse.Create(R, R, eigenvectors, false); _Einverse.Transpose(); Einverse = _Einverse.m; } else { if (nullptr == Einverse) { _Einverse.Create(R, R); Einverse = _Einverse.m; } else _Einverse.Create(R, R, Einverse, false); for (unsigned int i = 0; i < R; i++) { const double* Ecoli = eigenvectors[i]; for (unsigned int j = 0; j < R; j++) { _Einverse.m[j][i] = Ecoli[j]; } } } if (false == _Einverse.Invert(0.0)) return ON_SUBD_RETURN_ERROR(0); } } return R; } #endif