diff --git a/VolZmapGraphics.cpp b/VolZmapGraphics.cpp index 7f66c45..5dc1c13 100644 --- a/VolZmapGraphics.cpp +++ b/VolZmapGraphics.cpp @@ -292,6 +292,21 @@ CreateBigTriangleYZ( double dMinY, double dMaxY, double dMinZ, double dMaxZ, dou return trTria1.Validate( true) && trTria2.Validate( true) ; } +//---------------------------------------------------------------------------- +double +GetFanTotalNormalDeviation( Point3d ptP, AppliedVector Vector[7], int nVertNum) +{ + double dSum = 0 ; + for ( int m = 0 ; m < nVertNum ; ++ m) { + int l = ( m + 1) % nVertNum ; + Vector3d vtNlm = ( Vector[l].ptPApp - ptP) ^ ( Vector[m].ptPApp - ptP) ; + vtNlm.Normalize() ; + dSum += vtNlm * Vector[l].vtVec ; + dSum += vtNlm * Vector[m].vtVec ; + } + return dSum ; +} + //---------------------------------------------------------------------------- bool VolZmap::FindAdjComp( const vector& vVecVox, int nCurBlock, int nCurVox, int nCurComp, @@ -1998,6 +2013,20 @@ VolZmap::ExtMarchingCubes( int nBlock, VoxelContainer& vVox) const bool bExtConfirmed = true ; Vector3d vtNullSpace( dMatrixV( 0, 2), dMatrixV( 1, 2), dMatrixV( 2, 2)) ; if ( nFeatureType == EDGE && vtNullSpace.Normalize() && ! IsPointInsideVoxelApprox( i, j, k, ptSol)) { + double dBestDisp = 0 ; + double dMaxVal = - DBL_MAX ; + double dDisplH = 0.1 * m_dStep ; + double dDispStop = 0.5 * m_dStep ; + double dDisp = - 0.5 * m_dStep; + while ( dDisp < dDispStop + EPS_SMALL) { + double dCrVal = GetFanTotalNormalDeviation( ptSol + dDisp * vtNullSpace, CompoVert[nComp], nVertComp[nComp]) ; + if ( dCrVal > dMaxVal) { + dMaxVal = dCrVal ; + dBestDisp = dDisp ; + } + dDisp += dDisplH ; + } + ptSol += dBestDisp * vtNullSpace ; int nNewAdjVoxIJK[3] ; if ( GetPointVoxel( ptSol, nNewAdjVoxIJK[0], nNewAdjVoxIJK[1], nNewAdjVoxIJK[2])) { int nNewAdjVoxIndex = CalcIndex( nNewAdjVoxIJK[0], nNewAdjVoxIJK[1], nNewAdjVoxIJK[2]) ;