//---------------------------------------------------------------------------- // EgalTech 2023-2023 //---------------------------------------------------------------------------- // File : SurfAux.cpp Data : 09.08.23 Versione : // Contenuto : Implementazione di alcune funzioni di utilità per le Superfici. // // // // Modifiche : 09.08.23 DB Creazione modulo. // // //---------------------------------------------------------------------------- //--------------------------- Include ---------------------------------------- #include "stdafx.h" #include "CurveAux.h" #include "GeoConst.h" #include "CurveLine.h" #include "CurveArc.h" #include "CurveBezier.h" #include "CurveComposite.h" #include "/EgtDev/Include/EgtPointerOwner.h" #include "/EgtDev/Include/EGkSurf.h" #include "/EgtDev/Include/EGkSurfAux.h" #include "/EgtDev/Include/EGkSurfBezier.h" using namespace std ; //---------------------------------------------------------------------------- ISurf* NurbsToBezierSurface(const CNurbsSurfData& cnData) { //INTVECTOR vInt_sub( 10) ; //INTMATRIX vInt( 10, vInt_sub) ; //for ( int i = 0 ; i < 10 ; ++i ) { // for ( int j = 0 ; j < 10 ; ++j ) { // vInt[i][j] = i + 10 * j ; // } //} //vInt_sub.resize( 20) ; //vInt[0].resize( 20) ; // la superficie Nurbs deve essere in forma canonica if ( cnData.bPeriodicU || cnData.bPeriodicV || cnData.bExtraKnotes ) return nullptr ; // controllo sul numero dei nodi int nU = cnData.nCPU + cnData.nDegU - 1 ; int nV = cnData.nCPV + cnData.nDegV - 1 ; // controllo nodi e punti di controllo //if ( nU != int( cnData.vU.size()) || nV != int( cnData.vV.size()) || cnData.nCPU * cnData.nCPV != int( cnData.vCP.size())) if ( nU != int(cnData.vU.size()) || nV != int(cnData.vV.size())) { return nullptr ; } //// numero degli intervalli //int nInt = nU - 2 * cnData.nDeg + 1 ; // verifico le condizioni agli estremi sui nodi (i primi nDeg nodi e gli ultimi nDeg nodi devono essere uguali tra loro) bool bOk = true ; // direzione U for ( int i = 1 ; i < cnData.nDegU ; ++ i) { if ( abs( cnData.vU[i] - cnData.vU[0]) >= EPS_ZERO) bOk = false ; } for ( int i = 1 ; i < cnData.nDegU ; ++ i) { if ( abs( cnData.vU[nU - 1 - i] - cnData.vU[nU - 1]) >= EPS_ZERO) bOk = false ; } // direzione V for ( int i = 1 ; i < cnData.nDegV ; ++ i) { if ( abs( cnData.vV[i] - cnData.vV[0]) >= EPS_ZERO) bOk = false ; } for ( int i = 1 ; i < cnData.nDegV ; ++ i) { if ( abs( cnData.vV[nV - 1 - i] - cnData.vV[nV - 1]) >= EPS_ZERO) bOk = false ; } if ( ! bOk) return nullptr ; //// se 1 solo intervallo, la Nurbs è già una curva di Bezier //if ( nInt == 1) { // // creo la curva di Bezier // PtrOwner pCrvBez( CreateCurveBezier()) ; // if ( IsNull( pCrvBez)) // return nullptr ; // // la inizializzo // if ( ! pCrvBez->Init( cnData.nDeg, cnData.bRat)) // return nullptr ; // for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { // if ( ! cnData.bRat) { // if ( ! pCrvBez->SetControlPoint( i, cnData.vCP[i])) // return nullptr ; // } // else { // if ( ! pCrvBez->SetControlPoint( i, cnData.vCP[i], cnData.vW[i])) // return nullptr ; // } // } // // se non è una curva ma un punto, la invalido // if ( pCrvBez->IsAPoint()) // pCrvBez->Init( cnData.nDeg, cnData.bRat) ; // // restituisco la curva // return Release( pCrvBez) ; //} // algoritmo 5.7 del libro "The NURBS book"////////////////////////////////////////////////////////////////////////////////////////////////////////////// // creazione delle strips nella direzione U ( trasformo le curve iso con U costante in bezier) int a = cnData.nDegU - 1 ; int b = cnData.nDegU ; int nb = 0 ; // numero di strisce in U ( lunghezza con U costante) //PNTVECTOR vBC ; //vBC.resize( cnData.nCPV * cnData.nDegU) ; //for (int row = 0 ; row < cnData.nCPV ; ++row ) { // for ( int i = 0 ; i <= cnData.nDegU ; ++i ) { // vBC[nDegU*row + i] = ( cnData.vCP[nCPU*row + i]) ; // } //} vector vCPV( cnData.nCPV) ; vector< vector> mBC (cnData.nDegU + 1,vCPV ) ; vector< vector> mBC_next (cnData.nDegU - 1, vCPV) ; vector< vector> mPC_strip(cnData.nDegU + 1, vCPV) ; // matrice che verrà ingrandita e conterrà la superficie metà bezier e metà NURBS DBLVECTOR vV_W( cnData.nCPV) ; vector mW( cnData.nDegU + 1, vV_W) ; vector mW_next( cnData.nDegU - 1, vV_W) ; vector mW_strip( cnData.nDegU + 1, vV_W) ; DBLVECTOR vAlpha ; vAlpha.resize( cnData.nDegU - 1) ; if ( ! cnData.bRat ) { for ( int i = 0 ; i <= cnData.nDegU ; ++i ) { for ( int row = 0 ; row < cnData.nCPV ; ++row ) { mBC[i][row] = cnData.mCP[i][row] ; } } } else { for ( int i = 0 ; i <= cnData.nDegU ; ++i ) { for (int row = 0 ; row < cnData.nCPV ; ++ row) { mW[i][row] = cnData.mW[i][row] ; mBC[i][row] = cnData.mCP[i][row] * cnData.mW[i][row] ; } } } bool bRef = false ; while ( b < nU - 1) { // qui correggo un probabile errore, mettendo nU anziché nCPV, come indicato nell'algoritmo int i = b ; while ( b < nU - 1 && abs( cnData.vU[b+1] - cnData.vU[b]) < EPS_ZERO) ++ b ; int mult = b - i + 1 ; if ( mult < cnData.nDegU ) { bRef = true ; // calcolo numeratore e alpha double numer = cnData.vU[b] - cnData.vU[a] ; for ( int j = cnData.nDegU ; j > mult ; -- j) vAlpha[j-mult-1] = numer / ( cnData.vU[a+j] - cnData.vU[a]) ; int r = cnData.nDegU - mult ; for ( int j = 1 ; j <= cnData.nDegU - mult ; ++j ) { int save = r - j ; int s = mult + j ; //for ( int row = 0 ; row < cnData.nCPV ; ++row) { // for ( int k = cnData.nDegU ; k >= s ; --k ) { // vBC[nCPU*row + k] = vAlpha[k-s]*vBC[nCPU*row + k] + ( 1 - vAlfa[k-s]) * vBC[nCPU*row + k - 1] // } //} if ( ! cnData.bRat ) { for ( int k = cnData.nDegU ; k >= s ; --k ) { for ( int row = 0 ; row < cnData.nCPV ; ++row) { mBC[k][row] = vAlpha[k-s] * mBC[k][row] + ( 1 - vAlpha[k-s]) * mBC[k-1][row] ; } } } else { for ( int k = cnData.nDegU ; k >= s ; --k ) { for ( int row = 0 ; row < cnData.nCPV ; ++row) { mBC[k][row] = vAlpha[k-s] * mBC[k][row] + ( 1 - vAlpha[k-s]) * mBC[k-1][row] ; mW[k][row] = vAlpha[k-s] * mW[k][row] + ( 1 - vAlpha[k-s]) * mW[k-1][row] ; } } } if ( b < nU - 1 ) { for ( int row = 0 ; row < cnData.nCPV ; ++row) { mBC_next[save][row] = mBC[cnData.nDegU][row] ; } if ( cnData.bRat ) for ( int row = 0 ; row < cnData.nCPV ; ++row) { mW_next[save][row] = mW[cnData.nDegU][row] ; } } } mPC_strip.resize( cnData.nDegU * ( nb + 1) + 1 , vCPV) ; mW_strip.resize( cnData.nDegU * ( nb + 1) + 1, vV_W) ; if ( ! cnData.bRat) for ( int i = 0 ; i <= cnData.nDegU ; ++i) { for ( int row = 0 ; row < cnData.nCPV ; ++row ) { mPC_strip[i+ nb * cnData.nDegU][row] = mBC[i][row] ; } } else { for ( int i = 0 ; i <= cnData.nDegU ; ++i) { for ( int row = 0 ; row < cnData.nCPV ; ++row ) { mPC_strip[i+ nb * cnData.nDegU][row] = mBC[i][row]/mW[i][row] ; mW_strip[i+ nb * cnData.nDegU][row] = mW[i][row] ; } } } } ++ nb ; // ho finito di definire la patch di Bezier attuale e passo alla successiva // aggiorno mBC con i valori della prossima pezza di Bezier // corrisponde a nb = nb + 1 if ( ! cnData.bRat){ for (int i = 0 ; i < cnData.nDegU - 1 ; ++ i) { for ( int row = 0 ; row < cnData.nCPV ; ++row) { mBC[i][row] = mBC_next[i][row] ; } } } else { for (int i = 0 ; i < cnData.nDegU - 1 ; ++ i) { for ( int row = 0 ; row < cnData.nCPV ; ++row) { mBC[i][row] = mBC_next[i][row] ; mW[i][row] = mW_next[i][row] ; } } } if ( b < nU - 1 ) { for ( int i = cnData.nDegU - mult ; i <= cnData.nDegU ; ++ i) { for (int row = 0 ; row < cnData.nCPV ; ++ row ) { mBC[i][row] = cnData.mCP[b - cnData.nDegU + i + 1][row] ; } } if ( cnData.bRat ) { for ( int i = cnData.nDegU - mult ; i <= cnData.nDegU ; ++ i) { for (int row = 0 ; row < cnData.nCPV ; ++ row ) { mW[i][row] = cnData.mW[b - cnData.nDegU + i + 1][row] ; } } } a = b ; ++b ; } } // se non ho raffinato allora tutti i nodi avevano già molteplicità massima. Converto direttamente in Bezier la dir U int nCPU_ref ; // numero dei punti di controllo in U dopo il raffinamento if ( ! bRef ) { nCPU_ref = cnData.nCPU ; mPC_strip.resize( cnData.nCPU, vCPV) ; mW_strip.resize( cnData.nCPU, vV_W) ; if ( ! cnData.bRat) { for ( int i = 0 ; i < nCPU_ref ; ++i) { for ( int row = 0 ; row < cnData.nCPV ; ++ row) { mPC_strip[i][row] = cnData.mCP[i][row] ; } } } else { for ( int i = 0 ; i < nCPU_ref ; ++i) { for ( int row = 0 ; row < cnData.nCPV ; ++ row) { mPC_strip[i][row] = cnData.mCP[i][row] ; mW_strip[i][row] = cnData.mW[i][row] ; } } } // devo vedere quante patch ci stanno prendendo i punti che ci sono nb = (cnData.nCPU - 1) / cnData.nDegU ; } else nCPU_ref = cnData.nDegU * nb + 1 ; // numero dei punti di controllo in U dopo il raffinamento // ora ho ottenuto le strisce nDegU x nCPV // devo ripetere la procedura, sulla dir V, per ottenere le patch nDegU x nDegV a = cnData.nDegV - 1 ; b = cnData.nDegV ; int nc = 0 ; // numero di strisce in V ( lunghezza con V costante) vector vDegV(cnData.nDegV + 1) ; vector vDegV_1(cnData.nDegV - 1) ; vector< vector> m_BC1( nCPU_ref, vDegV) ; vector< vector> m_BC1_next( nCPU_ref, vDegV_1) ; DBLVECTOR vV1_W(cnData.nDegV + 1) ; DBLVECTOR vV2_W(cnData.nDegV - 1) ; vector mW1( nCPU_ref, vV1_W) ; vector mW1_next( nCPU_ref, vV2_W) ; DBLVECTOR vAlpha1( cnData.nDegV - 1) ; vector> mPC_tot( nCPU_ref, vDegV) ; vector mW_tot( nCPU_ref, vV1_W) ; if ( ! cnData.bRat ) { for ( int i = 0 ; i < nCPU_ref ; ++i ) { for ( int row = 0 ; row <= cnData.nDegV ; ++row ) { m_BC1[i][row] = mPC_strip[i][row] ; } } } else { for ( int i = 0 ; i < nCPU_ref ; ++i ) { for (int row = 0 ; row <= cnData.nDegV ; ++ row) { mW1[i][row] = mW_strip[i][row] ; m_BC1[i][row] = mPC_strip[i][row] * mW_strip[i][row] ; } } } bRef = false ; while ( b < nV - 1) { // qui correggo un probabile errore, mettendo nU anziché nCPV, come indicato nell'algoritmo int i = b ; while ( b < nV - 1 && abs( cnData.vV[b+1] - cnData.vV[b]) < EPS_ZERO) ++ b ; int mult = b - i + 1 ; if ( mult < cnData.nDegV ) { bRef = true ; // calcolo numeratore e alpha double numer = cnData.vV[b] - cnData.vV[a] ; for ( int j = cnData.nDegV ; j > mult ; -- j) vAlpha1[j-mult-1] = numer / ( cnData.vV[a+j] - cnData.vV[a]) ; int r = cnData.nDegV - mult ; for ( int j = 1 ; j <= cnData.nDegV - mult ; ++j ) { int save = r - j ; int s = mult + j ; //for ( int row = 0 ; row < cnData.nCPV ; ++row) { // for ( int k = cnData.nDegU ; k >= s ; --k ) { // vBC[nCPU*row + k] = vAlpha1[k-s]*vBC[nCPU*row + k] + ( 1 - vAlpha1[k-s]) * vBC[nCPU*row + k - 1] // } //} if ( ! cnData.bRat) { for ( int k = 0 ; k < nCPU_ref ; ++k) { for ( int row = cnData.nDegV ; row >= s ; --row ) { m_BC1[k][row] = vAlpha1[row-s] * m_BC1[k][row] + ( 1 - vAlpha1[row-s]) * m_BC1[k][row-1] ; } } } else { for ( int k = 0 ; k < nCPU_ref ; ++k) { for ( int row = cnData.nDegV ; row >= s ; --row ) { m_BC1[k][row] = vAlpha1[row-s] * m_BC1[k][row] + ( 1 - vAlpha1[row-s]) * m_BC1[k][row-1] ; mW1[k][row] = vAlpha1[row-s] * mW1[k][row] + ( 1 - vAlpha1[row-s]) * mW1[k][row-1] ; } } } if ( b < nV - 1 ) { if ( !cnData.bRat ){ for ( int i = 0 ; i < nCPU_ref ; ++i) { m_BC1_next[i][save] = m_BC1[i][cnData.nDegV] ; } } else { for ( int i = 0 ; i < nCPU_ref ; ++i) { m_BC1_next[i][save] = m_BC1[i][cnData.nDegV] ; mW1_next[save] = mW1[cnData.nDegV] ; } } } } } int nRef = cnData.nDegV * ( nc + 1) + 1 ; for ( int k = 0 ; k < nCPU_ref; ++k){ mPC_tot[k].resize( nRef) ; mW_tot[k].resize( nRef) ; } if ( ! cnData.bRat) for ( int i = 0 ; i < nCPU_ref ; ++i) { for ( int row = 0 ; row <= cnData.nDegV ; ++row ) { mPC_tot[i][row + nc * cnData.nDegV] = m_BC1[i][row] ; } } else { for ( int i = 0 ; i < nCPU_ref ; ++i) { for ( int row = 0 ; row <= cnData.nDegV ; ++row ) { mPC_tot[i][row + nc * cnData.nDegV] = m_BC1[i][row]/mW1[i][row] ; mW_tot[i][row + nc * cnData.nDegV] = mW1[i][row] ; } } } ++ nc ; // ho finito di definire la patch di Bezier attuale e passo alla successiva // aggiorno mBC con i valori della prossima pezza di Bezier // corrisponde a nc = nc + 1 if ( ! cnData.bRat){ for (int i = 0 ; i < nCPU_ref ; ++ i) { for ( int row = 0 ; row < cnData.nDegV - 1 ; ++row) { m_BC1[i][row] = m_BC1_next[i][row] ; } } } else { for (int i = 0 ; i < nCPU_ref ; ++ i) { for ( int row = 0 ; row < cnData.nDegV - 1 ; ++row) { m_BC1[i][row] = m_BC1_next[i][row] ; mW1[i][row] = mW1_next[i][row] ; } } } if ( b < nV - 1) { for (int i = 0 ; i < nCPU_ref ; ++ i ) { for ( int row = cnData.nDegV - mult ; row <= cnData.nDegV ; ++ row) { m_BC1[i][row] = cnData.mCP[i][b - cnData.nDegV + row + 1] ; } } if ( cnData.bRat ) { for (int i = 0 ; i < nCPU_ref ; ++ i ) { for ( int row = cnData.nDegV - mult ; row <= cnData.nDegV ; ++ row) { mW1[i][row] = cnData.mW[i][b - cnData.nDegV + row + 1] ; } } } a = b ; ++b ; } } // se non ho raffinato allora aggiungo direttamente alle matrici della superficie totale int nCPV_ref ; // numero dei punti di controllo in V dopo il raffinamento if ( ! bRef) { nCPV_ref = cnData.nCPV ; for ( int k = 0 ; k < nCPU_ref ; ++k){ mPC_tot[k].resize( cnData.nCPV) ; mW_tot[k].resize( cnData.nCPV) ; } if ( ! cnData.bRat) { for ( int i = 0 ; i < nCPU_ref ; ++i) { for ( int row = 0 ; row < nCPV_ref ; ++ row) { mPC_tot[i][row] = mPC_strip[i][row] ; } } } else { for ( int i = 0 ; i < nCPU_ref ; ++i) { for ( int row = 0 ; row < nCPV_ref ; ++ row) { mPC_tot[i][row] = mPC_strip[i][row] ; mW_tot[i][row] = mW_strip[i][row] ; } } } // devo vedere quante patch ci stanno prendendo i punti che ci sono nc = (cnData.nCPV - 1) / cnData.nDegV ; } else nCPV_ref = cnData.nDegV * nc + 1 ; // finalmente setto la superficie di bezier totale divisa in nb patch in U e nc patch in V PtrOwner pSrfBz( CreateSurfBezier()) ; if ( IsNull( pSrfBz)) return nullptr ; pSrfBz->Init(cnData.nDegU, cnData.nDegV, nb, nc, cnData.bRat) ; if ( !cnData.bRat ) { for ( int i = 0 ; i < nCPU_ref; ++ i) { for (int j = 0 ; j < nCPV_ref; ++j) { pSrfBz->SetControlPoint( i + nCPU_ref * j, mPC_tot[i][j]) ; } } } else { for ( int i = 0 ; i < nCPU_ref; ++ i) { for (int j = 0 ; j < nCPV_ref; ++j) { pSrfBz->SetControlPoint( i + nCPU_ref * j, mPC_tot[i][j] / mW_tot[i][j], mW_tot[i][j]) ; } } } return Release( pSrfBz) ; } //// algoritmo per le curve, da usare come reference/////////////////////////////////////////////////////////////////////////////// // // // // vettore dei punti di controllo della curva di Bezier // PNTVECTOR vBC ; // vBC.resize( cnData.nDeg + 1) ; // DBLVECTOR vBW ; // vBW.resize( cnData.nDeg + 1) ; // if ( ! cnData.bRat) { // for ( int i = 0 ; i <= cnData.nDeg ; ++ i) // vBC[i] = cnData.vCP[i] ; // } // else { // for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { // vBC[i] = cnData.vCP[i] * cnData.vW[i] ; // vBW[i] = cnData.vW[i] ; // } // } // // primi coefficienti della successiva // PNTVECTOR vNextBC ; // vNextBC.resize( cnData.nDeg - 1) ; // DBLVECTOR vNextBW ; // vNextBW.resize( cnData.nDeg - 1) ; // // ... // DBLVECTOR vAlfa ; // vAlfa.resize( cnData.nDeg - 1) ; // int a = cnData.nDeg - 1 ; // int b = cnData.nDeg ; // bool bPrevRejected = false ; // // ciclo // while ( b < nU - 1) { // int i = b ; // while ( b < nU - 1 && abs( cnData.vU[b+1] - cnData.vU[b]) < EPS_ZERO) // ++ b ; // int mult = min( b - i + 1, cnData.nDeg) ; // if ( mult < cnData.nDeg) { // // numeratore di alfa // double numer = cnData.vU[b] - cnData.vU[a] ; // // calcola e salva gli alfa // for ( int j = cnData.nDeg ; j > mult ; -- j) // vAlfa[j-mult-1] = numer / ( cnData.vU[a+j] - cnData.vU[a]) ; // // inserisco il nodo r volte // int r = cnData.nDeg - mult ; // for ( int j = 1 ; j <= r ; ++ j) { // int save = r - j ; // int s = mult + j ; // for ( int k = cnData.nDeg ; k >= s ; -- k) // vBC[k] = vAlfa[k-s] * vBC[k] + ( 1 - vAlfa[k-s]) * vBC[k-1] ; // if ( cnData.bRat) { // for ( int k = cnData.nDeg ; k >= s ; -- k) // vBW[k] = vAlfa[k-s] * vBW[k] + ( 1 - vAlfa[k-s]) * vBW[k-1] ; // } // if ( b < nU - 1) { // vNextBC[save] = vBC[cnData.nDeg] ; // if ( cnData.bRat) // vNextBW[save] = vBW[cnData.nDeg] ; // } // } // } // // // costruisco la curva di Bezier e la inserisco nella curva composita // PtrOwner pCrvBez( CreateCurveBezier()) ; // if ( IsNull( pCrvBez)) // return nullptr ; // // se precedente saltata // if ( bPrevRejected) { // // prendo l'ultimo punto della curva composita per garantire la continuità // Point3d ptEnd ; // if ( pCrvCompo->GetEndPoint( ptEnd)) // vBC[0] = ptEnd ; // } // // la inizializzo // if ( ! pCrvBez->Init( cnData.nDeg, cnData.bRat)) // return nullptr ; // if ( ! cnData.bRat) { // for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { // if ( ! pCrvBez->SetControlPoint( i, vBC[i])) // return nullptr ; // } // } // else { // for ( int i = 0 ; i <= cnData.nDeg ; ++ i) { // if ( ! pCrvBez->SetControlPoint( i, vBC[i] / vBW[i], vBW[i])) // return nullptr ; // } // } // // se è una vera curva, la aggiungo alla curva composita // if ( ! pCrvBez->IsAPoint()) { // if ( ! pCrvCompo->AddCurve( Release( pCrvBez))) // return nullptr ; // bPrevRejected = false ; // } // // altrimenti è un punto, la cancello // else { // pCrvBez.Reset() ; // bPrevRejected = true ; // } // // // inizializzazioni per la prossima curva di Bezier // if ( b < nU - 1) { // if ( ! cnData.bRat) { // for ( int i = 0 ; i < cnData.nDeg - 1 ; ++ i) // vBC[i] = vNextBC[i] ; // for ( int i = cnData.nDeg - mult ; i <= cnData.nDeg ; ++ i) // vBC[i] = cnData.vCP[b-cnData.nDeg+i+1] ; // } // else { // for ( int i = 0 ; i < cnData.nDeg - 1 ; ++ i) { // vBC[i] = vNextBC[i] ; // vBW[i] = vNextBW[i] ; // } // for ( int i = cnData.nDeg - mult ; i <= cnData.nDeg ; ++ i) { // vBC[i] = cnData.vCP[b-cnData.nDeg+i+1] * cnData.vW[b-cnData.nDeg+i+1] ; // vBW[i] = cnData.vW[b-cnData.nDeg+i+1] ; // } // } // a = b ; // ++ b ; // } // } // // // restituisco la curva composita // return Release( pCrvCompo) ; //}