Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 11 additions & 3 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,22 @@ PR_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cp
SR_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cpp SSDRecon.cpp
ST_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cpp SurfaceTrimmer.cpp

CFLAGS += -fopenmp -Wno-deprecated -Wno-write-strings -std=c++11
LFLAGS += -lgomp -lstdc++
CFLAGS += -Wno-deprecated -Wno-write-strings -std=c++11
LFLAGS += -lstdc++

#
# OpenMP
#
ifndef NO_OPENMP
CFLAGS += -fopenmp
LFLAGS += -lgomp
endif

CFLAGS_DEBUG = -DDEBUG -g3
LFLAGS_DEBUG =

CFLAGS_RELEASE = -O3 -DRELEASE -funroll-loops -ffast-math
LFLAGS_RELEASE = -O3
LFLAGS_RELEASE = -O3

SRC = Src/
BIN = Bin/Linux/
Expand Down
8 changes: 5 additions & 3 deletions Src/MAT.inl
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,14 @@ are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer. Redistributions in binary form must reproduce
the above copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the distribution.
in the documentation and/or other materials provided with the distribution.

Neither the name of the Johns Hopkins University nor the names of its contributors
may be used to endorse or promote products derived from this software without specific
prior written permission.
prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
Expand Down Expand Up @@ -146,7 +146,9 @@ void MinimalAreaTriangulation<Real>::GetTriangulation(const size_t& i,const size
if( j+1>=ii )
return;
ii=midPoint[i*eCount+j];
#ifdef BRUNO_LEVY_FIX // Always true when type of ii is size_t
if( ii>=0 )
#endif
{
tIndex.idx[0] = int( i );
tIndex.idx[1] = int( j );
Expand Down
22 changes: 11 additions & 11 deletions Src/MarchingCubes.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,14 @@ are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer. Redistributions in binary form must reproduce
the above copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the distribution.
in the documentation and/or other materials provided with the distribution.

Neither the name of the Johns Hopkins University nor the names of its contributors
may be used to endorse or promote products derived from this software without specific
prior written permission.
prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
Expand Down Expand Up @@ -145,9 +145,9 @@ bool Cube::IsEdgeCorner( int cIndex , int e )
FactorEdgeIndex( e , o , i , j );
switch( o )
{
case 0: return (cIndex && 2)==(i<<1) && (cIndex && 4)==(j<<2);
case 1: return (cIndex && 1)==(i<<0) && (cIndex && 4)==(j<<2);
case 2: return (cIndex && 4)==(i<<2) && (cIndex && 2)==(j<<1);
case 0: return (cIndex & 2)==(i<<1) && (cIndex & 4)==(j<<2);
case 1: return (cIndex & 1)==(i<<0) && (cIndex & 4)==(j<<2);
case 2: return (cIndex & 4)==(i<<2) && (cIndex & 2)==(j<<1);
default: return false;
}
}
Expand Down Expand Up @@ -309,7 +309,7 @@ bool MarchingSquares::HasEdgeRoots( unsigned char mcIndex , int edgeIndex )
{
int c1 , c2;
Square::EdgeCorners( edgeIndex , c1 , c2 );
return !(
return !(
( ( mcIndex&(1<<MarchingSquares::cornerMap[c1]) ) && ( mcIndex&(1<<MarchingSquares::cornerMap[c2])) )
||
(!( mcIndex&(1<<MarchingSquares::cornerMap[c1]) ) && !( mcIndex&(1<<MarchingSquares::cornerMap[c2])) )
Expand All @@ -325,21 +325,21 @@ bool MarchingSquares::HasEdgeRoots( unsigned char mcIndex , int edgeIndex )
const int MarchingSquares::edgeMask[1<<Square::CORNERS]=
{
0, // 0 -> -> ->
9, // 1 -> 0 -> (0,0) -> 0,3 -> 9
9, // 1 -> 0 -> (0,0) -> 0,3 -> 9
3, // 2 -> 1 -> (1,0) -> 0,1 -> 3
10, // 3 -> 0,1 -> (0,0) (1,0) -> 1,3 -> 10
12, // 4 -> 2 -> (0,1) -> 2,3 -> 12
5, // 5 -> 0,2 -> (0,0) (0,1) -> 0,2 -> 5
15, // 6 -> 1,2 -> (1,0) (0,1) -> 0,1,2,3 -> 15
6, // 7 -> 0,1,2 -> (0,0) (1,0) (0,1) -> 1,2 -> 6
6, // 8 -> 3 -> (1,1) -> 1,2 -> 6
15, // 9 -> 0,3 -> (0,0) (1,1) -> 0,1,2,3 -> 15
15, // 9 -> 0,3 -> (0,0) (1,1) -> 0,1,2,3 -> 15
5, // 10 -> 1,3 -> (1,0) (1,1) -> 0,2 -> 5
12, // 11 -> 0,1,3 -> (0,0) (1,0) (1,1) -> 2,3 -> 12
10, // 12 -> 2,3 -> (0,1) (1,1) -> 1,3 -> 10
3, // 13 -> 0,2,3 -> (0,0) (0,1) (1,1) -> 0,1 -> 3
9, // 14 -> 1,2,3 -> (1,0) (0,1) (1,1) -> 0,3 -> 9
0, // 15 -> 0,1,2,3 -> (0,0) (1,0) (0,1) (1,1) ->
0, // 15 -> 0,1,2,3 -> (0,0) (1,0) (0,1) (1,1) ->
};
#if NEW_ORDERING
/*
Expand Down Expand Up @@ -950,7 +950,7 @@ bool MarchingCubes::HasEdgeRoots( unsigned char mcIndex , int edgeIndex )
{
int c1 , c2;
Cube::EdgeCorners( edgeIndex , c1 , c2 );
return !(
return !(
( ( mcIndex&(1<<MarchingCubes::cornerMap[c1]) ) && ( mcIndex&(1<<MarchingCubes::cornerMap[c2])) )
||
(!( mcIndex&(1<<MarchingCubes::cornerMap[c1]) ) && !( mcIndex&(1<<MarchingCubes::cornerMap[c2])) )
Expand Down
59 changes: 32 additions & 27 deletions Src/MultiGridOctreeData.System.inl
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,14 @@ are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer. Redistributions in binary form must reproduce
the above copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the distribution.
in the documentation and/or other materials provided with the distribution.

Neither the name of the Johns Hopkins University nor the names of its contributors
may be used to endorse or promote products derived from this software without specific
prior written permission.
prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
Expand Down Expand Up @@ -45,7 +45,7 @@ struct _ConstraintCalculator_< Real , Degree , false >
{
#if POINT_DATA_RES
Real constraint = 0;
for( int c=0 ; c<PointData< Real , false >::SAMPLES ; c++ ) if( p[c].weight )
for( int c=0 ; c<PointData< Real , false >::SAMPLES ; c++ ) if( p[c].weight )
{
const Point3D< Real > q = p[c].position;
constraint += (Real)( px( q[0] ) * py( q[1] ) * pz( q[2] ) * p[c].weight * p[c].value );
Expand All @@ -60,7 +60,7 @@ struct _ConstraintCalculator_< Real , Degree , false >
{
#if POINT_DATA_RES
Real constraint = 0;
for( int c=0 ; c<PointData< Real , false >::SAMPLES ; c++ ) if( p[c].weight )
for( int c=0 ; c<PointData< Real , false >::SAMPLES ; c++ ) if( p[c].weight )
{
const Point3D< Real > q = p[c].position;
constraint += (Real)( px( q[0] ) * py( q[1] ) * pz( q[2] ) * p[c]._value );
Expand All @@ -84,7 +84,7 @@ struct _ConstraintCalculator_< Real , Degree , true >
{
#if POINT_DATA_RES
Real constraint = 0;
for( int c=0 ; c<PointData< Real , true >::SAMPLES ; c++ ) if( p[c].weight )
for( int c=0 ; c<PointData< Real , true >::SAMPLES ; c++ ) if( p[c].weight )
{
const Point3D< Real > q = p[c].position;
double _px = px( q[0] ) , _py = py( q[1] ) , _pz = pz( q[2] );
Expand All @@ -109,7 +109,7 @@ struct _ConstraintCalculator_< Real , Degree , true >
{
#if POINT_DATA_RES
Real constraint = 0;
for( int c=0 ; c<PointData< Real , true >::SAMPLES ; c++ ) if( p[c].weight )
for( int c=0 ; c<PointData< Real , true >::SAMPLES ; c++ ) if( p[c].weight )
{
const Point3D< Real > q = p[c].position;
double _px = px( q[0] ) , _py = py( q[1] ) , _pz = pz( q[2] );
Expand Down Expand Up @@ -178,7 +178,7 @@ double FEMSystemFunctor< 1 , BOUNDARY_FREE >::_integrate( const I& integrator ,
return
(
d00[0] * d00[1] * d00[2]
) * massWeight
) * massWeight
+
(
d11[0] * d00[1] * d00[2] +
Expand All @@ -196,7 +196,7 @@ double FEMSystemFunctor< 1 , BOUNDARY_NEUMANN >::_integrate( const I& integrator
return
(
d00[0] * d00[1] * d00[2]
) * massWeight
) * massWeight
+
(
d11[0] * d00[1] * d00[2] +
Expand All @@ -214,7 +214,7 @@ double FEMSystemFunctor< 1 , BOUNDARY_DIRICHLET >::_integrate( const I& integrat
return
(
d00[0] * d00[1] * d00[2]
) * massWeight
) * massWeight
+
(
d11[0] * d00[1] * d00[2] +
Expand All @@ -233,7 +233,7 @@ double FEMSystemFunctor< FEMDegree , BType >::_integrate( const I& integrator ,
return
(
d00[0] * d00[1] * d00[2]
) * massWeight
) * massWeight
+
(
d11[0] * d00[1] * d00[2] +
Expand All @@ -259,11 +259,11 @@ double FEMSFConstraintFunctor< SFDegree , SFBType , FEMDegree , FEMBType >::_int
double d00[] = D_DOT( 0 , 0 ) , d02[] = D_DOT( 0 , 2 ) , d20[] = D_DOT( 2 , 0 ) , d22[] = D_DOT( 2 , 2 ) , d11[] = D_DOT( 1 , 1 );
if( SFDegree==0 || FEMDegree==0 )
return d00[0] * d00[1] * d00[2] * massWeight;
else if( SFDegree<=1 || FEMDegree<=1 )
else if( SFDegree<=1 || FEMDegree<=1 )
return
(
d00[0] * d00[1] * d00[2]
) * massWeight
) * massWeight
+
(
d11[0] * d00[1] * d00[2] +
Expand All @@ -274,7 +274,7 @@ double FEMSFConstraintFunctor< SFDegree , SFBType , FEMDegree , FEMBType >::_int
return
(
d00[0] * d00[1] * d00[2]
) * massWeight
) * massWeight
+
(
d11[0] * d00[1] * d00[2] +
Expand Down Expand Up @@ -322,8 +322,8 @@ Point3D< double > FEMVFConstraintFunctor< VFDegree , VFBType , FEMDegree , FEMBT
+
Point3D< double >
(
d12[0] * d00[1] * d00[2] + d10[0] * ( d00[1] * d02[2] + d02[1] * d00[2] ) ,
d12[1] * d00[2] * d00[0] + d10[1] * ( d00[2] * d02[0] + d02[2] * d00[0] ) ,
d12[0] * d00[1] * d00[2] + d10[0] * ( d00[1] * d02[2] + d02[1] * d00[2] ) ,
d12[1] * d00[2] * d00[0] + d10[1] * ( d00[2] * d02[0] + d02[2] * d00[0] ) ,
d12[2] * d00[0] * d00[1] + d10[2] * ( d00[0] * d02[1] + d02[0] * d00[1] )
) * biLapWeight;
}
Expand Down Expand Up @@ -539,7 +539,7 @@ void Octree< Real >::_upSample( LocalDepth highDepth , DenseNodeData< C , FEMDeg
BSplineEvaluationData< FEMDegree , BType >::SetUpSampleEvaluator( upSampleEvaluator , lowDepth );
std::vector< DownSampleKey > neighborKeys( std::max< int >( 1 , threads ) );
for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( _localToGlobal( lowDepth ) );

static const int DownSampleSize = BSplineSupportSizes< FEMDegree >::DownSample0Size > BSplineSupportSizes< FEMDegree >::DownSample1Size ? BSplineSupportSizes< FEMDegree >::DownSample0Size : BSplineSupportSizes< FEMDegree >::DownSample1Size;
Stencil< double , DownSampleSize > downSampleStencils[ Cube::CORNERS ];
int lowCenter = ( 1<<lowDepth )>>1;
Expand All @@ -550,7 +550,7 @@ void Octree< Real >::_upSample( LocalDepth highDepth , DenseNodeData< C , FEMDeg
for( int ii=0 ; ii<BSplineSupportSizes< FEMDegree >::DownSampleSize[cx] ; ii++ )
for( int jj=0 ; jj<BSplineSupportSizes< FEMDegree >::DownSampleSize[cy] ; jj++ )
for( int kk=0 ; kk<BSplineSupportSizes< FEMDegree >::DownSampleSize[cz] ; kk++ )
downSampleStencils[c]( ii , jj , kk ) =
downSampleStencils[c]( ii , jj , kk ) =
upSampleEvaluator.value( lowCenter + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] , 2*lowCenter + cx ) *
upSampleEvaluator.value( lowCenter + jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] , 2*lowCenter + cy ) *
upSampleEvaluator.value( lowCenter + kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] , 2*lowCenter + cz ) ;
Expand Down Expand Up @@ -640,7 +640,7 @@ void Octree< Real >::_UpSample( LocalDepth highDepth , ConstPointer( C ) lowCoef
for( int ii=0 ; ii<BSplineSupportSizes< FEMDegree >::DownSampleSize[cx] ; ii++ )
for( int jj=0 ; jj<BSplineSupportSizes< FEMDegree >::DownSampleSize[cy] ; jj++ )
for( int kk=0 ; kk<BSplineSupportSizes< FEMDegree >::DownSampleSize[cz] ; kk++ )
downSampleStencils[c]( ii , jj , kk ) =
downSampleStencils[c]( ii , jj , kk ) =
upSampleEvaluator.value( lowCenter + ii + BSplineSupportSizes< FEMDegree >::DownSampleStart[cx] , 2*lowCenter + cx ) *
upSampleEvaluator.value( lowCenter + jj + BSplineSupportSizes< FEMDegree >::DownSampleStart[cy] , 2*lowCenter + cy ) *
upSampleEvaluator.value( lowCenter + kk + BSplineSupportSizes< FEMDegree >::DownSampleStart[cz] , 2*lowCenter + cz ) ;
Expand Down Expand Up @@ -855,7 +855,7 @@ Real Octree< Real >::_finerFunctionValue( Point3D< Real > p , const PointSupport
{
int fIdx[3];
functionIndex< FEMDegree , BType >( _node , fIdx );
pointValue +=
pointValue +=
bsData.baseBSplines[ fIdx[0] ][LeftSupportRadius-j]( p[0] ) *
bsData.baseBSplines[ fIdx[1] ][LeftSupportRadius-k]( p[1] ) *
bsData.baseBSplines[ fIdx[2] ][LeftSupportRadius-l]( p[2] ) *
Expand Down Expand Up @@ -918,15 +918,15 @@ void Octree< Real >::_setPointValuesFromCoarser( InterpolationInfo< HasGradients
c , *pData ,
_coarserFunctionValue( (*pData)[c].position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) ,
HasGradients ? _coarserFunctionGradient( (*pData)[c].position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) : Point3D< Real >() ,
interpolationInfo.valueWeight , interpolationInfo.gradientWeight
interpolationInfo.valueWeight , interpolationInfo.gradientWeight
);
#else // !POINT_DATA_RES
_ConstraintCalculator_< Real , FEMDegree , HasGradients >::_CalculateCoarser_
(
*pData ,
_coarserFunctionValue( pData->position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) ,
HasGradients ? _coarserFunctionGradient( pData->position , neighborKey , _sNodes.treeNodes[i] , bsData , upSampledCoefficients ) : Point3D< Real >() ,
interpolationInfo.valueWeight , interpolationInfo.gradientWeight
interpolationInfo.valueWeight , interpolationInfo.gradientWeight
);
#endif // POINT_DATA_RES
}
Expand Down Expand Up @@ -1118,7 +1118,7 @@ int Octree< Real >::_setMatrixRow( const FEMSystemFunctor& F , const Interpolati
if( interpolationInfo )
{
memset( pointValues , 0 , sizeof(Real) * OverlapSize * OverlapSize * OverlapSize );
// Iterate over all supported neighbors that could have a point constraint
// Iterate over all supported neighbors that could have a point constraint
for( int i=-LeftSupportRadius ; i<=RightSupportRadius ; i++ ) if( hasYZPoints[i+LeftSupportRadius] )
for( int j=-LeftSupportRadius ; j<=RightSupportRadius ; j++ ) if( hasZPoints[i+LeftSupportRadius][j+LeftSupportRadius] )
for( int k=-LeftSupportRadius ; k<=RightSupportRadius ; k++ )
Expand Down Expand Up @@ -1805,7 +1805,7 @@ void Octree< Real >::_addFEMConstraints( const FEMConstraintFunctor& F , const C
// splatted normals and compute the dot-product of the
// divergence of the normal field with all the basis functions.
// Within the same depth: set directly as a gather
// Coarser depths
// Coarser depths
maxDepth = std::min< LocalDepth >( maxDepth , _maxDepth );
DenseNodeData< Real , FEMDegree >* __constraints = new DenseNodeData< Real , FEMDegree >( _sNodesEnd(maxDepth-1) );
DenseNodeData< Real , FEMDegree >& _constraints = *__constraints;
Expand Down Expand Up @@ -1847,14 +1847,17 @@ void Octree< Real >::_addFEMConstraints( const FEMConstraintFunctor& F , const C
if( isValidFEMNode< CDegree , CBType >( _node ) )
{
const D* d = coefficients( _node );
if( d )
if( isInterior ) constraints[i] += _Dot( (D)stencil( x , y , z ) , *d );
if( d )
{
if( isInterior )
constraints[i] += _Dot( (D)stencil( x , y , z ) , *d );
else
{
LocalDepth _d ; LocalOffset _off;
_localDepthAndOffset( _node , _d , _off );
constraints[i] += _Dot( *d , (D)F.template integrate< false >( integrator , _off , off ) );
}
}
}
}
_SetParentOverlapBounds< CDegree , FEMDegree >( node , startX , endX , startY , endY , startZ , endZ );
Expand Down Expand Up @@ -2071,13 +2074,15 @@ double Octree< Real >::_dot( const DotFunctor& F , const InterpolationInfo< HasG
const TreeOctNode* _node = neighbors.neighbors[x][y][z];
const Real* _data2;
if( isValidFEMNode< FEMDegree2 , FEMBType2 >( _node ) && ( _data2=coefficients2( _node ) ) )
{
if( isInterior ) dot += (*_data1) * (*_data2 ) * stencil( x , y , z );
else
{
LocalDepth _d ; LocalOffset _off;
_localDepthAndOffset( _node , _d , _off );
dot += (*_data1) * (*_data2) * F.template integrate< false >( integrator , off , _off );
}
}
}
}
}
Expand Down Expand Up @@ -2243,7 +2248,7 @@ double Octree< Real >::_dot( const DotFunctor& F , const InterpolationInfo< HasG

const PointData< Real , HasGradients >& pData = *( (*iInfo)( _sNodes.treeNodes[i] ) );
#if POINT_DATA_RES
for( int c=0 ; c<PointData< Real , false >::SAMPLES ; c++ ) if( pData[c].weight )
for( int c=0 ; c<PointData< Real , false >::SAMPLES ; c++ ) if( pData[c].weight )
{
Point3D< Real > p = pData[c].position;
Real w = pData[c].weight;
Expand Down
Loading