From a7c1787fb9eb447e8d413b0bc5c2ec70f02096c5 Mon Sep 17 00:00:00 2001 From: Paolo Cignoni cignoni Date: Tue, 21 Oct 2014 16:57:55 +0000 Subject: [PATCH] Original version of the Misha code (version 6.11) --- .../filter_screened_poisson/Src/Geometry.cpp | 114 + .../filter_screened_poisson/Src/MAT.inl | 213 ++ .../filter_screened_poisson/Src/MemoryUsage.h | 200 ++ .../Src/MultiGridOctreeData.IsoSurface.inl | 1113 +++++++ .../Src/MultiGridOctreeData.h | 438 +++ .../filter_screened_poisson/Src/PPolynomial.h | 113 + .../Src/PPolynomial.inl | 431 +++ .../filter_screened_poisson/Src/Ply.h | 705 +++++ .../filter_screened_poisson/Src/PlyFile.cpp | 2727 +++++++++++++++++ .../filter_screened_poisson/Src/PlyFile.h | 221 ++ .../Src/PointStream.inl | 183 ++ .../Src/SparseMatrix.h | 218 ++ .../Src/SurfaceTrimmer.cpp | 426 +++ .../filter_screened_poisson/Src/Vector.inl | 302 ++ 14 files changed, 7404 insertions(+) create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/Geometry.cpp create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/MAT.inl create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/MemoryUsage.h create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.IsoSurface.inl create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.h create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.h create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.inl create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/Ply.h create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/PlyFile.cpp create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/PlyFile.h create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/PointStream.inl create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/SparseMatrix.h create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/SurfaceTrimmer.cpp create mode 100755 src/plugins_experimental/filter_screened_poisson/Src/Vector.inl diff --git a/src/plugins_experimental/filter_screened_poisson/Src/Geometry.cpp b/src/plugins_experimental/filter_screened_poisson/Src/Geometry.cpp new file mode 100755 index 000000000..afd03f107 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/Geometry.cpp @@ -0,0 +1,114 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ +#include "Geometry.h" +#include +#include + +/////////////////// +// CoredMeshData // +/////////////////// + +TriangulationEdge::TriangulationEdge(void){pIndex[0]=pIndex[1]=tIndex[0]=tIndex[1]=-1;} +TriangulationTriangle::TriangulationTriangle(void){eIndex[0]=eIndex[1]=eIndex[2]=-1;} + +/////////////////////////// +// BufferedReadWriteFile // +/////////////////////////// +BufferedReadWriteFile::BufferedReadWriteFile( char* fileName , int bufferSize ) +{ + _bufferIndex = 0; + _bufferSize = bufferSize; + if( fileName ) strcpy( _fileName , fileName ) , tempFile = false; +#ifdef _WIN32 + else strcpy( _fileName , _tempnam( "." , "foo" ) ) , tempFile = true; +#else // !_WIN32 + else strcpy( _fileName , tempnam( "." , "foo" ) ) , tempFile = true; +#endif // _WIN32 + _fp = fopen( _fileName , "w+b" ); + if( !_fp ) fprintf( stderr , "[ERROR] Failed to open file: %s\n" , _fileName ) , exit( 0 ); + _buffer = (char*) malloc( _bufferSize ); +} +BufferedReadWriteFile::~BufferedReadWriteFile( void ) +{ + free( _buffer ); + fclose( _fp ); + if( tempFile ) remove( _fileName ); +} +void BufferedReadWriteFile::reset( void ) +{ + if( _bufferIndex ) fwrite( _buffer , 1 , _bufferIndex , _fp ); + _bufferIndex = 0; + fseek( _fp , 0 , SEEK_SET ); + _bufferIndex = 0; + _bufferSize = fread( _buffer , 1 , _bufferSize , _fp ); +} +bool BufferedReadWriteFile::write( const void* data , size_t size ) +{ + if( !size ) return true; + char* _data = (char*) data; + size_t sz = _bufferSize - _bufferIndex; + while( sz<=size ) + { + memcpy( _buffer+_bufferIndex , _data , sz ); + fwrite( _buffer , 1 , _bufferSize , _fp ); + _data += sz; + size -= sz; + _bufferIndex = 0; + sz = _bufferSize; + } + if( size ) + { + memcpy( _buffer+_bufferIndex , _data , size ); + _bufferIndex += size; + } + return true; +} +bool BufferedReadWriteFile::read( void* data , size_t size ) +{ + if( !size ) return true; + char *_data = (char*) data; + size_t sz = _bufferSize - _bufferIndex; + while( sz<=size ) + { + if( size && !_bufferSize ) return false; + memcpy( _data , _buffer+_bufferIndex , sz ); + _bufferSize = fread( _buffer , 1 , _bufferSize , _fp ); + _data += sz; + size -= sz; + _bufferIndex = 0; + if( !size ) return true; + sz = _bufferSize; + } + if( size ) + { + if( !_bufferSize ) return false; + memcpy( _data , _buffer+_bufferIndex , size ); + _bufferIndex += size; + } + return true; +} \ No newline at end of file diff --git a/src/plugins_experimental/filter_screened_poisson/Src/MAT.inl b/src/plugins_experimental/filter_screened_poisson/Src/MAT.inl new file mode 100755 index 000000000..d6f2126e3 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/MAT.inl @@ -0,0 +1,213 @@ +/* +Copyright (c) 2007, Michael Kazhdan +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ +////////////////////////////// +// MinimalAreaTriangulation // +////////////////////////////// +template +MinimalAreaTriangulation::MinimalAreaTriangulation(void) +{ + bestTriangulation=NULL; + midPoint=NULL; +} +template +MinimalAreaTriangulation::~MinimalAreaTriangulation(void) +{ + if(bestTriangulation) + delete[] bestTriangulation; + bestTriangulation=NULL; + if(midPoint) + delete[] midPoint; + midPoint=NULL; +} +template +void MinimalAreaTriangulation::GetTriangulation(const std::vector >& vertices,std::vector& triangles) +{ + if(vertices.size()==3) + { + triangles.resize(1); + triangles[0].idx[0]=0; + triangles[0].idx[1]=1; + triangles[0].idx[2]=2; + return; + } + else if(vertices.size()==4) + { + TriangleIndex tIndex[2][2]; + Real area[2]; + + area[0]=area[1]=0; + triangles.resize(2); + + tIndex[0][0].idx[0]=0; + tIndex[0][0].idx[1]=1; + tIndex[0][0].idx[2]=2; + tIndex[0][1].idx[0]=2; + tIndex[0][1].idx[1]=3; + tIndex[0][1].idx[2]=0; + + tIndex[1][0].idx[0]=0; + tIndex[1][0].idx[1]=1; + tIndex[1][0].idx[2]=3; + tIndex[1][1].idx[0]=3; + tIndex[1][1].idx[1]=1; + tIndex[1][1].idx[2]=2; + + Point3D n,p1,p2; + for(int i=0;i<2;i++) + for(int j=0;j<2;j++) + { + p1=vertices[tIndex[i][j].idx[1]]-vertices[tIndex[i][j].idx[0]]; + p2=vertices[tIndex[i][j].idx[2]]-vertices[tIndex[i][j].idx[0]]; + CrossProduct(p1,p2,n); + area[i] += Real( Length(n) ); + } + if(area[0]>area[1]) + { + triangles[0]=tIndex[1][0]; + triangles[1]=tIndex[1][1]; + } + else + { + triangles[0]=tIndex[0][0]; + triangles[1]=tIndex[0][1]; + } + return; + } + if(bestTriangulation) + delete[] bestTriangulation; + if(midPoint) + delete[] midPoint; + bestTriangulation=NULL; + midPoint=NULL; + size_t eCount=vertices.size(); + bestTriangulation=new Real[eCount*eCount]; + midPoint=new int[eCount*eCount]; + for(size_t i=0;i +Real MinimalAreaTriangulation::GetArea(const std::vector >& vertices) +{ + if(bestTriangulation) + delete[] bestTriangulation; + if(midPoint) + delete[] midPoint; + bestTriangulation=NULL; + midPoint=NULL; + int eCount=vertices.size(); + bestTriangulation=new double[eCount*eCount]; + midPoint=new int[eCount*eCount]; + for(int i=0;i +void MinimalAreaTriangulation::GetTriangulation(const size_t& i,const size_t& j,const std::vector >& vertices,std::vector& triangles) +{ + TriangleIndex tIndex; + size_t eCount=vertices.size(); + size_t ii=i; + if(i=ii) + return; + ii=midPoint[i*eCount+j]; + if(ii>=0) + { + tIndex.idx[0] = int( i ); + tIndex.idx[1] = int( j ); + tIndex.idx[2] = int( ii ); + triangles.push_back(tIndex); + GetTriangulation(i,ii,vertices,triangles); + GetTriangulation(ii,j,vertices,triangles); + } +} + +template +Real MinimalAreaTriangulation::GetArea(const size_t& i,const size_t& j,const std::vector >& vertices) +{ + Real a=FLT_MAX,temp; + size_t eCount=vertices.size(); + size_t idx=i*eCount+j; + size_t ii=i; + if(i=ii) + { + bestTriangulation[idx]=0; + return 0; + } + if(midPoint[idx]!=-1) + return bestTriangulation[idx]; + int mid=-1; + for(size_t r=j+1;r p,p1,p2; + p1=vertices[i]-vertices[rr]; + p2=vertices[j]-vertices[rr]; + CrossProduct(p1,p2,p); + temp = Real( Length(p) ); + if(bestTriangulation[idx1]>=0) + { + temp+=bestTriangulation[idx1]; + if(temp>a) + continue; + if(bestTriangulation[idx2]>0) + temp+=bestTriangulation[idx2]; + else + temp+=GetArea(rr,j,vertices); + } + else + { + if(bestTriangulation[idx2]>=0) + temp+=bestTriangulation[idx2]; + else + temp+=GetArea(rr,j,vertices); + if(temp>a) + continue; + temp+=GetArea(i,rr,vertices); + } + + if(temp +class MemoryInfo +{ +public: + size_t TotalPhysicalMemory; + size_t FreePhysicalMemory; + size_t TotalSwapSpace; + size_t FreeSwapSpace; + size_t TotalVirtualAddressSpace; + size_t FreeVirtualAddressSpace; + size_t PageSize; + + void set(void){ + MEMORYSTATUSEX Mem; + SYSTEM_INFO Info; + ZeroMemory( &Mem, sizeof(Mem)); + ZeroMemory( &Info, sizeof(Info)); + Mem.dwLength = sizeof(Mem); + ::GlobalMemoryStatusEx( &Mem ); + ::GetSystemInfo( &Info ); + + TotalPhysicalMemory = (size_t)Mem.ullTotalPhys; + FreePhysicalMemory = (size_t)Mem.ullAvailPhys; + TotalSwapSpace = (size_t)Mem.ullTotalPageFile; + FreeSwapSpace = (size_t)Mem.ullAvailPageFile; + TotalVirtualAddressSpace = (size_t)Mem.ullTotalVirtual; + FreeVirtualAddressSpace = (size_t)Mem.ullAvailVirtual; + PageSize = (size_t)Info.dwPageSize; + } + size_t usage(void) const {return TotalVirtualAddressSpace-FreeVirtualAddressSpace;} + + static size_t Usage(void){ + MEMORY_BASIC_INFORMATION mbi; + size_t dwMemUsed = 0; + PVOID pvAddress = 0; + + + memset(&mbi, 0, sizeof(MEMORY_BASIC_INFORMATION)); + while(VirtualQuery(pvAddress, &mbi, sizeof(MEMORY_BASIC_INFORMATION)) == sizeof(MEMORY_BASIC_INFORMATION)){ + if(mbi.State == MEM_COMMIT && mbi.Type == MEM_PRIVATE){dwMemUsed += mbi.RegionSize;} + pvAddress = ((BYTE*)mbi.BaseAddress) + mbi.RegionSize; + } + return dwMemUsed; + } +}; + +#else // !WIN32 + +#ifndef __APPLE__ // Linux variants + +#include +#include + +class MemoryInfo +{ + public: + static size_t Usage(void) + { + FILE* f = fopen("/proc/self/stat","rb"); + + int d; + long ld; + unsigned long lu; + unsigned long long llu; + char* s; + char c; + + int pid; + unsigned long vm; + + int n = fscanf(f, "%d %s %c %d %d %d %d %d %lu %lu %lu %lu %lu %lu %lu %ld %ld %ld %ld %d %ld %llu %lu %ld %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %d %d %lu %lu" + ,&pid ,&s ,&c ,&d ,&d ,&d ,&d ,&d ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&ld ,&ld ,&ld ,&ld ,&d ,&ld ,&llu ,&vm ,&ld ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&d ,&d ,&lu ,&lu ); + + fclose(f); +/* +pid %d +comm %s +state %c +ppid %d +pgrp %d +session %d +tty_nr %d +tpgid %d +flags %lu +minflt %lu +cminflt %lu +majflt %lu +cmajflt %lu +utime %lu +stime %lu +cutime %ld +cstime %ld +priority %ld +nice %ld +0 %ld +itrealvalue %ld +starttime %lu +vsize %lu +rss %ld +rlim %lu +startcode %lu +endcode %lu +startstack %lu +kstkesp %lu +kstkeip %lu +signal %lu +blocked %lu +sigignore %lu +sigcatch %lu +wchan %lu +nswap %lu +cnswap %lu +exit_signal %d +processor %d +rt_priority %lu (since kernel 2.5.19) +policy %lu (since kernel 2.5.19) +*/ + return vm; + } + +}; +#else // __APPLE__: has no "/proc" pseudo-file system + +// Thanks to David O'Gwynn for providing this fix. +// This comes from a post by Michael Knight: +// +// http://miknight.blogspot.com/2005/11/resident-set-size-in-mac-os-x.html + +#include +#include +#include +#include +#include +#include +#include + +void getres(task_t task, unsigned long *rss, unsigned long *vs) +{ + struct task_basic_info t_info; + mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; + + task_info(task, TASK_BASIC_INFO, (task_info_t)&t_info, &t_info_count); + *rss = t_info.resident_size; + *vs = t_info.virtual_size; +} + +class MemoryInfo +{ + public: + static size_t Usage(void) + { + unsigned long rss, vs, psize; + task_t task = MACH_PORT_NULL; + + if (task_for_pid(current_task(), getpid(), &task) != KERN_SUCCESS) + abort(); + getres(task, &rss, &vs); + return rss; + } + +}; + +#endif // !__APPLE__ + +#endif // WIN32 + +#endif // MEMORY_USAGE_INCLUDE diff --git a/src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.IsoSurface.inl b/src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.IsoSurface.inl new file mode 100755 index 000000000..657c9da10 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.IsoSurface.inl @@ -0,0 +1,1113 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#include "Octree.h" +#include "time.h" +#include "MemoryUsage.h" +#include "MAT.h" + + + +template< class Real , int Degree > +template< class Vertex > +Octree< Real , Degree >::SliceValues< Vertex >::SliceValues( void ) +{ + _oldCCount = _oldECount = _oldFCount = _oldNCount = 0; + cornerValues = NullPointer< Real >() ; cornerNormals = NullPointer< Point3D< Real > >() ; cornerSet = NullPointer< char >(); + edgeKeys = NullPointer< long long >() ; edgeSet = NullPointer< char >(); + faceEdges = NullPointer< FaceEdges >() ; faceSet = NullPointer< char >(); + mcIndices = NullPointer< char >(); +} +template< class Real , int Degree > +template< class Vertex > +Octree< Real , Degree >::SliceValues< Vertex >::~SliceValues( void ) +{ + _oldCCount = _oldECount = _oldFCount = _oldNCount = 0; + FreePointer( cornerValues ) ; FreePointer( cornerNormals ) ; FreePointer( cornerSet ); + FreePointer( edgeKeys ) ; FreePointer( edgeSet ); + FreePointer( faceEdges ) ; FreePointer( faceSet ); + FreePointer( mcIndices ); +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::SliceValues< Vertex >::reset( bool nonLinearFit ) +{ + faceEdgeMap.clear() , edgeVertexMap.clear() , vertexPairMap.clear(); + + if( _oldNCount0 ) mcIndices = AllocPointer< char >( _oldNCount ); + } + if( _oldCCount0 ) + { + cornerValues = AllocPointer< Real >( _oldCCount ); + if( nonLinearFit ) cornerNormals = AllocPointer< Point3D< Real > >( _oldCCount ); + cornerSet = AllocPointer< char >( _oldCCount ); + } + } + if( _oldECount( _oldECount ); + edgeSet = AllocPointer< char >( _oldECount ); + } + if( _oldFCount( _oldFCount ); + faceSet = AllocPointer< char >( _oldFCount ); + } + + if( sliceData.cCount>0 ) memset( cornerSet , 0 , sizeof( char ) * sliceData.cCount ); + if( sliceData.eCount>0 ) memset( edgeSet , 0 , sizeof( char ) * sliceData.eCount ); + if( sliceData.fCount>0 ) memset( faceSet , 0 , sizeof( char ) * sliceData.fCount ); +} +template< class Real , int Degree > +template< class Vertex > +Octree< Real , Degree >::XSliceValues< Vertex >::XSliceValues( void ) +{ + _oldECount = _oldFCount = 0; + edgeKeys = NullPointer< long long >() ; edgeSet = NullPointer< char >(); + faceEdges = NullPointer< FaceEdges >() ; faceSet = NullPointer< char >(); +} +template< class Real , int Degree > +template< class Vertex > +Octree< Real , Degree >::XSliceValues< Vertex >::~XSliceValues( void ) +{ + _oldECount = _oldFCount = 0; + FreePointer( edgeKeys ) ; FreePointer( edgeSet ); + FreePointer( faceEdges ) ; FreePointer( faceSet ); +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::XSliceValues< Vertex >::reset( void ) +{ + faceEdgeMap.clear() , edgeVertexMap.clear() , vertexPairMap.clear(); + + if( _oldECount( _oldECount ); + edgeSet = AllocPointer< char >( _oldECount ); + } + if( _oldFCount( _oldFCount ); + faceSet = AllocPointer< char >( _oldFCount ); + } + if( xSliceData.eCount>0 ) memset( edgeSet , 0 , sizeof( char ) * xSliceData.eCount ); + if( xSliceData.fCount>0 ) memset( faceSet , 0 , sizeof( char ) * xSliceData.fCount ); +} + +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights , ConstPointer( Real ) solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit , bool addBarycenter , bool polygonMesh ) +{ + typename BSplineData< Degree >::template CornerEvaluator< 2 > evaluator; + _fData.setCornerEvaluator( evaluator , 0 , 0 , _boundaryType==0 ); + + int maxDepth = tree.maxDepth(); + + std::vector< Real > coarseSolution( _sNodes.nodeCount[maxDepth] , 0 ); +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodes.nodeCount[_minDepth] ; i<_sNodes.nodeCount[maxDepth] ; i++ ) coarseSolution[i] = solution[i]; + for( int d=_minDepth ; d vStencils( maxDepth+1 ); + std::vector< CornerNormalStencil > nStencils( maxDepth+1 ); + for( int d=_minDepth ; d<=maxDepth ; d++ ) + { + SetCornerEvaluationStencil ( evaluator , d , vStencils[d].stencil ); + SetCornerEvaluationStencils( evaluator , d , vStencils[d].stencils ); + SetCornerNormalEvaluationStencil ( evaluator , d , nStencils[d].stencil ); + SetCornerNormalEvaluationStencils( evaluator , d , nStencils[d].stencils ); + } + int vertexOffset = 0; + std::vector< SlabValues< Vertex > > slabValues( maxDepth+1 ); + + // Initialize the back slice + for( int d=maxDepth ; d>=_minDepth ; d-- ) + { + _sNodes.setSliceTableData ( slabValues[d].sliceValues(0).sliceData , d , 0 , threads ); + _sNodes.setSliceTableData ( slabValues[d].sliceValues(1).sliceData , d , 1 , threads ); + _sNodes.setXSliceTableData( slabValues[d].xSliceValues(0).xSliceData , d , 0 , threads ); + slabValues[d].sliceValues (0).reset( nonLinearFit ); + slabValues[d].sliceValues (1).reset( nonLinearFit ); + slabValues[d].xSliceValues(0).reset( ); + } + for( int d=maxDepth ; d>=_minDepth ; d-- ) + { + // Copy edges from finer + if( d=_minDepth ; d-- , o>>=1 ) + { + // Copy edges from finer (required to ensure we correctly track edge cancellations) + if( d=_minDepth ; d-- , o>>=1 ) + { + // Initialize for the next pass + if( o<(1< +Real Octree< Real , Degree >::GetIsoValue( ConstPointer( Real ) solution , const std::vector< Real >& centerWeights ) +{ + Real isoValue=0 , weightSum=0; + int maxDepth = tree.maxDepth(); + + typename BSplineData< Degree >::template CenterEvaluator< 1 > evaluator; + _fData.setCenterEvaluator( evaluator , 0 , 0 , _boundaryType==0 ); + std::vector< CenterValueStencil > vStencils( maxDepth+1 ); + for( int d=_minDepth ; d<=maxDepth ; d++ ) + { + SetCenterEvaluationStencil ( evaluator , d , vStencils[d].stencil ); + SetCenterEvaluationStencils( evaluator , d , vStencils[d].stencils ); + } + std::vector< Real > metSolution( _sNodes.nodeCount[maxDepth] , 0 ); + std::vector< Real > centerValues( _sNodes.nodeCount[maxDepth+1] ); +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodes.nodeCount[_minDepth] ; i<_sNodes.nodeCount[maxDepth] ; i++ ) metSolution[i] = solution[i]; + for( int d=_minDepth ; d=_minDepth ; d-- ) + { + typename TreeOctNode::ConstNeighborKey3 nKey; + nKey.set( d ); +#pragma omp parallel for num_threads( threads ) reduction( + : isoValue , weightSum ) firstprivate( nKey ) + for( int i=_sNodes.nodeCount[d] ; i<_sNodes.nodeCount[d+1] ; i++ ) + { + TreeOctNode* node = _sNodes.treeNodes[i]; + Real value = Real(0); + if( node->children ) + { + for( int c=0 ; cchildren[c].nodeData.nodeIndex ]; + value /= Cube::CORNERS; + } + else + { + nKey.getNeighbors( node ); + int c=0 , x , y , z; + if( node->parent ) c = int( node - node->parent->children ); + Cube::FactorCornerIndex( c , x , y , z ); + + int d , off[3]; + node->depthAndOffset( d , off ); + int o = _boundaryType==0 ? (1<<(d-2)) : 0; + int mn = 2+o , mx = (1<=mn && off[0]=mn && off[1]=mn && off[2]nodeData.nodeIndex ]; + if( w!=0 ) isoValue += value * w , weightSum += w; + } + } + if( _boundaryType==-1 ) return isoValue/weightSum - Real(0.5); + else return isoValue/weightSum; +} + +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::SetSliceIsoCorners( ConstPointer( Real ) solution , ConstPointer( Real ) coarseSolution , Real isoValue , int depth , int slice , std::vector< SlabValues< Vertex > >& slabValues , const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 > vStencil[8] , const Stencil< double , 3 > vStencils[8][8] , const Stencil< Point3D< double > , 3 > nStencil[8] , const Stencil< Point3D< double > , 3 > nStencils[8][8] , int threads ) +{ + if( slice>0 ) SetSliceIsoCorners( solution , coarseSolution , isoValue , depth , slice , 1 , slabValues , evaluator , vStencil , vStencils , nStencil , nStencils , threads ); + if( slice<(1< +template< class Vertex > +void Octree< Real , Degree >::SetSliceIsoCorners( ConstPointer( Real ) solution , ConstPointer( Real ) coarseSolution , Real isoValue , int depth , int slice , int z , std::vector< SlabValues< Vertex > >& slabValues , const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 > vStencil[8] , const Stencil< double , 3 > vStencils[8][8] , const Stencil< Point3D< double > , 3 > nStencil[8] , const Stencil< Point3D< double > , 3 > nStencils[8][8] , int threads ) +{ + typename Octree< Real , Degree >::template SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice ); + Real squareValues[ Square::CORNERS ]; + typename TreeOctNode::ConstNeighborKey3 nKey; + nKey.set( depth ); + { +#pragma omp parallel for num_threads( threads ) firstprivate( nKey , squareValues ) + for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z+1] ; i++ ) + { + TreeOctNode* leaf = _sNodes.treeNodes[i]; + if( !leaf->children ) + { + const typename SortedTreeNodes::SquareCornerIndices& cIndices = sValues.sliceData.cornerIndices( leaf ); + + bool isInterior; + { + int d , off[3]; + leaf->depthAndOffset( d , off ); + int o = _boundaryType==0 ? (1<<(d-2)) : 0; + int mn = 2+o , mx = (1<=mn && off[0]=mn && off[1]=mn && off[2] > p = getCornerValueAndNormal( nKey , leaf , cc , solution , coarseSolution , evaluator , vStencil[cc] , vStencils[cc] , nStencil[cc] , nStencils[cc] , isInterior ); + sValues.cornerValues[vIndex] = p.first , sValues.cornerNormals[vIndex] = p.second; + } + else sValues.cornerValues[vIndex] = getCornerValue( nKey , leaf , cc , solution , coarseSolution , evaluator , vStencil[cc] , vStencils[cc] , isInterior ); + sValues.cornerSet[vIndex] = 1; + } + squareValues[fc] = sValues.cornerValues[ vIndex ]; + TreeOctNode* node = leaf; + int _depth = depth , _slice = slice; + while( node->parent && (node-node->parent->children)==cc ) + { + node = node->parent , _depth-- , _slice >>= 1; + typename Octree< Real , Degree >::template SliceValues< Vertex >& _sValues = slabValues[_depth].sliceValues( _slice ); + const typename SortedTreeNodes::SquareCornerIndices& _cIndices = _sValues.sliceData.cornerIndices( node ); + int _vIndex = _cIndices[fc]; + _sValues.cornerValues[_vIndex] = sValues.cornerValues[vIndex]; + if( _sValues.cornerNormals ) _sValues.cornerNormals[_vIndex] = sValues.cornerNormals[vIndex]; + _sValues.cornerSet[_vIndex] = 1; + } + } + sValues.mcIndices[ i - sValues.sliceData.nodeOffset ] = MarchingSquares::GetIndex( squareValues , isoValue ); + } + } + } +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , Real isoValue , int depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + if( slice>0 ) SetSliceIsoVertices( kernelDensityWeight , isoValue , depth , slice , 1 , vOffset , mesh , slabValues , threads ); + if( slice<(1< +template< class Vertex > +void Octree< Real , Degree >::SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , Real isoValue , int depth , int slice , int z , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + typename Octree< Real , Degree >::template SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice ); + typename TreeOctNode::ConstNeighborKey3 nKey; + nKey.set( depth ); +#pragma omp parallel for num_threads( threads ) firstprivate( nKey ) + for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z+1] ; i++ ) + { + TreeOctNode* leaf = _sNodes.treeNodes[i]; + if( !leaf->children ) + { + int idx = i - sValues.sliceData.nodeOffset; + const typename SortedTreeNodes::SquareEdgeIndices& eIndices = sValues.sliceData.edgeIndices( leaf ); + if( MarchingSquares::HasRoots( sValues.mcIndices[idx] ) ) + { + nKey.getNeighbors( leaf ); + for( int e=0 ; e hashed_vertex; +#pragma omp critical (add_point_access) + { + if( !sValues.edgeSet[vIndex] ) + { + mesh.addOutOfCorePoint( vertex ); + sValues.edgeSet[ vIndex ] = 1; + sValues.edgeKeys[ vIndex ] = key; + sValues.edgeVertexMap[key] = hashed_vertex = std::pair< int , Vertex >( vOffset , vertex ); + vOffset++; + stillOwner = true; + } + } + if( stillOwner ) + { + // We only need to pass the iso-vertex down if the edge it lies on is adjacent to a coarser leaf + bool isNeeded; + switch( o ) + { + case 0: isNeeded = ( nKey.neighbors[depth].neighbors[1][2*y][1]==NULL || nKey.neighbors[depth].neighbors[1][2*y][2*z]==NULL || nKey.neighbors[depth].neighbors[1][1][2*z]==NULL ) ; break; + case 1: isNeeded = ( nKey.neighbors[depth].neighbors[2*y][1][1]==NULL || nKey.neighbors[depth].neighbors[2*y][1][2*z]==NULL || nKey.neighbors[depth].neighbors[1][1][2*z]==NULL ) ; break; + } + if( isNeeded ) + { + int f[2]; + Cube::FacesAdjacentToEdge( Cube::EdgeIndex( o , y , z ) , f[0] , f[1] ); + for( int k=0 ; k<2 ; k++ ) + { + TreeOctNode* node = leaf; + int _depth = depth , _slice = slice; + bool _isNeeded = isNeeded; + while( _isNeeded = node->parent && Cube::IsFaceCorner( (int)(node-node->parent->children) , f[k] ) ) + { + node = node->parent , _depth-- , _slice >>= 1; + typename Octree< Real , Degree >::template SliceValues< Vertex >& _sValues = slabValues[_depth].sliceValues( _slice ); +#pragma omp critical (add_coarser_point_access) + _sValues.edgeVertexMap[key] = hashed_vertex; + switch( o ) + { + case 0: _isNeeded = ( nKey.neighbors[_depth].neighbors[1][2*y][1]==NULL || nKey.neighbors[_depth].neighbors[1][2*y][2*z]==NULL || nKey.neighbors[_depth].neighbors[1][1][2*z]==NULL ) ; break; + case 1: _isNeeded = ( nKey.neighbors[_depth].neighbors[2*y][1][1]==NULL || nKey.neighbors[_depth].neighbors[2*y][1][2*z]==NULL || nKey.neighbors[_depth].neighbors[1][1][2*z]==NULL ) ; break; + } + } + } + } + } + } + } + } + } + } +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWeight , Real isoValue , int depth , int slab , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + typename Octree< Real , Degree >::template SliceValues< Vertex >& bValues = slabValues[depth].sliceValues ( slab ); + typename Octree< Real , Degree >::template SliceValues< Vertex >& fValues = slabValues[depth].sliceValues ( slab+1 ); + typename Octree< Real , Degree >::template XSliceValues< Vertex >& xValues = slabValues[depth].xSliceValues( slab ); + typename TreeOctNode::ConstNeighborKey3 nKey; + nKey.set( depth ); + +#pragma omp parallel for num_threads( threads ) firstprivate( nKey ) + for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab+1] ; i++ ) + { + TreeOctNode* leaf = _sNodes.treeNodes[i]; + if( !leaf->children ) + { + unsigned char mcIndex = ( bValues.mcIndices[ i - bValues.sliceData.nodeOffset ] ) | ( fValues.mcIndices[ i - fValues.sliceData.nodeOffset ] )<<4; + const typename SortedTreeNodes::SquareCornerIndices& eIndices = xValues.xSliceData.edgeIndices( leaf ); + if( MarchingCubes::HasRoots( mcIndex ) ) + { + nKey.getNeighbors( leaf ); + for( int x=0 ; x<2 ; x++ ) for( int y=0 ; y<2 ; y++ ) + { + int c = Square::CornerIndex( x , y ); + int e = Cube::EdgeIndex( 2 , x , y ); + if( MarchingCubes::HasEdgeRoots( mcIndex , e ) ) + { + int vIndex = eIndices[c]; + if( !xValues.edgeSet[vIndex] ) + { + Vertex vertex; + long long key = VertexData::EdgeIndex( leaf , e , _sNodes.maxDepth ); + GetIsoVertex( kernelDensityWeight , isoValue , nKey , leaf , c , bValues , fValues , vertex ); + vertex.point = vertex.point * _scale + _center; + bool stillOwner = false; + std::pair< int , Vertex > hashed_vertex; +#pragma omp critical (add_x_point_access) + { + if( !xValues.edgeSet[vIndex] ) + { + mesh.addOutOfCorePoint( vertex ); + xValues.edgeSet[ vIndex ] = 1; + xValues.edgeKeys[ vIndex ] = key; + xValues.edgeVertexMap[key] = hashed_vertex = std::pair< int , Vertex >( vOffset , vertex ); + stillOwner = true; + vOffset++; + } + } + if( stillOwner ) + { + // We only need to pass the iso-vertex down if the edge it lies on is adjacent to a coarser leaf + bool isNeeded = ( nKey.neighbors[depth].neighbors[2*x][1][1]==NULL || nKey.neighbors[depth].neighbors[2*x][2*y][1]==NULL || nKey.neighbors[depth].neighbors[1][2*y][1]==NULL ); + if( isNeeded ) + { + int f[2]; + Cube::FacesAdjacentToEdge( e , f[0] , f[1] ); + for( int k=0 ; k<2 ; k++ ) + { + TreeOctNode* node = leaf; + int _depth = depth , _slab = slab; + bool _isNeeded = isNeeded; + while( _isNeeded && node->parent && Cube::IsFaceCorner( (int)(node-node->parent->children) , f[k] ) ) + { + node = node->parent , _depth-- , _slab >>= 1; + typename Octree< Real , Degree >::template XSliceValues< Vertex >& _xValues = slabValues[_depth].xSliceValues( _slab ); +#pragma omp critical (add_x_coarser_point_access) + _xValues.edgeVertexMap[key] = hashed_vertex; + _isNeeded = ( nKey.neighbors[_depth].neighbors[2*x][1][1]==NULL || nKey.neighbors[_depth].neighbors[2*x][2*y][1]==NULL || nKey.neighbors[_depth].neighbors[1][2*y][1]==NULL ); + } + } + } + } + } + } + } + } + } + } +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::CopyFinerSliceIsoEdgeKeys( int depth , int slice , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + if( slice>0 ) CopyFinerSliceIsoEdgeKeys( depth , slice , 1 , slabValues , threads ); + if( slice<(1< +template< class Vertex > +void Octree< Real , Degree >::CopyFinerSliceIsoEdgeKeys( int depth , int slice , int z , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + SliceValues< Vertex >& pSliceValues = slabValues[depth ].sliceValues(slice ); + SliceValues< Vertex >& cSliceValues = slabValues[depth+1].sliceValues(slice<<1); + typename SortedTreeNodes::SliceTableData& pSliceData = pSliceValues.sliceData; + typename SortedTreeNodes::SliceTableData& cSliceData = cSliceValues.sliceData; +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodes.nodeCount[depth] + _sNodes.sliceOffsets[depth][slice-z] ; i<_sNodes.nodeCount[depth] + _sNodes.sliceOffsets[depth][slice-z+1] ; i++ ) + if( _sNodes.treeNodes[i]->children ) + { + typename SortedTreeNodes::SquareEdgeIndices& pIndices = pSliceData.edgeIndices( i ); + // Copy the edges that overlap the coarser edges + for( int orientation=0 ; orientation<2 ; orientation++ ) for( int y=0 ; y<2 ; y++ ) + { + int fe = Square::EdgeIndex( orientation , y ); + int pIndex = pIndices[fe]; + if( !pSliceValues.edgeSet[ pIndex ] ) + { + int ce = Cube::EdgeIndex( orientation , y , z ); + int c1 , c2; + switch( orientation ) + { + case 0: c1 = Cube::CornerIndex( 0 , y , z ) , c2 = Cube::CornerIndex( 1 , y , z ) ; break; + case 1: c1 = Cube::CornerIndex( y , 0 , z ) , c2 = Cube::CornerIndex( y , 1 , z ) ; break; + } + int cIndex1 = cSliceData.edgeIndices( _sNodes.treeNodes[i]->children + c1 )[fe]; + int cIndex2 = cSliceData.edgeIndices( _sNodes.treeNodes[i]->children + c2 )[fe]; + if( cSliceValues.edgeSet[cIndex1] != cSliceValues.edgeSet[cIndex2] ) + { + long long key; + if( cSliceValues.edgeSet[cIndex1] ) key = cSliceValues.edgeKeys[cIndex1]; + else key = cSliceValues.edgeKeys[cIndex2]; + std::pair< int , Vertex > vPair = cSliceValues.edgeVertexMap[key]; +#pragma omp critical ( copy_finer_edge_keys ) + pSliceValues.edgeVertexMap[key] = vPair; + pSliceValues.edgeKeys[pIndex] = key; + pSliceValues.edgeSet[pIndex] = 1; + } + else if( cSliceValues.edgeSet[cIndex1] && cSliceValues.edgeSet[cIndex2] ) + { + long long key1 = cSliceValues.edgeKeys[cIndex1] , key2 = cSliceValues.edgeKeys[cIndex2]; +#pragma omp critical ( set_edge_pairs ) + pSliceValues.vertexPairMap[ key1 ] = key2 , pSliceValues.vertexPairMap[ key2 ] = key1; + + const TreeOctNode* node = _sNodes.treeNodes[i]; + int _depth = depth , _slice = slice; + while( node->parent && Cube::IsEdgeCorner( (int)( node - node->parent->children ) , ce ) ) + { + node = node->parent , _depth-- , _slice >>= 1; + SliceValues< Vertex >& _pSliceValues = slabValues[_depth].sliceValues(_slice); +#pragma omp critical ( set_edge_pairs ) + _pSliceValues.vertexPairMap[ key1 ] = key2 , _pSliceValues.vertexPairMap[ key2 ] = key1; + } + } + } + } + } +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::CopyFinerXSliceIsoEdgeKeys( int depth , int slab , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + XSliceValues< Vertex >& pSliceValues = slabValues[depth ].xSliceValues(slab); + XSliceValues< Vertex >& cSliceValues0 = slabValues[depth+1].xSliceValues( (slab<<1)|0 ); + XSliceValues< Vertex >& cSliceValues1 = slabValues[depth+1].xSliceValues( (slab<<1)|1 ); + typename SortedTreeNodes::XSliceTableData& pSliceData = pSliceValues.xSliceData; + typename SortedTreeNodes::XSliceTableData& cSliceData0 = cSliceValues0.xSliceData; + typename SortedTreeNodes::XSliceTableData& cSliceData1 = cSliceValues1.xSliceData; +#pragma omp parallel for num_threads( threads ) + for( int i=_sNodes.nodeCount[depth] + _sNodes.sliceOffsets[depth][slab] ; i<_sNodes.nodeCount[depth] + _sNodes.sliceOffsets[depth][slab+1] ; i++ ) + if( _sNodes.treeNodes[i]->children ) + { + typename SortedTreeNodes::SquareCornerIndices& pIndices = pSliceData.edgeIndices( i ); + for( int x=0 ; x<2 ; x++ ) for( int y=0 ; y<2 ; y++ ) + { + int fc = Square::CornerIndex( x , y ); + int pIndex = pIndices[fc]; + if( !pSliceValues.edgeSet[pIndex] ) + { + int c0 = Cube::CornerIndex( x , y , 0 ) , c1 = Cube::CornerIndex( x , y , 1 ); + int cIndex0 = cSliceData0.edgeIndices( _sNodes.treeNodes[i]->children + c0 )[fc]; + int cIndex1 = cSliceData1.edgeIndices( _sNodes.treeNodes[i]->children + c1 )[fc]; + if( cSliceValues0.edgeSet[cIndex0] != cSliceValues1.edgeSet[cIndex1] ) + { + long long key; + std::pair< int , Vertex > vPair; + if( cSliceValues0.edgeSet[cIndex0] ) key = cSliceValues0.edgeKeys[cIndex0] , vPair = cSliceValues0.edgeVertexMap[key]; + else key = cSliceValues1.edgeKeys[cIndex1] , vPair = cSliceValues1.edgeVertexMap[key]; +#pragma omp critical ( copy_finer_x_edge_keys ) + pSliceValues.edgeVertexMap[key] = vPair; + pSliceValues.edgeKeys[ pIndex ] = key; + pSliceValues.edgeSet[ pIndex ] = 1; + } + else if( cSliceValues0.edgeSet[cIndex0] && cSliceValues1.edgeSet[cIndex1] ) + { + long long key0 = cSliceValues0.edgeKeys[cIndex0] , key1 = cSliceValues1.edgeKeys[cIndex1]; +#pragma omp critical ( set_x_edge_pairs ) + pSliceValues.vertexPairMap[ key0 ] = key1 , pSliceValues.vertexPairMap[ key1 ] = key0; + const TreeOctNode* node = _sNodes.treeNodes[i]; + int _depth = depth , _slab = slab , ce = Cube::CornerIndex( 2 , x , y ); + while( node->parent && Cube::IsEdgeCorner( (int)( node - node->parent->children ) , ce ) ) + { + node = node->parent , _depth-- , _slab>>= 1; + SliceValues< Vertex >& _pSliceValues = slabValues[_depth].sliceValues(_slab); +#pragma omp critical ( set_x_edge_pairs ) + _pSliceValues.vertexPairMap[ key0 ] = key1 , _pSliceValues.vertexPairMap[ key1 ] = key0; + } + } + } + } + } +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::SetSliceIsoEdges( int depth , int slice , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + if( slice>0 ) SetSliceIsoEdges( depth , slice , 1 , slabValues , threads ); + if( slice<(1< +template< class Vertex > +void Octree< Real , Degree >::SetSliceIsoEdges( int depth , int slice , int z , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + typename Octree< Real , Degree >::template SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice ); + int isoEdges[ 2 * MarchingSquares::MAX_EDGES ]; + typename TreeOctNode::ConstNeighborKey3 nKey; + nKey.set( depth ); +#pragma omp parallel for num_threads( threads ) firstprivate( nKey , isoEdges ) + for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slice-z+1] ; i++ ) + { + TreeOctNode* leaf = _sNodes.treeNodes[i]; + if( !leaf->children ) + { + int idx = i - sValues.sliceData.nodeOffset; + const typename SortedTreeNodes::SquareEdgeIndices& eIndices = sValues.sliceData.edgeIndices( leaf ); + const typename SortedTreeNodes::SquareFaceIndices& fIndices = sValues.sliceData.faceIndices( leaf ); + unsigned char mcIndex = sValues.mcIndices[idx]; + if( !sValues.faceSet[ fIndices[0] ] ) + { + nKey.getNeighbors( leaf ); + if( !nKey.neighbors[depth].neighbors[1][1][2*z] || !nKey.neighbors[depth].neighbors[1][1][2*z]->children ) + { + FaceEdges fe; + fe.count = MarchingSquares::AddEdgeIndices( mcIndex , isoEdges ); + for( int j=0 ; j edges; + edges.resize( fe.count ); + for( int j=0 ; jparent && Cube::IsFaceCorner( (int)(node-node->parent->children) , f ) ) + { + node = node->parent , _depth-- , _slice >>= 1; + if( nKey.neighbors[_depth].neighbors[1][1][2*z] && nKey.neighbors[_depth].neighbors[1][1][2*z]->children ) break; + long long key = VertexData::FaceIndex( node , f , _sNodes.maxDepth ); +#pragma omp critical( add_iso_edge_access ) + { + typename Octree< Real , Degree >::template SliceValues< Vertex >& _sValues = slabValues[_depth].sliceValues( _slice ); + typename hash_map< long long , std::vector< IsoEdge > >::iterator iter = _sValues.faceEdgeMap.find(key); + if( iter==_sValues.faceEdgeMap.end() ) _sValues.faceEdgeMap[key] = edges; + else for( int j=0 ; jsecond.push_back( fe.edges[j] ); + } + } + } + } + } + } +} +template< class Real , int Degree > +template< class Vertex > +void Octree< Real , Degree >::SetXSliceIsoEdges( int depth , int slab , std::vector< SlabValues< Vertex > >& slabValues , int threads ) +{ + typename Octree< Real , Degree >::template SliceValues< Vertex >& bValues = slabValues[depth].sliceValues ( slab ); + typename Octree< Real , Degree >::template SliceValues< Vertex >& fValues = slabValues[depth].sliceValues ( slab+1 ); + typename Octree< Real , Degree >::template XSliceValues< Vertex >& xValues = slabValues[depth].xSliceValues( slab ); + int isoEdges[ 2 * MarchingSquares::MAX_EDGES ]; + typename TreeOctNode::ConstNeighborKey3 nKey; + nKey.set( depth ); +#pragma omp parallel for num_threads( threads ) firstprivate( nKey , isoEdges ) + for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][slab+1] ; i++ ) + { + TreeOctNode* leaf = _sNodes.treeNodes[i]; + if( !leaf->children ) + { + const typename SortedTreeNodes::SquareCornerIndices& cIndices = xValues.xSliceData.edgeIndices( leaf ); + const typename SortedTreeNodes::SquareEdgeIndices& eIndices = xValues.xSliceData.faceIndices( leaf ); + unsigned char mcIndex = ( bValues.mcIndices[ i - bValues.sliceData.nodeOffset ] ) | ( fValues.mcIndices[ i - fValues.sliceData.nodeOffset ]<<4 ); + { + nKey.getNeighbors( leaf ); + for( int o=0 ; o<2 ; o++ ) for( int x=0 ; x<2 ; x++ ) + { + int e = Square::EdgeIndex( o , x ); + int f = Cube::FaceIndex( 1-o , x ); + unsigned char _mcIndex = MarchingCubes::GetFaceIndex( mcIndex , f ); + int xx = o==1 ? 2*x : 1 , yy = o==0 ? 2*x : 1 , zz = 1; + if( !xValues.faceSet[ eIndices[e] ] && ( !nKey.neighbors[depth].neighbors[xx][yy][zz] || !nKey.neighbors[depth].neighbors[xx][yy][zz]->children ) ) + { + FaceEdges fe; + fe.count = MarchingSquares::AddEdgeIndices( _mcIndex , isoEdges ); + for( int j=0 ; j::template SliceValues< Vertex >& sValues = (_x==0) ? bValues : fValues; + int idx = sValues.sliceData.edgeIndices(i)[ Square::EdgeIndex(o,x) ]; + if( !sValues.edgeSet[ idx ] ) fprintf( stderr , "[ERROR] Edge not set 5: %d / %d\n" , slab , 1< edges; + edges.resize( fe.count ); + for( int j=0 ; jparent && Cube::IsFaceCorner( (int)(node-node->parent->children) , f ) ) + { + node = node->parent , _depth-- , _slab >>= 1; + if( nKey.neighbors[_depth].neighbors[xx][yy][zz] && nKey.neighbors[_depth].neighbors[xx][yy][zz]->children ) break; + long long key = VertexData::FaceIndex( node , f , _sNodes.maxDepth ); +#pragma omp critical( add_x_iso_edge_access ) + { + typename Octree< Real , Degree >::template XSliceValues< Vertex >& _xValues = slabValues[_depth].xSliceValues( _slab ); + typename hash_map< long long , std::vector< IsoEdge > >::iterator iter = _xValues.faceEdgeMap.find(key); + if( iter==_xValues.faceEdgeMap.end() ) _xValues.faceEdgeMap[key] = edges; + else for( int j=0 ; jsecond.push_back( fe.edges[j] ); + } + } + } + } + } + } + } +} +template< class Real , int Degree > +template< class Vertex > +int Octree< Real , Degree >::SetIsoSurface( int depth , int offset , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , const XSliceValues< Vertex >& xValues , CoredMeshData< Vertex >& mesh , bool polygonMesh , bool addBarycenter , int& vOffset , int threads ) +{ + typename TreeOctNode::ConstNeighborKey3 nKey; + std::vector< std::pair< int , Vertex > > polygon; + std::vector< IsoEdge > edges; + nKey.set( depth ); +#pragma omp parallel for num_threads( threads ) firstprivate( nKey , edges ) + for( int i=_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][offset] ; i<_sNodes.nodeCount[depth]+_sNodes.sliceOffsets[depth][offset+1] ; i++ ) + { + TreeOctNode* leaf = _sNodes.treeNodes[i]; + if( !leaf->children ) + { + edges.clear(); + unsigned char mcIndex = ( bValues.mcIndices[ i - bValues.sliceData.nodeOffset ] ) | ( fValues.mcIndices[ i - fValues.sliceData.nodeOffset ]<<4 ); + // [WARNING] Just because the node looks empty doesn't mean it doesn't get eges from finer neighbors + { + // Gather the edges from the faces (with the correct orientation) + for( int f=0 ; f& sValues = (o==0) ? bValues : fValues; + int fIdx = sValues.sliceData.faceIndices(i)[0]; + if( sValues.faceSet[fIdx] ) + { + const FaceEdges& fe = sValues.faceEdges[ fIdx ]; + for( int j=0 ; j >::const_iterator iter = sValues.faceEdgeMap.find( key ); + if( iter!=sValues.faceEdgeMap.end() ) + { + const std::vector< IsoEdge >& _edges = iter->second; + for( int j=0 ; j<_edges.size() ; j++ ) edges.push_back( IsoEdge( _edges[j][flip] , _edges[j][1-flip] ) ); + } + else fprintf( stderr , "[ERROR] Invalid faces: %d %d %d\n" , i , d , o ) , exit( 0 ); + } + } + else + { + int fIdx = xValues.xSliceData.faceIndices(i)[ Square::EdgeIndex( 1-d , o ) ]; + if( xValues.faceSet[fIdx] ) + { + const FaceEdges& fe = xValues.faceEdges[ fIdx ]; + for( int j=0 ; j >::const_iterator iter = xValues.faceEdgeMap.find( key ); + if( iter!=xValues.faceEdgeMap.end() ) + { + const std::vector< IsoEdge >& _edges = iter->second; + for( int j=0 ; j<_edges.size() ; j++ ) edges.push_back( IsoEdge( _edges[j][flip] , _edges[j][1-flip] ) ); + } + else fprintf( stderr , "[ERROR] Invalid faces: %d %d %d\n" , i , d , o ) , exit( 0 ); + } + } + } + // Get the edge loops + std::vector< std::vector< long long > > loops; + while( edges.size() ) + { + loops.resize( loops.size()+1 ); + IsoEdge edge = edges.back(); + edges.pop_back(); + long long start = edge[0] , current = edge[1]; + while( current!=start ) + { + int idx; + for( idx=0 ; idx::const_iterator iter; + if ( (iter=bValues.vertexPairMap.find(current))!=bValues.vertexPairMap.end() ) loops.back().push_back( current ) , current = iter->second; + else if( (iter=fValues.vertexPairMap.find(current))!=fValues.vertexPairMap.end() ) loops.back().push_back( current ) , current = iter->second; + else if( (iter=xValues.vertexPairMap.find(current))!=xValues.vertexPairMap.end() ) loops.back().push_back( current ) , current = iter->second; + else fprintf( stderr , "[ERROR] Failed to close loop @ depth %d / %d (%d): %lld\n" , depth , _sNodes.maxDepth-1 , i , current ) , exit( 0 ); + } + else + { + loops.back().push_back( current ); + current = edges[idx][1]; + edges[idx] = edges.back() , edges.pop_back(); + } + } + loops.back().push_back( start ); + } + // Add the loops to the mesh + for( int j=0 ; j > polygon( loops[j].size() ); + for( int k=0 ; k >::const_iterator iter; + if ( ( iter=bValues.edgeVertexMap.find( key ) )!=bValues.edgeVertexMap.end() ) polygon[k] = iter->second; + else if( ( iter=fValues.edgeVertexMap.find( key ) )!=fValues.edgeVertexMap.end() ) polygon[k] = iter->second; + else if( ( iter=xValues.edgeVertexMap.find( key ) )!=xValues.edgeVertexMap.end() ) polygon[k] = iter->second; + else fprintf( stderr , "[ERROR] Couldn't find vertex in edge map\n" ) , exit( 0 ); + } + AddIsoPolygons( mesh , polygon , polygonMesh , addBarycenter , vOffset ); + } + } + } + } +} +template< class Real > void SetIsoVertexValue( PlyVertex< float >& vertex , Real value ){ ; } +template< class Real > void SetIsoVertexValue( PlyValueVertex< float >& vertex , Real value ){ vertex.value = float(value); } +template< class Real , int Degree > +template< class Vertex > +bool Octree< Real , Degree >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int edgeIndex , int z , const SliceValues< Vertex >& sValues , Vertex& vertex ) +{ + Point3D< Real > position; + int c0 , c1; + Square::EdgeCorners( edgeIndex , c0 , c1 ); + + Real x0 , x1; + Point3D< Real > n0 , n1; + const typename SortedTreeNodes::SquareCornerIndices& idx = sValues.sliceData.cornerIndices( node ); + x0 = sValues.cornerValues[idx[c0]] , x1 = sValues.cornerValues[idx[c1]]; + if( sValues.cornerNormals ) n0 = sValues.cornerNormals[idx[c0]] , n1 = sValues.cornerNormals[idx[c1]]; + + int o , y; + Square::FactorEdgeIndex( edgeIndex , o , y ); + + Point3D< Real > c; + Real center , width; + node->centerAndWidth( c , width ); + center = c[o]; + for( int i=0 ; i P; + P.coefficients[0] = x0; + P.coefficients[1] = dx0; + P.coefficients[2] = 3*(x1-x0)-dx1-2*dx0; + + Real averageRoot=0; + double roots[2]; + int rCount = 0 , rootCount = P.getSolutions( isoValue , roots , EPSILON ); + for( int i=0 ; i=0 && roots[i]<=1 ) averageRoot += Real( roots[i] ) , rCount++; + if( rCount && sValues.cornerNormals ) averageRoot /= rCount; + else averageRoot = Real( ( x0-isoValue ) / ( x0-x1 ) ); + if( averageRoot<0 || averageRoot>1 ) + { + fprintf( stderr , "[WARNING] Bad average root: %f\n" , averageRoot ); + fprintf( stderr , "\t(%f %f) , (%f %f) (%f)\n" , x0 , x1 , dx0 , dx1 , isoValue ); + if( averageRoot<0 ) averageRoot = 0; + if( averageRoot>1 ) averageRoot = 1; + } + position[o] = Real( center - width/2 + width*averageRoot ); + vertex.point = position; + if( kernelDensityWeights ) + { + Real depth , weight; + const TreeOctNode* temp = node; + while( temp->depth()>_splatDepth ) temp=temp->parent; + GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey3 , _samplesPerNode , depth , weight ); + SetIsoVertexValue( vertex , depth ); + } + return true; +} +template< class Real , int Degree > +template< class Vertex > +bool Octree< Real , Degree >::GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int cornerIndex , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , Vertex& vertex ) +{ + Point3D< Real > position; + + Real x0 , x1; + Point3D< Real > n0 , n1; + const typename SortedTreeNodes::SquareCornerIndices& idx0 = bValues.sliceData.cornerIndices( node ); + const typename SortedTreeNodes::SquareCornerIndices& idx1 = fValues.sliceData.cornerIndices( node ); + x0 = bValues.cornerValues[ idx0[cornerIndex] ] , x1 = fValues.cornerValues[ idx1[cornerIndex] ]; + if( bValues.cornerNormals || fValues.cornerNormals ) n0 = bValues.cornerNormals[ idx0[cornerIndex] ] , n1 = fValues.cornerNormals[ idx1[cornerIndex] ]; + + int x , y; + Square::FactorCornerIndex( cornerIndex , x , y ); + + Point3D< Real > c; + Real center , width; + node->centerAndWidth( c , width ); + center = c[2]; + for( int i=0 ; i P; + P.coefficients[0] = x0; + P.coefficients[1] = dx0; + P.coefficients[2] = 3*(x1-x0)-dx1-2*dx0; + + Real averageRoot=0; + double roots[2]; + int rCount = 0 , rootCount = P.getSolutions( isoValue , roots , EPSILON ); + for( int i=0 ; i=0 && roots[i]<=1 ) averageRoot += Real( roots[i] ) , rCount++; + if( rCount && bValues.cornerNormals && fValues.cornerNormals ) averageRoot /= rCount; + else averageRoot = Real( ( x0-isoValue ) / ( x0-x1 ) ); + if( averageRoot<0 || averageRoot>1 ) + { + fprintf( stderr , "[WARNING] Bad average root: %f\n" , averageRoot ); + fprintf( stderr , "\t(%f %f) , (%f %f) (%f)\n" , x0 , x1 , dx0 , dx1 , isoValue ); + if( averageRoot<0 ) averageRoot = 0; + if( averageRoot>1 ) averageRoot = 1; + } + position[2] = Real( center - width/2 + width*averageRoot ); + vertex.point = position; + if( kernelDensityWeights ) + { + Real depth , weight; + const TreeOctNode* temp = node; + while( temp->depth()>_splatDepth ) temp=temp->parent; + GetSampleDepthAndWeight( kernelDensityWeights , temp , position , neighborKey3 , _samplesPerNode , depth , weight ); + SetIsoVertexValue( vertex , depth ); + } + return true; +} + +template< class Real , int Degree > +template< class Vertex > +int Octree< Real , Degree >::AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< std::pair< int , Vertex > >& polygon , bool polygonMesh , bool addBarycenter , int& vOffset ) +{ + if( polygonMesh ) + { + std::vector< int > vertices( polygon.size() ); + for( int i=0 ; i<(int)polygon.size() ; i++ ) vertices[i] = polygon[polygon.size()-1-i].first; + mesh.addPolygon_s( vertices ); + return 1; + } + if( polygon.size()>3 ) + { + bool isCoplanar = false; + std::vector< int > triangle( 3 ); + + if( addBarycenter ) + for( int i=0 ; i<(int)polygon.size() ; i++ ) + for( int j=0 ; j MAT; + std::vector< Point3D< Real > > vertices; + std::vector< TriangleIndex > triangles; + vertices.resize( polygon.size() ); + // Add the points + for( int i=0 ; i<(int)polygon.size() ; i++ ) vertices[i] = polygon[i].second.point; + MAT.GetTriangulation( vertices , triangles ); + for( int i=0 ; i<(int)triangles.size() ; i++ ) + { + for( int j=0 ; j<3 ; j++ ) triangle[2-j] = polygon[ triangles[i].idx[j] ].first; + mesh.addPolygon_s( triangle ); + } + } + } + else if( polygon.size()==3 ) + { + std::vector< int > vertices( 3 ); + for( int i=0 ; i<3 ; i++ ) vertices[2-i] = polygon[i].first; + mesh.addPolygon_s( vertices ); + } + return (int)polygon.size()-2; +} diff --git a/src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.h b/src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.h new file mode 100755 index 000000000..98ef154db --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/MultiGridOctreeData.h @@ -0,0 +1,438 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#ifndef MULTI_GRID_OCTREE_DATA_INCLUDED +#define MULTI_GRID_OCTREE_DATA_INCLUDED + +#define NEW_CODE 1 + +#define MAX_MEMORY_GB 15 + +#define GRADIENT_DOMAIN_SOLUTION 1 // Given the constraint vector-field V(p), there are two ways to solve for the coefficients, x, of the indicator function + // with respect to the B-spline basis {B_i(p)} + // 1] Find x minimizing: + // || V(p) - \sum_i \nabla x_i B_i(p) ||^2 + // which is solved by the system A_1x = b_1 where: + // A_1[i,j] = < \nabla B_i(p) , \nabla B_j(p) > + // b_1[i] = < \nabla B_i(p) , V(p) > + // 2] Formulate this as a Poisson equation: + // \sum_i x_i \Delta B_i(p) = \nabla \cdot V(p) + // which is solved by the system A_2x = b_2 where: + // A_2[i,j] = - < \Delta B_i(p) , B_j(p) > + // b_2[i] = - < B_i(p) , \nabla \cdot V(p) > + // Although the two system matrices should be the same (assuming that the B_i satisfy dirichlet/neumann boundary conditions) + // the constraint vectors can differ when V does not satisfy the Neumann boundary conditions: + // A_1[i,j] = \int_R < \nabla B_i(p) , \nabla B_j(p) > + // = \int_R [ \nabla \cdot ( B_i(p) \nabla B_j(p) ) - B_i(p) \Delta B_j(p) ] + // = \int_dR < N(p) , B_i(p) \nabla B_j(p) > + A_2[i,j] + // and the first integral is zero if either f_i is zero on the boundary dR or the derivative of B_i across the boundary is zero. + // However, for the constraints we have: + // b_1(i) = \int_R < \nabla B_i(p) , V(p) > + // = \int_R [ \nabla \cdot ( B_i(p) V(p) ) - B_i(p) \nabla \cdot V(p) ] + // = \int_dR < N(p) , B_i(p) V(p) > + b_2[i] + // In particular, this implies that if the B_i satisfy the Neumann boundary conditions (rather than Dirichlet), + // and V is not zero across the boundary, then the two constraints are different. + // Forcing the < V(p) , N(p) > = 0 on the boundary, by killing off the component of the vector-field in the normal direction + // (FORCE_NEUMANN_FIELD), makes the two systems equal, and the value of this flag should be immaterial. + // Note that under interpretation 1, we have: + // \sum_i b_1(i) = < \nabla \sum_ i B_i(p) , V(p) > = 0 + // because the B_i's sum to one. However, in general, we could have + // \sum_i b_2(i) \neq 0. + // This could cause trouble because the constant functions are in the kernel of the matrix A, so CG will misbehave if the constraint + // has a non-zero DC term. (Again, forcing < V(p) , N(p) > = 0 along the boundary resolves this problem.) + +#define FORCE_NEUMANN_FIELD 1 // This flag forces the normal component across the boundary of the integration domain to be zero. + // This should be enabled if GRADIENT_DOMAIN_SOLUTION is not, so that CG doesn't run into trouble. + +#define ROBERTO_TOLDO_FIX 1 + +#if !FORCE_NEUMANN_FIELD +#pragma message( "[WARNING] Not zeroing out normal component on boundary" ) +#endif // !FORCE_NEUMANN_FIELD + +#include "Hash.h" +#include "BSplineData.h" + +class TreeNodeData +{ +public: + static int NodeCount; + int nodeIndex; + + TreeNodeData( void ); + ~TreeNodeData( void ); +}; + +class VertexData +{ + typedef OctNode< TreeNodeData > TreeOctNode; +public: + static const int VERTEX_COORDINATE_SHIFT = ( sizeof( long long ) * 8 ) / 3; + static long long EdgeIndex( const TreeOctNode* node , int eIndex , int maxDepth , int index[DIMENSION] ); + static long long EdgeIndex( const TreeOctNode* node , int eIndex , int maxDepth ); + static long long FaceIndex( const TreeOctNode* node , int fIndex , int maxDepth,int index[DIMENSION] ); + static long long FaceIndex( const TreeOctNode* node , int fIndex , int maxDepth ); + static long long CornerIndex( const TreeOctNode* node , int cIndex , int maxDepth , int index[DIMENSION] ); + static long long CornerIndex( const TreeOctNode* node , int cIndex , int maxDepth ); + static long long CenterIndex( const TreeOctNode* node , int maxDepth , int index[DIMENSION] ); + static long long CenterIndex( const TreeOctNode* node , int maxDepth ); + static long long CornerIndex( int depth , const int offSet[DIMENSION] , int cIndex , int maxDepth , int index[DIMENSION] ); + static long long CenterIndex( int depth , const int offSet[DIMENSION] , int maxDepth , int index[DIMENSION] ); + static long long CornerIndexKey( const int index[DIMENSION] ); +}; + +class SortedTreeNodes +{ + typedef OctNode< TreeNodeData > TreeOctNode; +protected: + void _sortByZCoordinate( void ); +public: + Pointer( TreeOctNode* ) treeNodes; + int *nodeCount; + int maxDepth; + SortedTreeNodes( void ); + ~SortedTreeNodes( void ); + void set( TreeOctNode& root , std::vector< int >* map ); + Pointer( Pointer( int ) ) sliceOffsets; + static int Slices( int depth ); + std::pair< int , int > sliceSpan( int depth , int off , int d ) const; + + template< int Indices > + struct _Indices + { + int idx[Indices]; + _Indices( void ){ memset( idx , -1 , sizeof( int ) * Indices ); } + int& operator[] ( int i ) { return idx[i]; } + const int& operator[] ( int i ) const { return idx[i]; } + }; + typedef _Indices< Square::CORNERS > SquareCornerIndices; + typedef _Indices< Square::EDGES > SquareEdgeIndices; + typedef _Indices< Square::FACES > SquareFaceIndices; + + struct SliceTableData + { + std::vector< SquareCornerIndices > cTable; + std::vector< SquareEdgeIndices > eTable; + std::vector< SquareFaceIndices > fTable; + int cCount , eCount , fCount , nodeOffset , nodeCount; + SliceTableData( void ){ fCount , eCount = cCount = 0; } + ~SliceTableData( void ){ clear(); } + void clear( void ) { cTable.clear() , eTable.clear() , fTable.clear() , fCount = eCount = cCount = 0; } + SquareCornerIndices& cornerIndices( const TreeOctNode* node ); + SquareCornerIndices& cornerIndices( int idx ); + const SquareCornerIndices& cornerIndices( const TreeOctNode* node ) const; + const SquareCornerIndices& cornerIndices( int idx ) const; + SquareEdgeIndices& edgeIndices( const TreeOctNode* node ); + SquareEdgeIndices& edgeIndices( int idx ); + const SquareEdgeIndices& edgeIndices( const TreeOctNode* node ) const; + const SquareEdgeIndices& edgeIndices( int idx ) const; + SquareFaceIndices& faceIndices( const TreeOctNode* node ); + SquareFaceIndices& faceIndices( int idx ); + const SquareFaceIndices& faceIndices( const TreeOctNode* node ) const; + const SquareFaceIndices& faceIndices( int idx ) const; + protected: + std::vector< int > _cMap , _eMap , _fMap; + friend class SortedTreeNodes; + }; + struct XSliceTableData + { + std::vector< SquareCornerIndices > eTable; + std::vector< SquareEdgeIndices > fTable; + int fCount , eCount , nodeOffset , nodeCount; + XSliceTableData( void ){ fCount , eCount = 0; } + ~XSliceTableData( void ){ clear(); } + void clear( void ) { fTable.clear() , eTable.clear() , fCount = eCount = 0; } + SquareCornerIndices& edgeIndices( const TreeOctNode* node ); + SquareCornerIndices& edgeIndices( int idx ); + const SquareCornerIndices& edgeIndices( const TreeOctNode* node ) const; + const SquareCornerIndices& edgeIndices( int idx ) const; + SquareEdgeIndices& faceIndices( const TreeOctNode* node ); + SquareEdgeIndices& faceIndices( int idx ); + const SquareEdgeIndices& faceIndices( const TreeOctNode* node ) const; + const SquareEdgeIndices& faceIndices( int idx ) const; + protected: + std::vector< int > _eMap , _fMap; + friend class SortedTreeNodes; + }; + void setSliceTableData ( SliceTableData& sData , int depth , int offset , int threads ) const; + void setXSliceTableData( XSliceTableData& sData , int depth , int offset , int threads ) const; +}; + + +template< class Real , int Degree > +class Octree +{ + typedef OctNode< TreeNodeData > TreeOctNode; + struct _PointData + { + Point3D< Real > position; + Real weightedCoarserValue; + Real weight; + _PointData( Point3D< Real > p=Point3D< Real >() , Real w=0 ) { position = p , weight = w , weightedCoarserValue = Real(0); } + }; +public: + struct NormalInfo + { + std::vector< int > normalIndices; + std::vector< Point3D< Real > > normals; + int normalIndex( const TreeOctNode* node ) const { return node->nodeData.nodeIndex>=normalIndices.size() ? -1 : normalIndices[ node->nodeData.nodeIndex ]; } + }; + struct PointInfo + { + std::vector< int > pointIndices; + std::vector< _PointData > points; + int pointIndex( const TreeOctNode* node ) const { return node->nodeData.nodeIndex>=pointIndices.size() ? -1 : pointIndices[ node->nodeData.nodeIndex ]; } + }; +protected: + SortedTreeNodes _sNodes; + Real _samplesPerNode; + int _splatDepth; + int _minDepth; + int _fullDepth; + bool _constrainValues; + int _boundaryType; + Real _scale; + Point3D< Real > _center; + std::vector< int > _pointCount; + Real _normalSmooth; + BSplineData< Degree > _fData; + + bool _InBounds( Point3D< Real > ) const; + + double GetLaplacian ( const typename BSplineData< Degree >::Integrator& integrator , int d , const int off1[3] , const int off2[3] , bool childParent ) const; + double GetDivergence1( const typename BSplineData< Degree >::Integrator& integrator , int d , const int off1[3] , const int off2[3] , bool childParent , const Point3D< Real >& normal1 ) const; + double GetDivergence2( const typename BSplineData< Degree >::Integrator& integrator , int d , const int off1[3] , const int off2[3] , bool childParent , const Point3D< Real >& normal2 ) const; + Point3D< double > GetDivergence1( const typename BSplineData< Degree >::Integrator& integrator , int d , const int off1[3] , const int off2[3] , bool childParent ) const; + Point3D< double > GetDivergence2( const typename BSplineData< Degree >::Integrator& integrator , int d , const int off1[3] , const int off2[3] , bool childParent ) const; + + template< class C , int N > struct Stencil{ C values[N][N][N]; }; + struct CenterValueStencil + { + Stencil< double , 3 > stencil; + Stencil< double , 3 > stencils[8]; + }; + struct CornerValueStencil + { + Stencil< double , 3 > stencil[8]; + Stencil< double , 3 > stencils[8][8]; + }; + struct CornerNormalStencil + { + Stencil< Point3D< double > , 3 > stencil[8]; + Stencil< Point3D< double > , 3 > stencils[8][8]; + }; + + void _setMultiColorIndices( int start , int end , std::vector< std::vector< int > >& indices ) const; + int _SolveSystemGS( PointInfo& pointInfo , int depth , const typename BSplineData< Degree >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual=false , double* bNorm2=NULL , double* inRNorm2=NULL , double* outRNorm2=NULL , bool forceSilent=false ); + int _SolveSystemCG( PointInfo& pointInfo , int depth , const typename BSplineData< Degree >::Integrator& integrator , const SortedTreeNodes& sNodes , Pointer( Real ) solution , Pointer( Real ) constraints , Pointer( Real ) metSolutionConstraints , int iters , bool coarseToFine , bool showResidual=false , double* bNorm2=NULL , double* inRNorm2=NULL , double* outRNorm2=NULL , double accuracy=0 ); + + int GetMatrixRowSize( const typename TreeOctNode::Neighbors5& neighbors5 , bool symmetric ) const; + int SetMatrixRow( const PointInfo& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , Pointer( MatrixEntry< Real > ) row , int offset , const typename BSplineData< Degree >::Integrator& integrator , const Stencil< double , 5 >& stencil , bool symmetric ) const; + + void SetDivergenceStencil ( int depth , const typename BSplineData< Degree >::Integrator& integrator , Stencil< Point3D< double > , 5 >& stencil , bool scatter ) const; + void SetDivergenceStencils( int depth , const typename BSplineData< Degree >::Integrator& integrator , Stencil< Point3D< double > , 5 > stencil[2][2][2] , bool scatter ) const; + void SetLaplacianStencil ( int depth , const typename BSplineData< Degree >::Integrator& integrator , Stencil< double , 5 >& stencil ) const; + void SetLaplacianStencils ( int depth , const typename BSplineData< Degree >::Integrator& integrator , Stencil< double , 5 > stencil[2][2][2] ) const; + void SetCenterEvaluationStencil ( const typename BSplineData< Degree >::template CenterEvaluator< 1 >& evaluator , int depth , Stencil< double , 3 >& stencil ) const; + void SetCenterEvaluationStencils( const typename BSplineData< Degree >::template CenterEvaluator< 1 >& evaluator , int depth , Stencil< double , 3 > stencil[8] ) const; + void SetCornerEvaluationStencil ( const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , int depth , Stencil< double , 3 > stencil [8] ) const; + void SetCornerEvaluationStencils( const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , int depth , Stencil< double , 3 > stencils[8][8] ) const; + void SetCornerNormalEvaluationStencil ( const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , int depth , Stencil< Point3D< double > , 3 > stencil [8] ) const; + void SetCornerNormalEvaluationStencils( const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , int depth , Stencil< Point3D< double > , 3 > stencils[8][8] ) const; + void SetCornerNormalEvaluationStencil ( const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , int depth , Stencil< Point3D< double > , 5 > stencil [8] ) const; + void SetCornerNormalEvaluationStencils( const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , int depth , Stencil< Point3D< double > , 5 > stencils[8][8] ) const; + + static void UpdateCoarserSupportBounds( const TreeOctNode* node , int& startX , int& endX , int& startY , int& endY , int& startZ , int& endZ ); + + void UpdateConstraintsFromCoarser( const PointInfo& pointInfo , const typename TreeOctNode::Neighbors5& neighbors5 , const typename TreeOctNode::Neighbors5& pNeighbors5 , TreeOctNode* node , Pointer( Real ) constraints , ConstPointer( Real ) metSolution , const typename BSplineData< Degree >::Integrator& integrator , const Stencil< double , 5 >& stencil ) const; + // Updates the constraints @(depth-1) based on the solution coefficients @(depth) + void UpdateConstraintsFromFiner( const typename BSplineData< Degree >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) fineSolution , Pointer( Real ) coarseConstraints ) const; + // Evaluate the points @(depth) using coefficients @(depth-1) + void SetPointValuesFromCoarser( PointInfo& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) coarseCoefficients ); + // Evalutes the solution @(depth) at the points @(depth-1) and updates the met constraints @(depth-1) + void SetPointConstraintsFromFiner( const PointInfo& pointInfo , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) finerCoefficients , Pointer( Real ) metConstraints) const; + Real _WeightedCoarserFunctionValue( const _PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) coarseCoefficients ) const; + Real _WeightedFinerFunctionValue ( const _PointData& pointData , const typename TreeOctNode::NeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) finerCoefficients ) const; + + // Down samples constraints @(depth) to constraints @(depth-1) + template< class C > void DownSample( int depth , const SortedTreeNodes& sNodes , ConstPointer( C ) fineConstraints , Pointer( C ) coarseConstraints ) const; + // Up samples solution @(depth-1) to solution @(depth) + template< class C > void UpSample ( int depth , const SortedTreeNodes& sNodes , ConstPointer( C ) coarseCoefficients , Pointer( C ) fineCoefficients ) const; + int GetSliceMatrixAndUpdateConstraints( const PointInfo& pointInfo , SparseMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< Degree >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine , int nStart , int nEnd ); + int GetMatrixAndUpdateConstraints( const PointInfo& pointInfo , SparseSymmetricMatrix< Real >& matrix , Pointer( Real ) constraints , const typename BSplineData< Degree >::Integrator& integrator , int depth , const SortedTreeNodes& sNodes , ConstPointer( Real ) metSolution , bool coarseToFine ); + + + int UpdateWeightContribution( std::vector< Real >& kernelDensityWeights , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , Real weight=Real(1.0) ); + Real GetSampleWeight( ConstPointer( Real ) kernelDensityWeight , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , int splatDepth ); + Real GetSampleWeight( ConstPointer( Real ) kernelDensityWeight , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey ); + void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , const TreeOctNode* node , const Point3D& position , typename TreeOctNode::ConstNeighborKey3& neighborKey , Real samplesPerNode , Real& depth , Real& weight ); + Real GetSampleWeight( ConstPointer( Real ) kernelDensityWeight , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey ); + void GetSampleDepthAndWeight( ConstPointer( Real ) kernelDensityWeight , TreeOctNode* node , const Point3D& position , typename TreeOctNode::NeighborKey3& neighborKey , Real samplesPerNode , Real& depth , Real& weight ); + int SplatOrientedPoint( ConstPointer( Real ) kernelDensityWeights , TreeOctNode* node , const Point3D& point , const Point3D< Real >& normal , NormalInfo& normalInfo , typename TreeOctNode::NeighborKey3& neighborKey ); + Real SplatOrientedPoint( ConstPointer( Real ) kernelDensityWeights , const Point3D& point , const Point3D& normal , NormalInfo& normalInfo , typename TreeOctNode::NeighborKey3& neighborKey , int kernelDepth , Real samplesPerNode , int minDepth , int maxDepth ); + + int HasNormals( TreeOctNode* node , const NormalInfo& normalInfo ); + + /////////////////////////// + // Iso-Surfacing Methods // + /////////////////////////// + struct IsoEdge + { + long long edges[2]; + IsoEdge( void ){ edges[0] = edges[1] = 0; } + IsoEdge( long long v1 , long long v2 ){ edges[0] = v1 , edges[1] = v2; } + long long& operator[]( int idx ){ return edges[idx]; } + const long long& operator[]( int idx ) const { return edges[idx]; } + }; + struct FaceEdges + { + IsoEdge edges[2]; + int count; + }; + template< class Vertex > + struct SliceValues + { + typename SortedTreeNodes::SliceTableData sliceData; + Pointer( Real ) cornerValues ; Pointer( Point3D< Real > ) cornerNormals ; Pointer( char ) cornerSet; + Pointer( long long ) edgeKeys ; Pointer( char ) edgeSet; + Pointer( FaceEdges ) faceEdges ; Pointer( char ) faceSet; + Pointer( char ) mcIndices; + hash_map< long long , std::vector< IsoEdge > > faceEdgeMap; + hash_map< long long , std::pair< int , Vertex > > edgeVertexMap; + hash_map< long long , long long > vertexPairMap; + + SliceValues( void ); + ~SliceValues( void ); + void reset( bool nonLinearFit ); + protected: + int _oldCCount , _oldECount , _oldFCount , _oldNCount; + }; + template< class Vertex > + struct XSliceValues + { + typename SortedTreeNodes::XSliceTableData xSliceData; + Pointer( long long ) edgeKeys ; Pointer( char ) edgeSet; + Pointer( FaceEdges ) faceEdges ; Pointer( char ) faceSet; + hash_map< long long , std::vector< IsoEdge > > faceEdgeMap; + hash_map< long long , std::pair< int , Vertex > > edgeVertexMap; + hash_map< long long , long long > vertexPairMap; + + XSliceValues( void ); + ~XSliceValues( void ); + void reset( void ); + protected: + int _oldECount , _oldFCount; + }; + template< class Vertex > + struct SlabValues + { + XSliceValues< Vertex > _xSliceValues[2]; + SliceValues< Vertex > _sliceValues[2]; + SliceValues< Vertex >& sliceValues( int idx ){ return _sliceValues[idx&1]; } + const SliceValues< Vertex >& sliceValues( int idx ) const { return _sliceValues[idx&1]; } + XSliceValues< Vertex >& xSliceValues( int idx ){ return _xSliceValues[idx&1]; } + const XSliceValues< Vertex >& xSliceValues( int idx ) const { return _xSliceValues[idx&1]; } + }; + template< class Vertex > + void SetSliceIsoCorners( ConstPointer( Real ) solution , ConstPointer( Real ) coarseSolution , Real isoValue , int depth , int slice , std::vector< SlabValues< Vertex > >& sValues , const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 > stencil[8] , const Stencil< double , 3 > stencils[8][8] , const Stencil< Point3D< double > , 3 > nStencil[8] , const Stencil< Point3D< double > , 3 > nStencils[8][8] , int threads ); + template< class Vertex > + void SetSliceIsoCorners( ConstPointer( Real ) solution , ConstPointer( Real ) coarseSolution , Real isoValue , int depth , int slice , int z , std::vector< SlabValues< Vertex > >& sValues , const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 > stencil[8] , const Stencil< double , 3 > stencils[8][8] , const Stencil< Point3D< double > , 3 > nStencil[8] , const Stencil< Point3D< double > , 3 > nStencils[8][8] , int threads ); + template< class Vertex > + void SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , Real isoValue , int depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); + template< class Vertex > + void SetSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , Real isoValue , int depth , int slice , int z , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); + template< class Vertex > + void SetXSliceIsoVertices( ConstPointer( Real ) kernelDensityWeights , Real isoValue , int depth , int slab , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< SlabValues< Vertex > >& sValues , int threads ); + template< class Vertex > + void CopyFinerSliceIsoEdgeKeys( int depth , int slice , std::vector< SlabValues< Vertex > >& sValues , int threads ); + template< class Vertex > + void CopyFinerSliceIsoEdgeKeys( int depth , int slice , int z , std::vector< SlabValues< Vertex > >& sValues , int threads ); + template< class Vertex > + void CopyFinerXSliceIsoEdgeKeys( int depth , int slab , std::vector< SlabValues< Vertex > >& sValues , int threads ); + template< class Vertex > + void SetSliceIsoEdges( int depth , int slice , std::vector< SlabValues< Vertex > >& slabValues , int threads ); + template< class Vertex > + void SetSliceIsoEdges( int depth , int slice , int z , std::vector< SlabValues< Vertex > >& slabValues , int threads ); + template< class Vertex > + void SetXSliceIsoEdges( int depth , int slice , std::vector< SlabValues< Vertex > >& slabValues , int threads ); + + template< class Vertex > + int SetIsoSurface( int depth , int offset , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , const XSliceValues< Vertex >& xValues , CoredMeshData< Vertex >& mesh , bool polygonMesh , bool addBarycenter , int& vOffset , int threads ); + + template< class Vertex > + static int AddIsoPolygons( CoredMeshData< Vertex >& mesh , std::vector< std::pair< int , Vertex > >& polygon , bool polygonMesh , bool addBarycenter , int& vOffset ); + + template< class Vertex > + bool GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int edgeIndex , int z , const SliceValues< Vertex >& sValues , Vertex& vertex ); + template< class Vertex > + bool GetIsoVertex( ConstPointer( Real ) kernelDensityWeights , Real isoValue , typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int cornerIndex , const SliceValues< Vertex >& bValues , const SliceValues< Vertex >& fValues , Vertex& vertex ); + + + //////////////////////// + // Evaluation Methods // + //////////////////////// + Real getCornerValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int corner , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 > stencils[8] , bool isInterior ) const; + Point3D< Real > getCornerNormal( const typename TreeOctNode::ConstNeighbors5& neighbors5 , const typename TreeOctNode::ConstNeighbors5& pNeighbors5 , const TreeOctNode* node , int corner , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , const Stencil< Point3D< double > , 5 >& nStencil , const Stencil< Point3D< double > , 5 > nStencils[8] , bool isInterior ) const; + std::pair< Real , Point3D< Real > > getCornerValueAndNormal( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , int corner , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< Degree >::template CornerEvaluator< 2 >& evaluator , const Stencil< double , 3 >& vStencil , const Stencil< double , 3 > vStencils[8] , const Stencil< Point3D< double > , 3 >& nStencil , const Stencil< Point3D< double > , 3 > nStencils[8] , bool isInterior ) const; + Real getCenterValue( const typename TreeOctNode::ConstNeighborKey3& neighborKey3 , const TreeOctNode* node , ConstPointer( Real ) solution , ConstPointer( Real ) metSolution , const typename BSplineData< Degree >::template CenterEvaluator< 1 >& evaluator , const Stencil< double , 3 >& stencil , const Stencil< double , 3 >& pStencil , bool isInterior ) const; + + static bool _IsInset( const TreeOctNode* node ); + static bool _IsInsetSupported( const TreeOctNode* node ); + + void refineBoundary( std::vector< int >* map ); +public: + int threads; + static double maxMemoryUsage; + TreeOctNode tree; + + static double MemoryUsage( void ); + Octree( void ); + + void MakeComplete( std::vector< int >* map=NULL ); + void Finalize( std::vector< int >* map=NULL ); + void ClipTree( const NormalInfo& normalInfo ); + + Real Evaluate( ConstPointer( Real ) coefficients , Point3D< Real > p , const BSplineData< Degree >* fData=NULL ) const; + Pointer( Real ) Evaluate( ConstPointer( Real ) coefficients , int& res , Real isoValue=0.f , int depth=-1 ); + template< class PointReal > + int SetTree( char* fileName , int minDepth , int maxDepth , int fullDepth , int splatDepth , Real samplesPerNode , + Real scaleFactor , bool useConfidence , bool useNormalWeight , Real constraintWeight , int adaptiveExponent , + PointInfo& pointInfo , NormalInfo& normalInfo , std::vector< Real >& kernelDensityWeights , std::vector< Real >& centerWeights , + int boundaryType=BSplineElements< Degree >::NONE , XForm4x4< Real > xForm=XForm4x4< Real >::Identity , bool makeComplete=false ); + Pointer( Real ) SetLaplacianConstraints( const NormalInfo& normalInfo ); + Pointer( Real ) SolveSystem( PointInfo& pointInfo , Pointer( Real ) constraints , bool showResidual , int iters , int maxSolveDepth , int cgDepth=0 , double cgAccuracy=0 ); + + Real GetIsoValue( ConstPointer( Real ) solution , const std::vector< Real >& centerWeights ); + template< class Vertex > + void GetMCIsoSurface( ConstPointer( Real ) kernelDensityWeights , ConstPointer( Real ) solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit=true , bool addBarycenter=false , bool polygonMesh=false ); +}; + +#include "MultiGridOctreeData.inl" +#include "MultiGridOctreeData.SortedTreeNodes.inl" +#include "MultiGridOctreeData.IsoSurface.inl" +#endif // MULTI_GRID_OCTREE_DATA_INCLUDED diff --git a/src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.h b/src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.h new file mode 100755 index 000000000..fc0651460 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.h @@ -0,0 +1,113 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#ifndef P_POLYNOMIAL_INCLUDED +#define P_POLYNOMIAL_INCLUDED +#include +#include "Polynomial.h" + +template +class StartingPolynomial{ +public: + Polynomial p; + double start; + + template + StartingPolynomial operator * (const StartingPolynomial& p) const; + StartingPolynomial scale(double s) const; + StartingPolynomial shift(double t) const; + int operator < (const StartingPolynomial& sp) const; + static int Compare(const void* v1,const void* v2); +}; + +template +class PPolynomial +{ +public: + size_t polyCount; + StartingPolynomial* polys; + + PPolynomial(void); + PPolynomial(const PPolynomial& p); + ~PPolynomial(void); + + PPolynomial& operator = (const PPolynomial& p); + + int size(void) const; + + void set( size_t size ); + // Note: this method will sort the elements in sps + void set( StartingPolynomial* sps , int count ); + void reset( size_t newSize ); + + + double operator()( double t ) const; + double integral( double tMin , double tMax ) const; + double Integral( void ) const; + + template + PPolynomial& operator = (const PPolynomial& p); + + PPolynomial operator + (const PPolynomial& p) const; + PPolynomial operator - (const PPolynomial& p) const; + + template + PPolynomial operator * (const Polynomial& p) const; + + template + PPolynomial operator * (const PPolynomial& p) const; + + + PPolynomial& operator += ( double s ); + PPolynomial& operator -= ( double s ); + PPolynomial& operator *= ( double s ); + PPolynomial& operator /= ( double s ); + PPolynomial operator + ( double s ) const; + PPolynomial operator - ( double s ) const; + PPolynomial operator * ( double s ) const; + PPolynomial operator / ( double s ) const; + + PPolynomial& addScaled(const PPolynomial& poly,double scale); + + PPolynomial scale( double s ) const; + PPolynomial shift( double t ) const; + + PPolynomial< Degree-1 > derivative(void) const; + PPolynomial< Degree+1 > integral(void) const; + + void getSolutions(double c,std::vector& roots,double EPS,double min=-DBL_MAX,double max=DBL_MAX) const; + + void printnl( void ) const; + + PPolynomial< Degree+1 > MovingAverage( double radius ) const; + static PPolynomial BSpline( double radius=0.5 ); + + void write( FILE* fp , int samples , double min , double max ) const; +}; +#include "PPolynomial.inl" +#endif // P_POLYNOMIAL_INCLUDED diff --git a/src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.inl b/src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.inl new file mode 100755 index 000000000..24b53fc83 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/PPolynomial.inl @@ -0,0 +1,431 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#include "Factor.h" + +//////////////////////// +// StartingPolynomial // +//////////////////////// +template +template +StartingPolynomial StartingPolynomial::operator * (const StartingPolynomial& p) const{ + StartingPolynomial sp; + if(start>p.start){sp.start=start;} + else{sp.start=p.start;} + sp.p=this->p*p.p; + return sp; +} +template +StartingPolynomial StartingPolynomial::scale(double s) const{ + StartingPolynomial q; + q.start=start*s; + q.p=p.scale(s); + return q; +} +template +StartingPolynomial StartingPolynomial::shift(double s) const{ + StartingPolynomial q; + q.start=start+s; + q.p=p.shift(s); + return q; +} + + +template +int StartingPolynomial::operator < (const StartingPolynomial& sp) const{ + if(start +int StartingPolynomial::Compare(const void* v1,const void* v2){ + double d=((StartingPolynomial*)(v1))->start-((StartingPolynomial*)(v2))->start; + if (d<0) {return -1;} + else if (d>0) {return 1;} + else {return 0;} +} + +///////////////// +// PPolynomial // +///////////////// +template +PPolynomial::PPolynomial(void){ + polyCount=0; + polys=NULL; +} +template +PPolynomial::PPolynomial(const PPolynomial& p){ + polyCount=0; + polys=NULL; + set(p.polyCount); + memcpy(polys,p.polys,sizeof(StartingPolynomial)*p.polyCount); +} + +template +PPolynomial::~PPolynomial(void){ + if(polyCount){free(polys);} + polyCount=0; + polys=NULL; +} +template +void PPolynomial::set(StartingPolynomial* sps,int count){ + int i,c=0; + set(count); + qsort(sps,count,sizeof(StartingPolynomial),StartingPolynomial::Compare); + for( i=0 ; i +int PPolynomial::size(void) const{return int(sizeof(StartingPolynomial)*polyCount);} + +template +void PPolynomial::set( size_t size ) +{ + if(polyCount){free(polys);} + polyCount=0; + polys=NULL; + polyCount=size; + if(size){ + polys=(StartingPolynomial*)malloc(sizeof(StartingPolynomial)*size); + memset(polys,0,sizeof(StartingPolynomial)*size); + } +} +template +void PPolynomial::reset( size_t newSize ) +{ + polyCount=newSize; + polys=(StartingPolynomial*)realloc(polys,sizeof(StartingPolynomial)*newSize); +} + +template +PPolynomial& PPolynomial::operator = (const PPolynomial& p){ + set(p.polyCount); + memcpy(polys,p.polys,sizeof(StartingPolynomial)*p.polyCount); + return *this; +} + +template +template +PPolynomial& PPolynomial::operator = (const PPolynomial& p){ + set(p.polyCount); + for(int i=0;i +double PPolynomial::operator ()( double t ) const +{ + double v=0; + for( int i=0 ; ipolys[i].start ; i++ ) v+=polys[i].p(t); + return v; +} + +template +double PPolynomial::integral( double tMin , double tMax ) const +{ + int m=1; + double start,end,s,v=0; + start=tMin; + end=tMax; + if(tMin>tMax){ + m=-1; + start=tMax; + end=tMin; + } + for(int i=0;i +double PPolynomial::Integral(void) const{return integral(polys[0].start,polys[polyCount-1].start);} +template +PPolynomial PPolynomial::operator + (const PPolynomial& p) const{ + PPolynomial q; + int i,j; + size_t idx=0; + q.set(polyCount+p.polyCount); + i=j=-1; + + while(idx=int(p.polyCount)-1) {q.polys[idx]= polys[++i];} + else if (i>=int( polyCount)-1) {q.polys[idx]=p.polys[++j];} + else if(polys[i+1].start +PPolynomial PPolynomial::operator - (const PPolynomial& p) const{ + PPolynomial q; + int i,j; + size_t idx=0; + q.set(polyCount+p.polyCount); + i=j=-1; + + while(idx=int(p.polyCount)-1) {q.polys[idx]= polys[++i];} + else if (i>=int( polyCount)-1) {q.polys[idx].start=p.polys[++j].start;q.polys[idx].p=p.polys[j].p*(-1.0);} + else if(polys[i+1].start +PPolynomial& PPolynomial::addScaled(const PPolynomial& p,double scale){ + int i,j; + StartingPolynomial* oldPolys=polys; + size_t idx=0,cnt=0,oldPolyCount=polyCount; + polyCount=0; + polys=NULL; + set(oldPolyCount+p.polyCount); + i=j=-1; + while(cnt=int( p.polyCount)-1) {polys[idx]=oldPolys[++i];} + else if (i>=int(oldPolyCount)-1) {polys[idx].start= p.polys[++j].start;polys[idx].p=p.polys[j].p*scale;} + else if (oldPolys[i+1].start +template +PPolynomial PPolynomial::operator * (const PPolynomial& p) const{ + PPolynomial q; + StartingPolynomial *sp; + int i,j,spCount=int(polyCount*p.polyCount); + + sp=(StartingPolynomial*)malloc(sizeof(StartingPolynomial)*spCount); + for(i=0;i +template +PPolynomial PPolynomial::operator * (const Polynomial& p) const{ + PPolynomial q; + q.set(polyCount); + for(int i=0;i +PPolynomial PPolynomial::scale( double s ) const +{ + PPolynomial q; + q.set(polyCount); + for(size_t i=0;i +PPolynomial PPolynomial::shift( double s ) const +{ + PPolynomial q; + q.set(polyCount); + for(size_t i=0;i +PPolynomial PPolynomial::derivative(void) const{ + PPolynomial q; + q.set(polyCount); + for(size_t i=0;i +PPolynomial PPolynomial::integral(void) const{ + int i; + PPolynomial q; + q.set(polyCount); + for(i=0;i +PPolynomial& PPolynomial::operator += ( double s ) {polys[0].p+=s;} +template +PPolynomial& PPolynomial::operator -= ( double s ) {polys[0].p-=s;} +template +PPolynomial& PPolynomial::operator *= ( double s ) +{ + for(int i=0;i +PPolynomial& PPolynomial::operator /= ( double s ) +{ + for(size_t i=0;i +PPolynomial PPolynomial::operator + ( double s ) const +{ + PPolynomial q=*this; + q+=s; + return q; +} +template +PPolynomial PPolynomial::operator - ( double s ) const +{ + PPolynomial q=*this; + q-=s; + return q; +} +template +PPolynomial PPolynomial::operator * ( double s ) const +{ + PPolynomial q=*this; + q*=s; + return q; +} +template +PPolynomial PPolynomial::operator / ( double s ) const +{ + PPolynomial q=*this; + q/=s; + return q; +} + +template +void PPolynomial::printnl(void) const{ + Polynomial p; + + if(!polyCount){ + Polynomial p; + printf("[-Infinity,Infinity]\n"); + } + else{ + for(size_t i=0;i +PPolynomial< 0 > PPolynomial< 0 >::BSpline( double radius ) +{ + PPolynomial q; + q.set(2); + + q.polys[0].start=-radius; + q.polys[1].start= radius; + + q.polys[0].p.coefficients[0]= 1.0; + q.polys[1].p.coefficients[0]=-1.0; + return q; +} +template< int Degree > +PPolynomial< Degree > PPolynomial::BSpline( double radius ) +{ + return PPolynomial< Degree-1 >::BSpline().MovingAverage( radius ); +} +template +PPolynomial PPolynomial::MovingAverage( double radius ) const +{ + PPolynomial A; + Polynomial p; + StartingPolynomial* sps; + + sps=(StartingPolynomial*)malloc(sizeof(StartingPolynomial)*polyCount*2); + + for(int i=0;i +void PPolynomial::getSolutions(double c,std::vector& roots,double EPS,double min,double max) const{ + Polynomial p; + std::vector tempRoots; + + p.setZero(); + for(size_t i=0;imax){break;} + if(ipolys[i].start && (i+1==polyCount || tempRoots[j]<=polys[i+1].start)){ + if(tempRoots[j]>min && tempRoots[j] +void PPolynomial::write(FILE* fp,int samples,double min,double max) const{ + fwrite(&samples,sizeof(int),1,fp); + for(int i=0;i +#include +#include +#include + +#define PLY_ASCII 1 /* ascii PLY file */ +#define PLY_BINARY_BE 2 /* binary PLY file, big endian */ +#define PLY_BINARY_LE 3 /* binary PLY file, little endian */ +#define PLY_BINARY_NATIVE 4 /* binary PLY file, same endianness as current architecture */ + +#define PLY_OKAY 0 /* ply routine worked okay */ +#define PLY_ERROR -1 /* error in ply routine */ + + /* scalar data types supported by PLY format */ + +#define PLY_START_TYPE 0 +#define PLY_CHAR 1 +#define PLY_SHORT 2 +#define PLY_INT 3 +#define PLY_UCHAR 4 +#define PLY_USHORT 5 +#define PLY_UINT 6 +#define PLY_FLOAT 7 +#define PLY_DOUBLE 8 +#define PLY_INT_8 9 +#define PLY_UINT_8 10 +#define PLY_INT_16 11 +#define PLY_UINT_16 12 +#define PLY_INT_32 13 +#define PLY_UINT_32 14 +#define PLY_FLOAT_32 15 +#define PLY_FLOAT_64 16 + +#define PLY_END_TYPE 17 + +#define PLY_SCALAR 0 +#define PLY_LIST 1 + + +typedef struct PlyProperty { /* description of a property */ + + char *name; /* property name */ + int external_type; /* file's data type */ + int internal_type; /* program's data type */ + int offset; /* offset bytes of prop in a struct */ + + int is_list; /* 1 = list, 0 = scalar */ + int count_external; /* file's count type */ + int count_internal; /* program's count type */ + int count_offset; /* offset byte for list count */ + +} PlyProperty; + +typedef struct PlyElement { /* description of an element */ + char *name; /* element name */ + int num; /* number of elements in this object */ + int size; /* size of element (bytes) or -1 if variable */ + int nprops; /* number of properties for this element */ + PlyProperty **props; /* list of properties in the file */ + char *store_prop; /* flags: property wanted by user? */ + int other_offset; /* offset to un-asked-for props, or -1 if none*/ + int other_size; /* size of other_props structure */ +} PlyElement; + +typedef struct PlyOtherProp { /* describes other properties in an element */ + char *name; /* element name */ + int size; /* size of other_props */ + int nprops; /* number of properties in other_props */ + PlyProperty **props; /* list of properties in other_props */ +} PlyOtherProp; + +typedef struct OtherData { /* for storing other_props for an other element */ + void *other_props; +} OtherData; + +typedef struct OtherElem { /* data for one "other" element */ + char *elem_name; /* names of other elements */ + int elem_count; /* count of instances of each element */ + OtherData **other_data; /* actual property data for the elements */ + PlyOtherProp *other_props; /* description of the property data */ +} OtherElem; + +typedef struct PlyOtherElems { /* "other" elements, not interpreted by user */ + int num_elems; /* number of other elements */ + OtherElem *other_list; /* list of data for other elements */ +} PlyOtherElems; + +typedef struct PlyFile { /* description of PLY file */ + FILE *fp; /* file pointer */ + int file_type; /* ascii or binary */ + float version; /* version number of file */ + int nelems; /* number of elements of object */ + PlyElement **elems; /* list of elements */ + int num_comments; /* number of comments */ + char **comments; /* list of comments */ + int num_obj_info; /* number of items of object information */ + char **obj_info; /* list of object info items */ + PlyElement *which_elem; /* which element we're currently writing */ + PlyOtherElems *other_elems; /* "other" elements from a PLY file */ +} PlyFile; + + /* memory allocation */ +extern char *my_alloc(); +#define myalloc(mem_size) my_alloc((mem_size), __LINE__, __FILE__) + +#ifndef ALLOCN +#define REALLOCN(PTR,TYPE,OLD_N,NEW_N) \ +{ \ + if ((OLD_N) == 0) \ +{ ALLOCN((PTR),TYPE,(NEW_N));} \ + else \ +{ \ + (PTR) = (TYPE *)realloc((PTR),(NEW_N)*sizeof(TYPE)); \ + if (((PTR) == NULL) && ((NEW_N) != 0)) \ +{ \ + fprintf(stderr, "Memory reallocation failed on line %d in %s\n", \ + __LINE__, __FILE__); \ + fprintf(stderr, " tried to reallocate %d->%d\n", \ + (OLD_N), (NEW_N)); \ + exit(-1); \ +} \ + if ((NEW_N)>(OLD_N)) \ + memset((char *)(PTR)+(OLD_N)*sizeof(TYPE), 0, \ + ((NEW_N)-(OLD_N))*sizeof(TYPE)); \ +} \ +} + +#define ALLOCN(PTR,TYPE,N) \ +{ (PTR) = (TYPE *) calloc(((unsigned)(N)),sizeof(TYPE));\ + if ((PTR) == NULL) { \ + fprintf(stderr, "Memory allocation failed on line %d in %s\n", \ + __LINE__, __FILE__); \ + exit(-1); \ + } \ +} + + +#define FREE(PTR) { free((PTR)); (PTR) = NULL; } +#endif + + +/*** delcaration of routines ***/ + +extern PlyFile *ply_write(FILE *, int, const char **, int); +extern PlyFile *ply_open_for_writing(char *, int, const char **, int, float *); +extern void ply_describe_element(PlyFile *, char *, int, int, PlyProperty *); +extern void ply_describe_property(PlyFile *, char *, PlyProperty *); +extern void ply_element_count(PlyFile *, char *, int); +extern void ply_header_complete(PlyFile *); +extern void ply_put_element_setup(PlyFile *, char *); +extern void ply_put_element(PlyFile *, void *); +extern void ply_put_comment(PlyFile *, char *); +extern void ply_put_obj_info(PlyFile *, char *); +extern PlyFile *ply_read(FILE *, int *, char ***); +extern PlyFile *ply_open_for_reading( char *, int *, char ***, int *, float *); +extern PlyProperty **ply_get_element_description(PlyFile *, char *, int*, int*); +extern void ply_get_element_setup( PlyFile *, char *, int, PlyProperty *); +extern int ply_get_property(PlyFile *, char *, PlyProperty *); +extern PlyOtherProp *ply_get_other_properties(PlyFile *, char *, int); +extern void ply_get_element(PlyFile *, void *); +extern char **ply_get_comments(PlyFile *, int *); +extern char **ply_get_obj_info(PlyFile *, int *); +extern void ply_close(PlyFile *); +extern void ply_get_info(PlyFile *, float *, int *); +extern PlyOtherElems *ply_get_other_element (PlyFile *, char *, int); +extern void ply_describe_other_elements ( PlyFile *, PlyOtherElems *); +extern void ply_put_other_elements (PlyFile *); +extern void ply_free_other_elements (PlyOtherElems *); +extern void ply_describe_other_properties(PlyFile *, PlyOtherProp *, int); + +extern int equal_strings(const char *, const char *); + +#ifdef __cplusplus +} +#endif +#include "Geometry.h" +#include + +typedef struct PlyFace +{ + unsigned char nr_vertices; + int *vertices; + int segment; +} PlyFace; +static PlyProperty face_props[] = +{ + { "vertex_indices" , PLY_INT , PLY_INT , offsetof(PlyFace,vertices) , 1 , PLY_UCHAR, PLY_UCHAR , offsetof(PlyFace,nr_vertices) }, +}; + +template< class Real > +class PlyVertex +{ +public: + const static int Components=3; + static PlyProperty Properties[]; + + Point3D< Real > point; + + PlyVertex( void ) { ; } + PlyVertex( Point3D< Real > p ) { point=p; } + PlyVertex operator + ( PlyVertex p ) const { return PlyVertex( point+p.point ); } + PlyVertex operator - ( PlyVertex p ) const { return PlyVertex( point-p.point ); } + template< class _Real > PlyVertex operator * ( _Real s ) const { return PlyVertex( point*s ); } + template< class _Real > PlyVertex operator / ( _Real s ) const { return PlyVertex( point/s ); } + PlyVertex& operator += ( PlyVertex p ) { point += p.point ; return *this; } + PlyVertex& operator -= ( PlyVertex p ) { point -= p.point ; return *this; } + template< class _Real > PlyVertex& operator *= ( _Real s ) { point *= s ; return *this; } + template< class _Real > PlyVertex& operator /= ( _Real s ) { point /= s ; return *this; } +}; +template< class Real , class _Real > PlyVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyVertex< Real > v ) { return PlyVertex< Real >( xForm * v.point ); } +template<> +PlyProperty PlyVertex< float >::Properties[]= +{ + {"x", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyVertex,point.coords[2])), 0, 0, 0, 0} +}; +template<> +PlyProperty PlyVertex< double >::Properties[]= +{ + {"x", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyVertex,point.coords[2])), 0, 0, 0, 0} +}; +template< class Real > +class PlyValueVertex +{ +public: + const static int Components=4; + static PlyProperty Properties[]; + + Point3D point; + Real value; + + PlyValueVertex( void ) : value( Real(0) ) { ; } + PlyValueVertex( Point3D< Real > p , Real v ) : point(p) , value(v) { ; } + PlyValueVertex operator + ( PlyValueVertex p ) const { return PlyValueVertex( point+p.point , value+p.value ); } + PlyValueVertex operator - ( PlyValueVertex p ) const { return PlyValueVertex( point-p.value , value-p.value ); } + template< class _Real > PlyValueVertex operator * ( _Real s ) const { return PlyValueVertex( point*s , Real(value*s) ); } + template< class _Real > PlyValueVertex operator / ( _Real s ) const { return PlyValueVertex( point/s , Real(value/s) ); } + PlyValueVertex& operator += ( PlyValueVertex p ) { point += p.point , value += p.value ; return *this; } + PlyValueVertex& operator -= ( PlyValueVertex p ) { point -= p.point , value -= p.value ; return *this; } + template< class _Real > PlyValueVertex& operator *= ( _Real s ) { point *= s , value *= Real(s) ; return *this; } + template< class _Real > PlyValueVertex& operator /= ( _Real s ) { point /= s , value /= Real(s) ; return *this; } +}; +template< class Real , class _Real > PlyValueVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyValueVertex< Real > v ) { return PlyValueVertex< Real >( xForm * v.point , v.value ); } +template< > +PlyProperty PlyValueVertex< float >::Properties[]= +{ + {"x", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyValueVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyValueVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyValueVertex,point.coords[2])), 0, 0, 0, 0}, + {"value", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyValueVertex,value)), 0, 0, 0, 0} +}; +template< > +PlyProperty PlyValueVertex< double >::Properties[]= +{ + {"x", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyValueVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyValueVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyValueVertex,point.coords[2])), 0, 0, 0, 0}, + {"value", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyValueVertex,value)), 0, 0, 0, 0} +}; +template< class Real > +class PlyOrientedVertex +{ +public: + const static int Components=6; + static PlyProperty Properties[]; + + Point3D point , normal; + + PlyOrientedVertex( void ) { ; } + PlyOrientedVertex( Point3D< Real > p , Point3D< Real > n ) : point(p) , normal(n) { ; } + PlyOrientedVertex operator + ( PlyOrientedVertex p ) const { return PlyOrientedVertex( point+p.point , normal+p.normal ); } + PlyOrientedVertex operator - ( PlyOrientedVertex p ) const { return PlyOrientedVertex( point-p.value , normal-p.normal ); } + template< class _Real > PlyOrientedVertex operator * ( _Real s ) const { return PlyOrientedVertex( point*s , normal*s ); } + template< class _Real > PlyOrientedVertex operator / ( _Real s ) const { return PlyOrientedVertex( point/s , normal/s ); } + PlyOrientedVertex& operator += ( PlyOrientedVertex p ) { point += p.point , normal += p.normal ; return *this; } + PlyOrientedVertex& operator -= ( PlyOrientedVertex p ) { point -= p.point , normal -= p.normal ; return *this; } + template< class _Real > PlyOrientedVertex& operator *= ( _Real s ) { point *= s , normal *= s ; return *this; } + template< class _Real > PlyOrientedVertex& operator /= ( _Real s ) { point /= s , normal /= s ; return *this; } +}; +template< class Real , class _Real > PlyOrientedVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyOrientedVertex< Real > v ) { return PlyOrientedVertex< Real >( xForm * v.point , xForm.inverse().transpose() * v.normal ); } +template<> +PlyProperty PlyOrientedVertex< float >::Properties[]= +{ + {"x", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyOrientedVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyOrientedVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyOrientedVertex,point.coords[2])), 0, 0, 0, 0}, + {"nx", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyOrientedVertex,normal.coords[0])), 0, 0, 0, 0}, + {"ny", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyOrientedVertex,normal.coords[1])), 0, 0, 0, 0}, + {"nz", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyOrientedVertex,normal.coords[2])), 0, 0, 0, 0} +}; +template<> +PlyProperty PlyOrientedVertex< double >::Properties[]= +{ + {"x", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyOrientedVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyOrientedVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyOrientedVertex,point.coords[2])), 0, 0, 0, 0}, + {"nx", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyOrientedVertex,normal.coords[0])), 0, 0, 0, 0}, + {"ny", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyOrientedVertex,normal.coords[1])), 0, 0, 0, 0}, + {"nz", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyOrientedVertex,normal.coords[2])), 0, 0, 0, 0} +}; +template< class Real > +class PlyColorVertex +{ +public: + const static int Components=6; + static PlyProperty Properties[]; + + Point3D point; + unsigned char color[3]; + + operator Point3D& () {return point;} + operator const Point3D& () const {return point;} + PlyColorVertex(void) {point.coords[0]=point.coords[1]=point.coords[2]=0,color[0]=color[1]=color[2]=0;} + PlyColorVertex(const Point3D& p) {point=p;} +}; +template<> +PlyProperty PlyColorVertex< float >::Properties[]= +{ + {"x", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyColorVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyColorVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_FLOAT, PLY_FLOAT, int(offsetof(PlyColorVertex,point.coords[2])), 0, 0, 0, 0}, + {"red", PLY_UCHAR, PLY_UCHAR, int(offsetof(PlyColorVertex,color[0])), 0, 0, 0, 0}, + {"green", PLY_UCHAR, PLY_UCHAR, int(offsetof(PlyColorVertex,color[1])), 0, 0, 0, 0}, + {"blue", PLY_UCHAR, PLY_UCHAR, int(offsetof(PlyColorVertex,color[2])), 0, 0, 0, 0} +}; +template<> +PlyProperty PlyColorVertex< double >::Properties[]= +{ + {"x", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyColorVertex,point.coords[0])), 0, 0, 0, 0}, + {"y", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyColorVertex,point.coords[1])), 0, 0, 0, 0}, + {"z", PLY_DOUBLE, PLY_DOUBLE, int(offsetof(PlyColorVertex,point.coords[2])), 0, 0, 0, 0}, + {"red", PLY_UCHAR, PLY_UCHAR, int(offsetof(PlyColorVertex,color[0])), 0, 0, 0, 0}, + {"green", PLY_UCHAR, PLY_UCHAR, int(offsetof(PlyColorVertex,color[1])), 0, 0, 0, 0}, + {"blue", PLY_UCHAR, PLY_UCHAR, int(offsetof(PlyColorVertex,color[2])), 0, 0, 0, 0} +}; + +template< class Vertex , class Real > +int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , const Point3D< float >& translate , float scale , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() ); + +template< class Vertex , class Real > +int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() ); + +template +int PlyReadPolygons(char* fileName, + std::vector& vertices,std::vector >& polygons, + PlyProperty* properties,int propertyNum, + int& file_type, + char*** comments=NULL,int* commentNum=NULL , bool* readFlags=NULL ); + +template +int PlyWritePolygons(char* fileName, + const std::vector& vertices,const std::vector >& polygons, + PlyProperty* properties,int propertyNum, + int file_type, + char** comments=NULL,const int& commentNum=0); + +template +int PlyWritePolygons(char* fileName, + const std::vector& vertices , const std::vector< std::vector< int > >& polygons, + PlyProperty* properties,int propertyNum, + int file_type, + char** comments,const int& commentNum) +{ + int nr_vertices=int(vertices.size()); + int nr_faces=int(polygons.size()); + float version; + const char *elem_names[] = { "vertex" , "face" }; + PlyFile *ply = ply_open_for_writing( fileName , 2 , elem_names , file_type , &version ); + if (!ply){return 0;} + + // + // describe vertex and face properties + // + ply_element_count(ply, "vertex", nr_vertices); + for(int i=0;imaxFaceVerts) + { + delete[] ply_face.vertices; + maxFaceVerts=int(polygons[i].size()); + ply_face.vertices=new int[maxFaceVerts]; + } + ply_face.nr_vertices=int(polygons[i].size()); + for(int j=0;j +int PlyReadPolygons(char* fileName, + std::vector& vertices , std::vector >& polygons , + PlyProperty* properties , int propertyNum , + int& file_type , + char*** comments , int* commentNum , bool* readFlags ) +{ + int nr_elems; + char **elist; + float version; + int i,j,k; + PlyFile* ply; + char* elem_name; + int num_elems; + int nr_props; + PlyProperty** plist; + PlyFace ply_face; + + ply = ply_open_for_reading(fileName, &nr_elems, &elist, &file_type, &version); + if(!ply) return 0; + + if(comments) + { + (*comments)=new char*[*commentNum+ply->num_comments]; + for(int i=0;inum_comments;i++) + (*comments)[i]=_strdup(ply->comments[i]); + *commentNum=ply->num_comments; + } + + for (i=0; i < nr_elems; i++) { + elem_name = elist[i]; + plist = ply_get_element_description(ply, elem_name, &num_elems, &nr_props); + if(!plist) + { + for(i=0;ielems[i]->name); + free(ply->elems[i]->store_prop); + for(j=0;jelems[i]->nprops;j++){ + free(ply->elems[i]->props[j]->name); + free(ply->elems[i]->props[j]); + } + free(ply->elems[i]->props); + } + for(i=0;ielems[i]);} + free(ply->elems); + for(i=0;inum_comments;i++){free(ply->comments[i]);} + free(ply->comments); + for(i=0;inum_obj_info;i++){free(ply->obj_info[i]);} + free(ply->obj_info); + ply_free_other_elements (ply->other_elems); + + for(i=0;iname); + free(plist[j]); + } + free(plist); + } // for each type of element + + for(i=0;ielems[i]->name); + free(ply->elems[i]->store_prop); + for(j=0;jelems[i]->nprops;j++){ + free(ply->elems[i]->props[j]->name); + free(ply->elems[i]->props[j]); + } + if(ply->elems[i]->props && ply->elems[i]->nprops){free(ply->elems[i]->props);} + } + for(i=0;ielems[i]);} + free(ply->elems); + for(i=0;inum_comments;i++){free(ply->comments[i]);} + free(ply->comments); + for(i=0;inum_obj_info;i++){free(ply->obj_info[i]);} + free(ply->obj_info); + ply_free_other_elements (ply->other_elems); + + + for(i=0;i +int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , const Point3D& translate , float scale , char** comments , int commentNum , XForm4x4< Real > xForm ) +{ + int i; + int nr_vertices=int(mesh->outOfCorePointCount()+mesh->inCorePoints.size()); + int nr_faces=mesh->polygonCount(); + float version; + const char *elem_names[] = { "vertex" , "face" }; + PlyFile *ply = ply_open_for_writing( fileName , 2 , elem_names , file_type , &version ); + if( !ply ) return 0; + + mesh->resetIterator(); + + // + // describe vertex and face properties + // + ply_element_count( ply , "vertex" , nr_vertices ); + for( int i=0 ; iinCorePoints.size() ) ; i++ ) + { + Vertex vertex = xForm * ( mesh->inCorePoints[i] * scale + translate ); + ply_put_element(ply, (void *) &vertex); + } + for( i=0; ioutOfCorePointCount() ; i++ ) + { + Vertex vertex; + mesh->nextOutOfCorePoint( vertex ); + vertex = xForm * ( vertex * scale +translate ); + ply_put_element(ply, (void *) &vertex); + } // for, write vertices + + // write faces + std::vector< CoredVertexIndex > polygon; + ply_put_element_setup( ply , "face" ); + for( i=0 ; inextPolygon( polygon ); + ply_face.nr_vertices = int( polygon.size() ); + ply_face.vertices = new int[ polygon.size() ]; + for( int i=0 ; iinCorePoints.size() ); + ply_put_element( ply, (void *) &ply_face ); + delete[] ply_face.vertices; + } // for, write faces + + ply_close( ply ); + return 1; +} +template< class Vertex , class Real > +int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , char** comments , int commentNum , XForm4x4< Real > xForm ) +{ + int i; + int nr_vertices=int(mesh->outOfCorePointCount()+mesh->inCorePoints.size()); + int nr_faces=mesh->polygonCount(); + float version; + const char *elem_names[] = { "vertex" , "face" }; + PlyFile *ply = ply_open_for_writing( fileName , 2 , elem_names , file_type , &version ); + if( !ply ) return 0; + + mesh->resetIterator(); + + // + // describe vertex and face properties + // + ply_element_count( ply , "vertex" , nr_vertices ); + for( int i=0 ; iinCorePoints.size() ) ; i++ ) + { + Vertex vertex = xForm * mesh->inCorePoints[i]; + ply_put_element(ply, (void *) &vertex); + } + for( i=0; ioutOfCorePointCount() ; i++ ) + { + Vertex vertex; + mesh->nextOutOfCorePoint( vertex ); + vertex = xForm * ( vertex ); + ply_put_element(ply, (void *) &vertex); + } // for, write vertices + + // write faces + std::vector< CoredVertexIndex > polygon; + ply_put_element_setup( ply , "face" ); + for( i=0 ; inextPolygon( polygon ); + ply_face.nr_vertices = int( polygon.size() ); + ply_face.vertices = new int[ polygon.size() ]; + for( int i=0 ; iinCorePoints.size() ); + ply_put_element( ply, (void *) &ply_face ); + delete[] ply_face.vertices; + } // for, write faces + + ply_close( ply ); + return 1; +} +inline int PlyDefaultFileType(void){return PLY_ASCII;} + +#endif /* !__PLY_H__ */ diff --git a/src/plugins_experimental/filter_screened_poisson/Src/PlyFile.cpp b/src/plugins_experimental/filter_screened_poisson/Src/PlyFile.cpp new file mode 100755 index 000000000..d3e70f8e2 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/PlyFile.cpp @@ -0,0 +1,2727 @@ +/* + + The interface routines for reading and writing PLY polygon files. + + Greg Turk, February 1994 + + --------------------------------------------------------------- + + A PLY file contains a single polygonal _object_. + + An object is composed of lists of _elements_. Typical elements are + vertices, faces, edges and materials. + + Each type of element for a given object has one or more _properties_ + associated with the element type. For instance, a vertex element may + have as properties the floating-point values x,y,z and the three unsigned + chars representing red, green and blue. + + --------------------------------------------------------------- + + Copyright (c) 1994 The Board of Trustees of The Leland Stanford + Junior University. All rights reserved. + + Permission to use, copy, modify and distribute this software and its + documentation for any purpose is hereby granted without fee, provided + that the above copyright notice and this permission notice appear in + all copies of this software and that you do not sell the software. + + THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND, + EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY + WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. + +*/ + +#include +#include +#include +#include +#include "PlyFile.h" + +char *type_names[] = { + "invalid", + "char", + "short", + "int", + "uchar", + "ushort", + "uint", + "float", + "double", + + "int8", // character 1 + "uint8", // unsigned character 1 + "int16", // short integer 2 + "uint16", // unsigned short integer 2 + "int32", // integer 4 + "uint32", // unsigned integer 4 + "float32", // single-precision float 4 + "float64", // double-precision float 8 + +}; + +int ply_type_size[] = { + 0, + 1, + 2, + 4, + 1, + 2, + 4, + 4, + 8, + 1, + 1, + 2, + 2, + 4, + 4, + 8 +}; + +typedef union +{ + int int_value; + char byte_values[sizeof(int)]; +} endian_test_type; + + +static int native_binary_type = -1; +static int types_checked = 0; + +#define NO_OTHER_PROPS -1 + +#define DONT_STORE_PROP 0 +#define STORE_PROP 1 + +#define OTHER_PROP 0 +#define NAMED_PROP 1 + + +/* returns 1 if strings are equal, 0 if not */ +int equal_strings(const char *, const char *); + +/* find an element in a plyfile's list */ +PlyElement *find_element(PlyFile *, char *); + +/* find a property in an element's list */ +PlyProperty *find_property(PlyElement *, char *, int *); + +/* write to a file the word describing a PLY file data type */ +void write_scalar_type (FILE *, int); + +/* read a line from a file and break it up into separate words */ +char **get_words(FILE *, int *, char **); +char **old_get_words(FILE *, int *); + +/* write an item to a file */ +void write_binary_item(FILE *, int, int, unsigned int, double, int); +void write_ascii_item(FILE *, int, unsigned int, double, int); +double old_write_ascii_item(FILE *, char *, int); + +/* add information to a PLY file descriptor */ +void add_element(PlyFile *, char **); +void add_property(PlyFile *, char **); +void add_comment(PlyFile *, char *); +void add_obj_info(PlyFile *, char *); + +/* copy a property */ +void copy_property(PlyProperty *, PlyProperty *); + +/* store a value into where a pointer and a type specify */ +void store_item(char *, int, int, unsigned int, double); + +/* return the value of a stored item */ +void get_stored_item( void *, int, int *, unsigned int *, double *); + +/* return the value stored in an item, given ptr to it and its type */ +double get_item_value(char *, int); + +/* get binary or ascii item and store it according to ptr and type */ +void get_ascii_item(char *, int, int *, unsigned int *, double *); +void get_binary_item(FILE *, int, int, int *, unsigned int *, double *); + +/* get a bunch of elements from a file */ +void ascii_get_element(PlyFile *, char *); +void binary_get_element(PlyFile *, char *); + +/* memory allocation */ +char *my_alloc(int, int, char *); + +/* byte ordering */ +void get_native_binary_type(); +void swap_bytes(char *, int); + +void check_types(); + +/*************/ +/* Writing */ +/*************/ + + +/****************************************************************************** +Given a file pointer, get ready to write PLY data to the file. + + Entry: + fp - the given file pointer + nelems - number of elements in object + elem_names - list of element names + file_type - file type, either ascii or binary + + Exit: + returns a pointer to a PlyFile, used to refer to this file, or NULL if error +******************************************************************************/ + +PlyFile *ply_write( + FILE *fp, + int nelems, + const char **elem_names, + int file_type + ) +{ + int i; + PlyFile *plyfile; + PlyElement *elem; + + /* check for NULL file pointer */ + if (fp == NULL) + return (NULL); + + if (native_binary_type == -1) + get_native_binary_type(); + if (!types_checked) + check_types(); + + /* create a record for this object */ + + plyfile = (PlyFile *) myalloc (sizeof (PlyFile)); + if (file_type == PLY_BINARY_NATIVE) + plyfile->file_type = native_binary_type; + else + plyfile->file_type = file_type; + plyfile->num_comments = 0; + plyfile->num_obj_info = 0; + plyfile->nelems = nelems; + plyfile->version = 1.0; + plyfile->fp = fp; + plyfile->other_elems = NULL; + + /* tuck aside the names of the elements */ + + plyfile->elems = (PlyElement **) myalloc (sizeof (PlyElement *) * nelems); + for (i = 0; i < nelems; i++) { + elem = (PlyElement *) myalloc (sizeof (PlyElement)); + plyfile->elems[i] = elem; + elem->name = _strdup (elem_names[i]); + elem->num = 0; + elem->nprops = 0; + } + + /* return pointer to the file descriptor */ + return (plyfile); +} + + +/****************************************************************************** +Open a polygon file for writing. + + Entry: + filename - name of file to read from + nelems - number of elements in object + elem_names - list of element names + file_type - file type, either ascii or binary + + Exit: + version - version number of PLY file + returns a file identifier, used to refer to this file, or NULL if error +******************************************************************************/ + +PlyFile *ply_open_for_writing( + char *filename, + int nelems, + const char **elem_names, + int file_type, + float *version + ) +{ + PlyFile *plyfile; + char *name; + FILE *fp; + + /* tack on the extension .ply, if necessary */ + + name = (char *) myalloc (int(sizeof (char) * (strlen (filename)) + 5)); + strcpy (name, filename); + if (strlen (name) < 4 || + strcmp (name + strlen (name) - 4, ".ply") != 0) + strcat (name, ".ply"); + + /* open the file for writing */ + + fp = fopen (name, "wb"); + if (fp == NULL) { + return (NULL); + } + + /* create the actual PlyFile structure */ + + plyfile = ply_write (fp, nelems, elem_names, file_type); + if (plyfile == NULL) + return (NULL); + + /* say what PLY file version number we're writing */ + *version = plyfile->version; + + /* return pointer to the file descriptor */ + return (plyfile); +} + + +/****************************************************************************** +Describe an element, including its properties and how many will be written +to the file. + + Entry: + plyfile - file identifier + elem_name - name of element that information is being specified about + nelems - number of elements of this type to be written + nprops - number of properties contained in the element + prop_list - list of properties +******************************************************************************/ + +void ply_describe_element( + PlyFile *plyfile, + char *elem_name, + int nelems, + int nprops, + PlyProperty *prop_list + ) +{ + int i; + PlyElement *elem; + PlyProperty *prop; + + /* look for appropriate element */ + elem = find_element (plyfile, elem_name); + if (elem == NULL) { + fprintf(stderr,"ply_describe_element: can't find element '%s'\n",elem_name); + exit (-1); + } + + elem->num = nelems; + + /* copy the list of properties */ + + elem->nprops = nprops; + elem->props = (PlyProperty **) myalloc (sizeof (PlyProperty *) * nprops); + elem->store_prop = (char *) myalloc (sizeof (char) * nprops); + + for (i = 0; i < nprops; i++) { + prop = (PlyProperty *) myalloc (sizeof (PlyProperty)); + elem->props[i] = prop; + elem->store_prop[i] = NAMED_PROP; + copy_property (prop, &prop_list[i]); + } +} + + +/****************************************************************************** +Describe a property of an element. + + Entry: + plyfile - file identifier + elem_name - name of element that information is being specified about + prop - the new property +******************************************************************************/ + +void ply_describe_property( + PlyFile *plyfile, + char *elem_name, + PlyProperty *prop + ) +{ + PlyElement *elem; + PlyProperty *elem_prop; + + /* look for appropriate element */ + elem = find_element (plyfile, elem_name); + if (elem == NULL) { + fprintf(stderr, "ply_describe_property: can't find element '%s'\n", + elem_name); + return; + } + + /* create room for new property */ + + if (elem->nprops == 0) { + elem->props = (PlyProperty **) myalloc (sizeof (PlyProperty *)); + elem->store_prop = (char *) myalloc (sizeof (char)); + elem->nprops = 1; + } + else { + elem->nprops++; + elem->props = (PlyProperty **) + realloc (elem->props, sizeof (PlyProperty *) * elem->nprops); + elem->store_prop = (char *) + realloc (elem->store_prop, sizeof (char) * elem->nprops); + } + + /* copy the new property */ + + elem_prop = (PlyProperty *) myalloc (sizeof (PlyProperty)); + elem->props[elem->nprops - 1] = elem_prop; + elem->store_prop[elem->nprops - 1] = NAMED_PROP; + copy_property (elem_prop, prop); +} + + +/****************************************************************************** +Describe what the "other" properties are that are to be stored, and where +they are in an element. +******************************************************************************/ + +void ply_describe_other_properties( + PlyFile *plyfile, + PlyOtherProp *other, + int offset + ) +{ + int i; + PlyElement *elem; + PlyProperty *prop; + + /* look for appropriate element */ + elem = find_element (plyfile, other->name); + if (elem == NULL) { + fprintf(stderr, "ply_describe_other_properties: can't find element '%s'\n", + other->name); + return; + } + + /* create room for other properties */ + + if (elem->nprops == 0) { + elem->props = (PlyProperty **) + myalloc (sizeof (PlyProperty *) * other->nprops); + elem->store_prop = (char *) myalloc (sizeof (char) * other->nprops); + elem->nprops = 0; + } + else { + int newsize; + newsize = elem->nprops + other->nprops; + elem->props = (PlyProperty **) + realloc (elem->props, sizeof (PlyProperty *) * newsize); + elem->store_prop = (char *) + realloc (elem->store_prop, sizeof (char) * newsize); + } + + /* copy the other properties */ + + for (i = 0; i < other->nprops; i++) { + prop = (PlyProperty *) myalloc (sizeof (PlyProperty)); + copy_property (prop, other->props[i]); + elem->props[elem->nprops] = prop; + elem->store_prop[elem->nprops] = OTHER_PROP; + elem->nprops++; + } + + /* save other info about other properties */ + elem->other_size = other->size; + elem->other_offset = offset; +} + + +/****************************************************************************** +State how many of a given element will be written. + + Entry: + plyfile - file identifier + elem_name - name of element that information is being specified about + nelems - number of elements of this type to be written +******************************************************************************/ + +void ply_element_count( + PlyFile *plyfile, + char *elem_name, + int nelems + ) +{ + PlyElement *elem; + + /* look for appropriate element */ + elem = find_element (plyfile, elem_name); + if (elem == NULL) { + fprintf(stderr,"ply_element_count: can't find element '%s'\n",elem_name); + exit (-1); + } + + elem->num = nelems; +} + + +/****************************************************************************** +Signal that we've described everything a PLY file's header and that the +header should be written to the file. + + Entry: + plyfile - file identifier +******************************************************************************/ + +void ply_header_complete(PlyFile *plyfile) +{ + int i,j; + FILE *fp = plyfile->fp; + PlyElement *elem; + PlyProperty *prop; + + fprintf (fp, "ply\n"); + + switch (plyfile->file_type) { + case PLY_ASCII: + fprintf (fp, "format ascii 1.0\n"); + break; + case PLY_BINARY_BE: + fprintf (fp, "format binary_big_endian 1.0\n"); + break; + case PLY_BINARY_LE: + fprintf (fp, "format binary_little_endian 1.0\n"); + break; + default: + fprintf (stderr, "ply_header_complete: bad file type = %d\n", + plyfile->file_type); + exit (-1); + } + + /* write out the comments */ + + for (i = 0; i < plyfile->num_comments; i++) + fprintf (fp, "comment %s\n", plyfile->comments[i]); + + /* write out object information */ + + for (i = 0; i < plyfile->num_obj_info; i++) + fprintf (fp, "obj_info %s\n", plyfile->obj_info[i]); + + /* write out information about each element */ + + for (i = 0; i < plyfile->nelems; i++) { + + elem = plyfile->elems[i]; + fprintf (fp, "element %s %d\n", elem->name, elem->num); + + /* write out each property */ + for (j = 0; j < elem->nprops; j++) { + prop = elem->props[j]; + if (prop->is_list) { + fprintf (fp, "property list "); + write_scalar_type (fp, prop->count_external); + fprintf (fp, " "); + write_scalar_type (fp, prop->external_type); + fprintf (fp, " %s\n", prop->name); + } + else { + fprintf (fp, "property "); + write_scalar_type (fp, prop->external_type); + fprintf (fp, " %s\n", prop->name); + } + } + } + + fprintf (fp, "end_header\n"); +} + + +/****************************************************************************** +Specify which elements are going to be written. This should be called +before a call to the routine ply_put_element(). + + Entry: + plyfile - file identifier + elem_name - name of element we're talking about +******************************************************************************/ + +void ply_put_element_setup(PlyFile *plyfile, char *elem_name) +{ + PlyElement *elem; + + elem = find_element (plyfile, elem_name); + if (elem == NULL) { + fprintf(stderr, "ply_elements_setup: can't find element '%s'\n", elem_name); + exit (-1); + } + + plyfile->which_elem = elem; +} + + +/****************************************************************************** +Write an element to the file. This routine assumes that we're +writing the type of element specified in the last call to the routine +ply_put_element_setup(). + + Entry: + plyfile - file identifier + elem_ptr - pointer to the element +******************************************************************************/ + +void ply_put_element(PlyFile *plyfile, void *elem_ptr) +{ + int j,k; + FILE *fp = plyfile->fp; + PlyElement *elem; + PlyProperty *prop; + char *elem_data,*item; + char **item_ptr; + int list_count; + int item_size; + int int_val; + unsigned int uint_val; + double double_val; + char **other_ptr; + + elem = plyfile->which_elem; + elem_data = (char *)elem_ptr; + other_ptr = (char **) (((char *) elem_ptr) + elem->other_offset); + + /* write out either to an ascii or binary file */ + + if (plyfile->file_type == PLY_ASCII) { + + /* write an ascii file */ + + /* write out each property of the element */ + for (j = 0; j < elem->nprops; j++) { + prop = elem->props[j]; + if (elem->store_prop[j] == OTHER_PROP) + elem_data = *other_ptr; + else + elem_data = (char *)elem_ptr; + if (prop->is_list) { + item = elem_data + prop->count_offset; + get_stored_item ((void *) item, prop->count_internal, + &int_val, &uint_val, &double_val); + write_ascii_item (fp, int_val, uint_val, double_val, + prop->count_external); + list_count = uint_val; + item_ptr = (char **) (elem_data + prop->offset); + item = item_ptr[0]; + item_size = ply_type_size[prop->internal_type]; + for (k = 0; k < list_count; k++) { + get_stored_item ((void *) item, prop->internal_type, + &int_val, &uint_val, &double_val); + write_ascii_item (fp, int_val, uint_val, double_val, + prop->external_type); + item += item_size; + } + } + else { + item = elem_data + prop->offset; + get_stored_item ((void *) item, prop->internal_type, + &int_val, &uint_val, &double_val); + write_ascii_item (fp, int_val, uint_val, double_val, + prop->external_type); + } + } + + fprintf (fp, "\n"); + } + else { + + /* write a binary file */ + + /* write out each property of the element */ + for (j = 0; j < elem->nprops; j++) { + prop = elem->props[j]; + if (elem->store_prop[j] == OTHER_PROP) + elem_data = *other_ptr; + else + elem_data = (char *)elem_ptr; + if (prop->is_list) { + item = elem_data + prop->count_offset; + item_size = ply_type_size[prop->count_internal]; + get_stored_item ((void *) item, prop->count_internal, + &int_val, &uint_val, &double_val); + write_binary_item (fp, plyfile->file_type, int_val, uint_val, + double_val, prop->count_external); + list_count = uint_val; + item_ptr = (char **) (elem_data + prop->offset); + item = item_ptr[0]; + item_size = ply_type_size[prop->internal_type]; + for (k = 0; k < list_count; k++) { + get_stored_item ((void *) item, prop->internal_type, + &int_val, &uint_val, &double_val); + write_binary_item (fp, plyfile->file_type, int_val, uint_val, + double_val, prop->external_type); + item += item_size; + } + } + else { + item = elem_data + prop->offset; + item_size = ply_type_size[prop->internal_type]; + get_stored_item ((void *) item, prop->internal_type, + &int_val, &uint_val, &double_val); + write_binary_item (fp, plyfile->file_type, int_val, uint_val, + double_val, prop->external_type); + } + } + + } +} + + +/****************************************************************************** +Specify a comment that will be written in the header. + + Entry: + plyfile - file identifier + comment - the comment to be written + ******************************************************************************/ + + void ply_put_comment(PlyFile *plyfile, char *comment) + { + /* (re)allocate space for new comment */ + if (plyfile->num_comments == 0) + plyfile->comments = (char **) myalloc (sizeof (char *)); + else + plyfile->comments = (char **) realloc (plyfile->comments, + sizeof (char *) * (plyfile->num_comments + 1)); + + /* add comment to list */ + plyfile->comments[plyfile->num_comments] = _strdup (comment); + plyfile->num_comments++; + } + + + /****************************************************************************** + Specify a piece of object information (arbitrary text) that will be written + in the header. + + Entry: + plyfile - file identifier + obj_info - the text information to be written + ******************************************************************************/ + + void ply_put_obj_info(PlyFile *plyfile, char *obj_info) + { + /* (re)allocate space for new info */ + if (plyfile->num_obj_info == 0) + plyfile->obj_info = (char **) myalloc (sizeof (char *)); + else + plyfile->obj_info = (char **) realloc (plyfile->obj_info, + sizeof (char *) * (plyfile->num_obj_info + 1)); + + /* add info to list */ + plyfile->obj_info[plyfile->num_obj_info] = _strdup (obj_info); + plyfile->num_obj_info++; + } + + + + + + + + /*************/ + /* Reading */ + /*************/ + + + + /****************************************************************************** + Given a file pointer, get ready to read PLY data from the file. + + Entry: + fp - the given file pointer + + Exit: + nelems - number of elements in object + elem_names - list of element names + returns a pointer to a PlyFile, used to refer to this file, or NULL if error + ******************************************************************************/ + + PlyFile *ply_read(FILE *fp, int *nelems, char ***elem_names) + { + int i,j; + PlyFile *plyfile; + int nwords; + char **words; + char **elist; + PlyElement *elem; + char *orig_line; + /* check for NULL file pointer */ + if (fp == NULL) + return (NULL); + + if (native_binary_type == -1) + get_native_binary_type(); + if (!types_checked) + check_types(); + /* create record for this object */ + + plyfile = (PlyFile *) myalloc (sizeof (PlyFile)); + plyfile->nelems = 0; + plyfile->comments = NULL; + plyfile->num_comments = 0; + plyfile->obj_info = NULL; + plyfile->num_obj_info = 0; + plyfile->fp = fp; + plyfile->other_elems = NULL; + + /* read and parse the file's header */ + + words = get_words (plyfile->fp, &nwords, &orig_line); + if (!words || !equal_strings (words[0], "ply")) + { + if (words) + free(words); + return (NULL); + } + while (words) { + /* parse words */ + + if (equal_strings (words[0], "format")) { + if (nwords != 3) { + free(words); + return (NULL); + } + if (equal_strings (words[1], "ascii")) + plyfile->file_type = PLY_ASCII; + else if (equal_strings (words[1], "binary_big_endian")) + plyfile->file_type = PLY_BINARY_BE; + else if (equal_strings (words[1], "binary_little_endian")) + plyfile->file_type = PLY_BINARY_LE; + else { + free(words); + return (NULL); + } + plyfile->version = (float)atof (words[2]); + } + else if (equal_strings (words[0], "element")) + add_element (plyfile, words); + else if (equal_strings (words[0], "property")) + add_property (plyfile, words); + else if (equal_strings (words[0], "comment")) + add_comment (plyfile, orig_line); + else if (equal_strings (words[0], "obj_info")) + add_obj_info (plyfile, orig_line); + else if (equal_strings (words[0], "end_header")) { + free(words); + break; + } + + /* free up words space */ + free (words); + + words = get_words (plyfile->fp, &nwords, &orig_line); + } + + /* create tags for each property of each element, to be used */ + /* later to say whether or not to store each property for the user */ + + for (i = 0; i < plyfile->nelems; i++) { + elem = plyfile->elems[i]; + elem->store_prop = (char *) myalloc (sizeof (char) * elem->nprops); + for (j = 0; j < elem->nprops; j++) + elem->store_prop[j] = DONT_STORE_PROP; + elem->other_offset = NO_OTHER_PROPS; /* no "other" props by default */ + } + + /* set return values about the elements */ + + elist = (char **) myalloc (sizeof (char *) * plyfile->nelems); + for (i = 0; i < plyfile->nelems; i++) + elist[i] = _strdup (plyfile->elems[i]->name); + + *elem_names = elist; + *nelems = plyfile->nelems; + + /* return a pointer to the file's information */ + + return (plyfile); +} + + +/****************************************************************************** +Open a polygon file for reading. + + Entry: + filename - name of file to read from + + Exit: + nelems - number of elements in object + elem_names - list of element names + file_type - file type, either ascii or binary + version - version number of PLY file + returns a file identifier, used to refer to this file, or NULL if error + ******************************************************************************/ + + PlyFile *ply_open_for_reading( + char *filename, + int *nelems, + char ***elem_names, + int *file_type, + float *version + ) + { + FILE *fp; + PlyFile *plyfile; + char *name; + + /* tack on the extension .ply, if necessary */ + + name = (char *) myalloc (int(sizeof (char) * (strlen (filename) + 5))); + strcpy (name, filename); + if (strlen (name) < 4 || + strcmp (name + strlen (name) - 4, ".ply") != 0) + strcat (name, ".ply"); + + /* open the file for reading */ + + fp = fopen (name, "rb"); + if (fp == NULL) + return (NULL); + + /* create the PlyFile data structure */ + + plyfile = ply_read (fp, nelems, elem_names); + + /* determine the file type and version */ + + *file_type = plyfile->file_type; + *version = plyfile->version; + + /* return a pointer to the file's information */ + + free(name); + return (plyfile); + } + + + /****************************************************************************** + Get information about a particular element. + + Entry: + plyfile - file identifier + elem_name - name of element to get information about + + Exit: + nelems - number of elements of this type in the file + nprops - number of properties + returns a list of properties, or NULL if the file doesn't contain that elem + ******************************************************************************/ + + PlyProperty **ply_get_element_description( + PlyFile *plyfile, + char *elem_name, + int *nelems, + int *nprops + ) + { + int i; + PlyElement *elem; + PlyProperty *prop; + PlyProperty **prop_list; + + /* find information about the element */ + elem = find_element (plyfile, elem_name); + if (elem == NULL) + return (NULL); + + *nelems = elem->num; + *nprops = elem->nprops; + + /* make a copy of the element's property list */ + prop_list = (PlyProperty **) myalloc (sizeof (PlyProperty *) * elem->nprops); + for (i = 0; i < elem->nprops; i++) { + prop = (PlyProperty *) myalloc (sizeof (PlyProperty)); + copy_property (prop, elem->props[i]); + prop_list[i] = prop; + } + + /* return this duplicate property list */ + return (prop_list); + } + + + /****************************************************************************** + Specify which properties of an element are to be returned. This should be + called before a call to the routine ply_get_element(). + + Entry: + plyfile - file identifier + elem_name - which element we're talking about + nprops - number of properties + prop_list - list of properties + ******************************************************************************/ + + void ply_get_element_setup( + PlyFile *plyfile, + char *elem_name, + int nprops, + PlyProperty *prop_list + ) + { + int i; + PlyElement *elem; + PlyProperty *prop; + int index; + + /* find information about the element */ + elem = find_element (plyfile, elem_name); + plyfile->which_elem = elem; + + /* deposit the property information into the element's description */ + for (i = 0; i < nprops; i++) { + + /* look for actual property */ + prop = find_property (elem, prop_list[i].name, &index); + if (prop == NULL) { + fprintf (stderr, "Warning: Can't find property '%s' in element '%s'\n", + prop_list[i].name, elem_name); + continue; + } + + /* store its description */ + prop->internal_type = prop_list[i].internal_type; + prop->offset = prop_list[i].offset; + prop->count_internal = prop_list[i].count_internal; + prop->count_offset = prop_list[i].count_offset; + + /* specify that the user wants this property */ + elem->store_prop[index] = STORE_PROP; + } + } + + + /****************************************************************************** + Specify a property of an element that is to be returned. This should be + called (usually multiple times) before a call to the routine ply_get_element(). + This routine should be used in preference to the less flexible old routine + called ply_get_element_setup(). + + Entry: + plyfile - file identifier + elem_name - which element we're talking about + prop - property to add to those that will be returned + ******************************************************************************/ + + int ply_get_property( + PlyFile *plyfile, + char *elem_name, + PlyProperty *prop + ) + { + PlyElement *elem; + PlyProperty *prop_ptr; + int index; + + /* find information about the element */ + elem = find_element (plyfile, elem_name); + plyfile->which_elem = elem; + + /* deposit the property information into the element's description */ + + prop_ptr = find_property (elem, prop->name, &index); + if (prop_ptr == NULL) { +// fprintf (stderr, "Warning: Can't find property '%s' in element '%s'\n", +// prop->name, elem_name); +// return; + return 0; + } + prop_ptr->internal_type = prop->internal_type; + prop_ptr->offset = prop->offset; + prop_ptr->count_internal = prop->count_internal; + prop_ptr->count_offset = prop->count_offset; + + /* specify that the user wants this property */ + elem->store_prop[index] = STORE_PROP; + return 1; + } + + + /****************************************************************************** + Read one element from the file. This routine assumes that we're reading + the type of element specified in the last call to the routine + ply_get_element_setup(). + + Entry: + plyfile - file identifier + elem_ptr - pointer to location where the element information should be put + ******************************************************************************/ + + void ply_get_element(PlyFile *plyfile, void *elem_ptr) + { + if (plyfile->file_type == PLY_ASCII) + ascii_get_element (plyfile, (char *) elem_ptr); + else + binary_get_element (plyfile, (char *) elem_ptr); + } + + + /****************************************************************************** + Extract the comments from the header information of a PLY file. + + Entry: + plyfile - file identifier + + Exit: + num_comments - number of comments returned + returns a pointer to a list of comments + ******************************************************************************/ + + char **ply_get_comments(PlyFile *plyfile, int *num_comments) + { + *num_comments = plyfile->num_comments; + return (plyfile->comments); + } + + + /****************************************************************************** + Extract the object information (arbitrary text) from the header information + of a PLY file. + + Entry: + plyfile - file identifier + + Exit: + num_obj_info - number of lines of text information returned + returns a pointer to a list of object info lines + ******************************************************************************/ + + char **ply_get_obj_info(PlyFile *plyfile, int *num_obj_info) + { + *num_obj_info = plyfile->num_obj_info; + return (plyfile->obj_info); + } + + + /****************************************************************************** + Make ready for "other" properties of an element-- those properties that + the user has not explicitly asked for, but that are to be stashed away + in a special structure to be carried along with the element's other + information. + + Entry: + plyfile - file identifier + elem - element for which we want to save away other properties + ******************************************************************************/ + + void setup_other_props(PlyElement *elem) + { + int i; + PlyProperty *prop; + int size = 0; + int type_size; + + /* Examine each property in decreasing order of size. */ + /* We do this so that all data types will be aligned by */ + /* word, half-word, or whatever within the structure. */ + + for (type_size = 8; type_size > 0; type_size /= 2) { + + /* add up the space taken by each property, and save this information */ + /* away in the property descriptor */ + + for (i = 0; i < elem->nprops; i++) { + + /* don't bother with properties we've been asked to store explicitly */ + if (elem->store_prop[i]) + continue; + + prop = elem->props[i]; + + /* internal types will be same as external */ + prop->internal_type = prop->external_type; + prop->count_internal = prop->count_external; + + /* check list case */ + if (prop->is_list) { + + /* pointer to list */ + if (type_size == sizeof (void *)) { + prop->offset = size; + size += sizeof (void *); /* always use size of a pointer here */ + } + + /* count of number of list elements */ + if (type_size == ply_type_size[prop->count_external]) { + prop->count_offset = size; + size += ply_type_size[prop->count_external]; + } + } + /* not list */ + else if (type_size == ply_type_size[prop->external_type]) { + prop->offset = size; + size += ply_type_size[prop->external_type]; + } + } + + } + + /* save the size for the other_props structure */ + elem->other_size = size; + } + + + /****************************************************************************** + Specify that we want the "other" properties of an element to be tucked + away within the user's structure. The user needn't be concerned for how + these properties are stored. + + Entry: + plyfile - file identifier + elem_name - name of element that we want to store other_props in + offset - offset to where other_props will be stored inside user's structure + + Exit: + returns pointer to structure containing description of other_props + ******************************************************************************/ + + PlyOtherProp *ply_get_other_properties( + PlyFile *plyfile, + char *elem_name, + int offset + ) + { + int i; + PlyElement *elem; + PlyOtherProp *other; + PlyProperty *prop; + int nprops; + + /* find information about the element */ + elem = find_element (plyfile, elem_name); + if (elem == NULL) { + fprintf (stderr, "ply_get_other_properties: Can't find element '%s'\n", + elem_name); + return (NULL); + } + + /* remember that this is the "current" element */ + plyfile->which_elem = elem; + + /* save the offset to where to store the other_props */ + elem->other_offset = offset; + + /* place the appropriate pointers, etc. in the element's property list */ + setup_other_props (elem); + + /* create structure for describing other_props */ + other = (PlyOtherProp *) myalloc (sizeof (PlyOtherProp)); + other->name = _strdup (elem_name); + other->size = elem->other_size; + other->props = (PlyProperty **) myalloc (sizeof(PlyProperty) * elem->nprops); + + /* save descriptions of each "other" property */ + nprops = 0; + for (i = 0; i < elem->nprops; i++) { + if (elem->store_prop[i]) + continue; + prop = (PlyProperty *) myalloc (sizeof (PlyProperty)); + copy_property (prop, elem->props[i]); + other->props[nprops] = prop; + nprops++; + } + other->nprops = nprops; + + /* set other_offset pointer appropriately if there are NO other properties */ + if (other->nprops == 0) { + elem->other_offset = NO_OTHER_PROPS; + } + + /* return structure */ + return (other); + } + + + + + /*************************/ + /* Other Element Stuff */ + /*************************/ + + + + + /****************************************************************************** + Grab all the data for an element that a user does not want to explicitly + read in. + + Entry: + plyfile - pointer to file + elem_name - name of element whose data is to be read in + elem_count - number of instances of this element stored in the file + + Exit: + returns pointer to ALL the "other" element data for this PLY file + ******************************************************************************/ + + PlyOtherElems *ply_get_other_element ( + PlyFile *plyfile, + char *elem_name, + int elem_count + ) + { + int i; + PlyElement *elem; + PlyOtherElems *other_elems; + OtherElem *other; + + /* look for appropriate element */ + elem = find_element (plyfile, elem_name); + if (elem == NULL) { + fprintf (stderr, + "ply_get_other_element: can't find element '%s'\n", elem_name); + exit (-1); + } + + /* create room for the new "other" element, initializing the */ + /* other data structure if necessary */ + + if (plyfile->other_elems == NULL) { + plyfile->other_elems = (PlyOtherElems *) myalloc (sizeof (PlyOtherElems)); + other_elems = plyfile->other_elems; + other_elems->other_list = (OtherElem *) myalloc (sizeof (OtherElem)); + other = &(other_elems->other_list[0]); + other_elems->num_elems = 1; + } + else { + other_elems = plyfile->other_elems; + other_elems->other_list = (OtherElem *) realloc (other_elems->other_list, + sizeof (OtherElem) * other_elems->num_elems + 1); + other = &(other_elems->other_list[other_elems->num_elems]); + other_elems->num_elems++; + } + + /* count of element instances in file */ + other->elem_count = elem_count; + + /* save name of element */ + other->elem_name = _strdup (elem_name); + + /* create a list to hold all the current elements */ + other->other_data = (OtherData **) + malloc (sizeof (OtherData *) * other->elem_count); + + /* set up for getting elements */ + other->other_props = ply_get_other_properties (plyfile, elem_name, + offsetof(OtherData,other_props)); + + /* grab all these elements */ + for (i = 0; i < other->elem_count; i++) { + /* grab and element from the file */ + other->other_data[i] = (OtherData *) malloc (sizeof (OtherData)); + ply_get_element (plyfile, (void *) other->other_data[i]); + } + + /* return pointer to the other elements data */ + return (other_elems); + } + + + /****************************************************************************** + Pass along a pointer to "other" elements that we want to save in a given + PLY file. These other elements were presumably read from another PLY file. + + Entry: + plyfile - file pointer in which to store this other element info + other_elems - info about other elements that we want to store + ******************************************************************************/ + + void ply_describe_other_elements ( + PlyFile *plyfile, + PlyOtherElems *other_elems + ) + { + int i; + OtherElem *other; + PlyElement *elem; + + /* ignore this call if there is no other element */ + if (other_elems == NULL) + return; + + /* save pointer to this information */ + plyfile->other_elems = other_elems; + + /* describe the other properties of this element */ + /* store them in the main element list as elements with + only other properties */ + + REALLOCN(plyfile->elems, PlyElement *, + plyfile->nelems, plyfile->nelems + other_elems->num_elems); + for (i = 0; i < other_elems->num_elems; i++) { + other = &(other_elems->other_list[i]); + elem = (PlyElement *) myalloc (sizeof (PlyElement)); + plyfile->elems[plyfile->nelems++] = elem; + elem->name = _strdup (other->elem_name); + elem->num = other->elem_count; + elem->nprops = 0; + ply_describe_other_properties (plyfile, other->other_props, + offsetof(OtherData,other_props)); + } + } + + + /****************************************************************************** + Write out the "other" elements specified for this PLY file. + + Entry: + plyfile - pointer to PLY file to write out other elements for + ******************************************************************************/ + + void ply_put_other_elements (PlyFile *plyfile) + { + int i,j; + OtherElem *other; + + /* make sure we have other elements to write */ + if (plyfile->other_elems == NULL) + return; + + /* write out the data for each "other" element */ + + for (i = 0; i < plyfile->other_elems->num_elems; i++) { + + other = &(plyfile->other_elems->other_list[i]); + ply_put_element_setup (plyfile, other->elem_name); + + /* write out each instance of the current element */ + for (j = 0; j < other->elem_count; j++) + ply_put_element (plyfile, (void *) other->other_data[j]); + } + } + + + /****************************************************************************** + Free up storage used by an "other" elements data structure. + + Entry: + other_elems - data structure to free up + ******************************************************************************/ + + void ply_free_other_elements (PlyOtherElems *other_elems) + { + other_elems = other_elems; + } + + + + /*******************/ + /* Miscellaneous */ + /*******************/ + + + + /****************************************************************************** + Close a PLY file. + + Entry: + plyfile - identifier of file to close + ******************************************************************************/ + + void ply_close(PlyFile *plyfile) + { + fclose (plyfile->fp); + + /* free up memory associated with the PLY file */ + free (plyfile); + } + + + /****************************************************************************** + Get version number and file type of a PlyFile. + + Entry: + ply - pointer to PLY file + + Exit: + version - version of the file + file_type - PLY_ASCII, PLY_BINARY_BE, or PLY_BINARY_LE + ******************************************************************************/ + + void ply_get_info(PlyFile *ply, float *version, int *file_type) + { + if (ply == NULL) + return; + + *version = ply->version; + *file_type = ply->file_type; + } + + + /****************************************************************************** + Compare two strings. Returns 1 if they are the same, 0 if not. + ******************************************************************************/ + + int equal_strings(const char *s1, const char *s2) + { + + while (*s1 && *s2) + if (*s1++ != *s2++) + return (0); + + if (*s1 != *s2) + return (0); + else + return (1); + } + + + /****************************************************************************** + Find an element from the element list of a given PLY object. + + Entry: + plyfile - file id for PLY file + element - name of element we're looking for + + Exit: + returns the element, or NULL if not found + ******************************************************************************/ + + PlyElement *find_element(PlyFile *plyfile, char *element) + { + int i; + + for (i = 0; i < plyfile->nelems; i++) + if (equal_strings (element, plyfile->elems[i]->name)) + return (plyfile->elems[i]); + + return (NULL); + } + + + /****************************************************************************** + Find a property in the list of properties of a given element. + + Entry: + elem - pointer to element in which we want to find the property + prop_name - name of property to find + + Exit: + index - index to position in list + returns a pointer to the property, or NULL if not found + ******************************************************************************/ + + PlyProperty *find_property(PlyElement *elem, char *prop_name, int *index) + { + int i; + + for (i = 0; i < elem->nprops; i++) + if (equal_strings (prop_name, elem->props[i]->name)) { + *index = i; + return (elem->props[i]); + } + + *index = -1; + return (NULL); + } + + + /****************************************************************************** + Read an element from an ascii file. + + Entry: + plyfile - file identifier + elem_ptr - pointer to element + ******************************************************************************/ + + void ascii_get_element(PlyFile *plyfile, char *elem_ptr) + { + int j,k; + PlyElement *elem; + PlyProperty *prop; + char **words; + int nwords; + int which_word; + char *elem_data,*item=NULL; + char *item_ptr; + int item_size; + int int_val; + unsigned int uint_val; + double double_val; + int list_count; + int store_it; + char **store_array; + char *orig_line; + char *other_data=NULL; + int other_flag; + + /* the kind of element we're reading currently */ + elem = plyfile->which_elem; + + /* do we need to setup for other_props? */ + + if (elem->other_offset != NO_OTHER_PROPS) { + char **ptr; + other_flag = 1; + /* make room for other_props */ + other_data = (char *) myalloc (elem->other_size); + /* store pointer in user's structure to the other_props */ + ptr = (char **) (elem_ptr + elem->other_offset); + *ptr = other_data; + } + else + other_flag = 0; + + /* read in the element */ + + words = get_words (plyfile->fp, &nwords, &orig_line); + if (words == NULL) { + fprintf (stderr, "ply_get_element: unexpected end of file\n"); + exit (-1); + } + + which_word = 0; + + for (j = 0; j < elem->nprops; j++) { + + prop = elem->props[j]; + store_it = (elem->store_prop[j] | other_flag); + + /* store either in the user's structure or in other_props */ + if (elem->store_prop[j]) + elem_data = elem_ptr; + else + elem_data = other_data; + + if (prop->is_list) { /* a list */ + + /* get and store the number of items in the list */ + get_ascii_item (words[which_word++], prop->count_external, + &int_val, &uint_val, &double_val); + if (store_it) { + item = elem_data + prop->count_offset; + store_item(item, prop->count_internal, int_val, uint_val, double_val); + } + + /* allocate space for an array of items and store a ptr to the array */ + list_count = int_val; + item_size = ply_type_size[prop->internal_type]; + store_array = (char **) (elem_data + prop->offset); + + if (list_count == 0) { + if (store_it) + *store_array = NULL; + } + else { + if (store_it) { + item_ptr = (char *) myalloc (sizeof (char) * item_size * list_count); + item = item_ptr; + *store_array = item_ptr; + } + + /* read items and store them into the array */ + for (k = 0; k < list_count; k++) { + get_ascii_item (words[which_word++], prop->external_type, + &int_val, &uint_val, &double_val); + if (store_it) { + store_item (item, prop->internal_type, + int_val, uint_val, double_val); + item += item_size; + } + } + } + + } + else { /* not a list */ + get_ascii_item (words[which_word++], prop->external_type, + &int_val, &uint_val, &double_val); + if (store_it) { + item = elem_data + prop->offset; + store_item (item, prop->internal_type, int_val, uint_val, double_val); + } + } + + } + + free (words); +} + + +/****************************************************************************** +Read an element from a binary file. + + Entry: + plyfile - file identifier + elem_ptr - pointer to an element + ******************************************************************************/ + + void binary_get_element(PlyFile *plyfile, char *elem_ptr) + { + int j,k; + PlyElement *elem; + PlyProperty *prop; + FILE *fp = plyfile->fp; + char *elem_data,*item=NULL; + char *item_ptr; + int item_size; + int int_val; + unsigned int uint_val; + double double_val; + int list_count; + int store_it; + char **store_array; + char *other_data=NULL; + int other_flag; + + /* the kind of element we're reading currently */ + elem = plyfile->which_elem; + + /* do we need to setup for other_props? */ + + if (elem->other_offset != NO_OTHER_PROPS) { + char **ptr; + other_flag = 1; + /* make room for other_props */ + other_data = (char *) myalloc (elem->other_size); + /* store pointer in user's structure to the other_props */ + ptr = (char **) (elem_ptr + elem->other_offset); + *ptr = other_data; + } + else + other_flag = 0; + + /* read in a number of elements */ + + for (j = 0; j < elem->nprops; j++) { + + prop = elem->props[j]; + store_it = (elem->store_prop[j] | other_flag); + + /* store either in the user's structure or in other_props */ + if (elem->store_prop[j]) + elem_data = elem_ptr; + else + elem_data = other_data; + + if (prop->is_list) { /* a list */ + + /* get and store the number of items in the list */ + get_binary_item (fp, plyfile->file_type, prop->count_external, + &int_val, &uint_val, &double_val); + if (store_it) { + item = elem_data + prop->count_offset; + store_item(item, prop->count_internal, int_val, uint_val, double_val); + } + + /* allocate space for an array of items and store a ptr to the array */ + list_count = int_val; + item_size = ply_type_size[prop->internal_type]; + store_array = (char **) (elem_data + prop->offset); + if (list_count == 0) { + if (store_it) + *store_array = NULL; + } + else { + if (store_it) { + item_ptr = (char *) myalloc (sizeof (char) * item_size * list_count); + item = item_ptr; + *store_array = item_ptr; + } + + /* read items and store them into the array */ + for (k = 0; k < list_count; k++) { + get_binary_item (fp, plyfile->file_type, prop->external_type, + &int_val, &uint_val, &double_val); + if (store_it) { + store_item (item, prop->internal_type, + int_val, uint_val, double_val); + item += item_size; + } + } + } + + } + else { /* not a list */ + get_binary_item (fp, plyfile->file_type, prop->external_type, + &int_val, &uint_val, &double_val); + if (store_it) { + item = elem_data + prop->offset; + store_item (item, prop->internal_type, int_val, uint_val, double_val); + } + } + + } + } + + + /****************************************************************************** + Write to a file the word that represents a PLY data type. + + Entry: + fp - file pointer + code - code for type + ******************************************************************************/ + + void write_scalar_type (FILE *fp, int code) + { + /* make sure this is a valid code */ + + if (code <= PLY_START_TYPE || code >= PLY_END_TYPE) { + fprintf (stderr, "write_scalar_type: bad data code = %d\n", code); + exit (-1); + } + + /* write the code to a file */ + + fprintf (fp, "%s", type_names[code]); + } + + /****************************************************************************** + Reverse the order in an array of bytes. This is the conversion from big + endian to little endian and vice versa + + Entry: + bytes - array of bytes to reverse (in place) + num_bytes - number of bytes in array + ******************************************************************************/ + + void swap_bytes(char *bytes, int num_bytes) + { + int i; + char temp; + + for (i=0; i < num_bytes/2; i++) + { + temp = bytes[i]; + bytes[i] = bytes[(num_bytes-1)-i]; + bytes[(num_bytes-1)-i] = temp; + } + } + + /****************************************************************************** + Find out if this machine is big endian or little endian + + Exit: + set global variable, native_binary_type = + either PLY_BINARY_BE or PLY_BINARY_LE + + ******************************************************************************/ + + void get_native_binary_type() + { + endian_test_type test; + + test.int_value = 0; + test.int_value = 1; + if (test.byte_values[0] == 1) + native_binary_type = PLY_BINARY_LE; + else if (test.byte_values[sizeof(int)-1] == 1) + native_binary_type = PLY_BINARY_BE; + else + { + fprintf(stderr, "ply: Couldn't determine machine endianness.\n"); + fprintf(stderr, "ply: Exiting...\n"); + exit(1); + } + } + + /****************************************************************************** + Verify that all the native types are the sizes we need + + + ******************************************************************************/ + + void check_types() + { + if ((ply_type_size[PLY_CHAR] != sizeof(char)) || + (ply_type_size[PLY_SHORT] != sizeof(short)) || + (ply_type_size[PLY_INT] != sizeof(int)) || + (ply_type_size[PLY_UCHAR] != sizeof(unsigned char)) || + (ply_type_size[PLY_USHORT] != sizeof(unsigned short)) || + (ply_type_size[PLY_UINT] != sizeof(unsigned int)) || + (ply_type_size[PLY_FLOAT] != sizeof(float)) || + (ply_type_size[PLY_DOUBLE] != sizeof(double))) + { + fprintf(stderr, "ply: Type sizes do not match built-in types\n"); + fprintf(stderr, "ply: Exiting...\n"); + exit(1); + } + + types_checked = 1; + } + + /****************************************************************************** + Get a text line from a file and break it up into words. + + IMPORTANT: The calling routine call "free" on the returned pointer once + finished with it. + + Entry: + fp - file to read from + + Exit: + nwords - number of words returned + orig_line - the original line of characters + returns a list of words from the line, or NULL if end-of-file + ******************************************************************************/ + + char **get_words(FILE *fp, int *nwords, char **orig_line) + { +#define BIG_STRING 4096 + static char str[BIG_STRING]; + static char str_copy[BIG_STRING]; + char **words; + int max_words = 10; + int num_words = 0; + char *ptr,*ptr2; + char *result; + + words = (char **) myalloc (sizeof (char *) * max_words); + + /* read in a line */ + result = fgets (str, BIG_STRING, fp); + if (result == NULL) { + *nwords = 0; + *orig_line = NULL; + return (NULL); + } + /* convert line-feed and tabs into spaces */ + /* (this guarentees that there will be a space before the */ + /* null character at the end of the string) */ + + str[BIG_STRING-2] = ' '; + str[BIG_STRING-1] = '\0'; + + for (ptr = str, ptr2 = str_copy; *ptr != '\0'; ptr++, ptr2++) { + *ptr2 = *ptr; + // Added line here to manage carriage returns + if (*ptr == '\t' || *ptr == '\r') { + *ptr = ' '; + *ptr2 = ' '; + } + else if (*ptr == '\n') { + *ptr = ' '; + *ptr2 = '\0'; + break; + } + } + + /* find the words in the line */ + + ptr = str; + while (*ptr != '\0') { + + /* jump over leading spaces */ + while (*ptr == ' ') + ptr++; + + /* break if we reach the end */ + if (*ptr == '\0') + break; + + /* save pointer to beginning of word */ + if (num_words >= max_words) { + max_words += 10; + words = (char **) realloc (words, sizeof (char *) * max_words); + } + words[num_words++] = ptr; + + /* jump over non-spaces */ + while (*ptr != ' ') + ptr++; + + /* place a null character here to mark the end of the word */ + *ptr++ = '\0'; + } + + /* return the list of words */ + *nwords = num_words; + *orig_line = str_copy; + return (words); + } + + + /****************************************************************************** + Return the value of an item, given a pointer to it and its type. + + Entry: + item - pointer to item + type - data type that "item" points to + + Exit: + returns a double-precision float that contains the value of the item + ******************************************************************************/ + + double get_item_value(char *item, int type) + { + unsigned char *puchar; + char *pchar; + short int *pshort; + unsigned short int *pushort; + int *pint; + unsigned int *puint; + float *pfloat; + double *pdouble; + int int_value; + unsigned int uint_value; + double double_value; + + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + pchar = (char *) item; + int_value = *pchar; + return ((double) int_value); + case PLY_UCHAR: + case PLY_UINT_8: + puchar = (unsigned char *) item; + int_value = *puchar; + return ((double) int_value); + case PLY_SHORT: + case PLY_INT_16: + pshort = (short int *) item; + int_value = *pshort; + return ((double) int_value); + case PLY_USHORT: + case PLY_UINT_16: + pushort = (unsigned short int *) item; + int_value = *pushort; + return ((double) int_value); + case PLY_INT: + case PLY_INT_32: + pint = (int *) item; + int_value = *pint; + return ((double) int_value); + case PLY_UINT: + case PLY_UINT_32: + puint = (unsigned int *) item; + uint_value = *puint; + return ((double) uint_value); + case PLY_FLOAT: + case PLY_FLOAT_32: + pfloat = (float *) item; + double_value = *pfloat; + return (double_value); + case PLY_DOUBLE: + case PLY_FLOAT_64: + pdouble = (double *) item; + double_value = *pdouble; + return (double_value); + default: + fprintf (stderr, "get_item_value: bad type = %d\n", type); + exit (-1); + } + } + + + /****************************************************************************** + Write out an item to a file as raw binary bytes. + + Entry: + fp - file to write to + int_val - integer version of item + uint_val - unsigned integer version of item + double_val - double-precision float version of item + type - data type to write out + ******************************************************************************/ + + void write_binary_item( + FILE *fp, + int file_type, + int int_val, + unsigned int uint_val, + double double_val, + int type + ) + { + unsigned char uchar_val; + char char_val; + unsigned short ushort_val; + short short_val; + float float_val; + void *value; + + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + char_val = char(int_val); + value = &char_val; + break; + case PLY_SHORT: + case PLY_INT_16: + short_val = short(int_val); + value = &short_val; + break; + case PLY_INT: + case PLY_INT_32: + value = &int_val; + break; + case PLY_UCHAR: + case PLY_UINT_8: + uchar_val = (unsigned char)(uint_val); + value = &uchar_val; + break; + case PLY_USHORT: + case PLY_UINT_16: + ushort_val = (unsigned short)(uint_val); + value = &ushort_val; + break; + case PLY_UINT: + case PLY_UINT_32: + value = &uint_val; + break; + case PLY_FLOAT: + case PLY_FLOAT_32: + float_val = (float)double_val; + value = &float_val; + break; + case PLY_DOUBLE: + case PLY_FLOAT_64: + value = &double_val; + break; + default: + fprintf (stderr, "write_binary_item: bad type = %d\n", type); + exit (-1); + } + + + if ((file_type != native_binary_type) && (ply_type_size[type] > 1)) + swap_bytes((char *)value, ply_type_size[type]); + + if (fwrite (value, ply_type_size[type], 1, fp) != 1) + { + fprintf(stderr, "PLY ERROR: fwrite() failed -- aborting.\n"); + exit(1); + } + } + + + /****************************************************************************** + Write out an item to a file as ascii characters. + + Entry: + fp - file to write to + int_val - integer version of item + uint_val - unsigned integer version of item + double_val - double-precision float version of item + type - data type to write out + ******************************************************************************/ + + void write_ascii_item( + FILE *fp, + int int_val, + unsigned int uint_val, + double double_val, + int type + ) + { + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + case PLY_SHORT: + case PLY_INT_16: + case PLY_INT: + case PLY_INT_32: + if (fprintf (fp, "%d ", int_val) <= 0) + { + fprintf(stderr, "PLY ERROR: fprintf() failed -- aborting.\n"); + exit(1); + } + break; + case PLY_UCHAR: + case PLY_UINT_8: + case PLY_USHORT: + case PLY_UINT_16: + case PLY_UINT: + case PLY_UINT_32: + + if (fprintf (fp, "%u ", uint_val) <= 0) + { + fprintf(stderr, "PLY ERROR: fprintf() failed -- aborting.\n"); + exit(1); + } + break; + case PLY_FLOAT: + case PLY_FLOAT_32: + case PLY_DOUBLE: + case PLY_FLOAT_64: + if (fprintf (fp, "%g ", double_val) <= 0) + { + fprintf(stderr, "PLY ERROR: fprintf() failed -- aborting.\n"); + exit(1); + } + break; + default: + fprintf (stderr, "write_ascii_item: bad type = %d\n", type); + exit (-1); + } + } + + + /****************************************************************************** + Write out an item to a file as ascii characters. + + Entry: + fp - file to write to + item - pointer to item to write + type - data type that "item" points to + + Exit: + returns a double-precision float that contains the value of the written item + ******************************************************************************/ + + double old_write_ascii_item(FILE *fp, char *item, int type) + { + unsigned char *puchar; + char *pchar; + short int *pshort; + unsigned short int *pushort; + int *pint; + unsigned int *puint; + float *pfloat; + double *pdouble; + int int_value; + unsigned int uint_value; + double double_value; + + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + pchar = (char *) item; + int_value = *pchar; + fprintf (fp, "%d ", int_value); + return ((double) int_value); + case PLY_UCHAR: + case PLY_UINT_8: + puchar = (unsigned char *) item; + int_value = *puchar; + fprintf (fp, "%d ", int_value); + return ((double) int_value); + case PLY_SHORT: + case PLY_INT_16: + pshort = (short int *) item; + int_value = *pshort; + fprintf (fp, "%d ", int_value); + return ((double) int_value); + case PLY_USHORT: + case PLY_UINT_16: + pushort = (unsigned short int *) item; + int_value = *pushort; + fprintf (fp, "%d ", int_value); + return ((double) int_value); + case PLY_INT: + case PLY_INT_32: + pint = (int *) item; + int_value = *pint; + fprintf (fp, "%d ", int_value); + return ((double) int_value); + case PLY_UINT: + case PLY_UINT_32: + puint = (unsigned int *) item; + uint_value = *puint; + fprintf (fp, "%u ", uint_value); + return ((double) uint_value); + case PLY_FLOAT: + case PLY_FLOAT_32: + pfloat = (float *) item; + double_value = *pfloat; + fprintf (fp, "%g ", double_value); + return (double_value); + case PLY_DOUBLE: + case PLY_FLOAT_64: + pdouble = (double *) item; + double_value = *pdouble; + fprintf (fp, "%g ", double_value); + return (double_value); + default: + fprintf (stderr, "old_write_ascii_item: bad type = %d\n", type); + exit (-1); + } + } + + + /****************************************************************************** + Get the value of an item that is in memory, and place the result + into an integer, an unsigned integer and a double. + + Entry: + ptr - pointer to the item + type - data type supposedly in the item + + Exit: + int_val - integer value + uint_val - unsigned integer value + double_val - double-precision floating point value + ******************************************************************************/ + + void get_stored_item( + void *ptr, + int type, + int *int_val, + unsigned int *uint_val, + double *double_val + ) + { + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + *int_val = *((char *) ptr); + *uint_val = *int_val; + *double_val = *int_val; + break; + case PLY_UCHAR: + case PLY_UINT_8: + *uint_val = *((unsigned char *) ptr); + *int_val = *uint_val; + *double_val = *uint_val; + break; + case PLY_SHORT: + case PLY_INT_16: + *int_val = *((short int *) ptr); + *uint_val = *int_val; + *double_val = *int_val; + break; + case PLY_USHORT: + case PLY_UINT_16: + *uint_val = *((unsigned short int *) ptr); + *int_val = *uint_val; + *double_val = *uint_val; + break; + case PLY_INT: + case PLY_INT_32: + *int_val = *((int *) ptr); + *uint_val = *int_val; + *double_val = *int_val; + break; + case PLY_UINT: + case PLY_UINT_32: + *uint_val = *((unsigned int *) ptr); + *int_val = *uint_val; + *double_val = *uint_val; + break; + case PLY_FLOAT: + case PLY_FLOAT_32: + *double_val = *((float *) ptr); + *int_val = (int) *double_val; + *uint_val = (unsigned int) *double_val; + break; + case PLY_DOUBLE: + case PLY_FLOAT_64: + *double_val = *((double *) ptr); + *int_val = (int) *double_val; + *uint_val = (unsigned int) *double_val; + break; + default: + fprintf (stderr, "get_stored_item: bad type = %d\n", type); + exit (-1); + } + } + + + /****************************************************************************** + Get the value of an item from a binary file, and place the result + into an integer, an unsigned integer and a double. + + Entry: + fp - file to get item from + type - data type supposedly in the word + + Exit: + int_val - integer value + uint_val - unsigned integer value + double_val - double-precision floating point value + ******************************************************************************/ + + void get_binary_item( + FILE *fp, + int file_type, + int type, + int *int_val, + unsigned int *uint_val, + double *double_val + ) + { + char c[8]; + void *ptr; + + ptr = (void *) c; + + if (fread (ptr, ply_type_size[type], 1, fp) != 1) + { + fprintf(stderr, "PLY ERROR: fread() failed -- aborting.\n"); + exit(1); + } + + + if ((file_type != native_binary_type) && (ply_type_size[type] > 1)) + swap_bytes((char *)ptr, ply_type_size[type]); + + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + *int_val = *((char *) ptr); + *uint_val = *int_val; + *double_val = *int_val; + break; + case PLY_UCHAR: + case PLY_UINT_8: + *uint_val = *((unsigned char *) ptr); + *int_val = *uint_val; + *double_val = *uint_val; + break; + case PLY_SHORT: + case PLY_INT_16: + *int_val = *((short int *) ptr); + *uint_val = *int_val; + *double_val = *int_val; + break; + case PLY_USHORT: + case PLY_UINT_16: + *uint_val = *((unsigned short int *) ptr); + *int_val = *uint_val; + *double_val = *uint_val; + break; + case PLY_INT: + case PLY_INT_32: + *int_val = *((int *) ptr); + *uint_val = *int_val; + *double_val = *int_val; + break; + case PLY_UINT: + case PLY_UINT_32: + *uint_val = *((unsigned int *) ptr); + *int_val = *uint_val; + *double_val = *uint_val; + break; + case PLY_FLOAT: + case PLY_FLOAT_32: + *double_val = *((float *) ptr); + *int_val = (int) *double_val; + *uint_val = (unsigned int) *double_val; + break; + case PLY_DOUBLE: + case PLY_FLOAT_64: + *double_val = *((double *) ptr); + *int_val = (int) *double_val; + *uint_val = (unsigned int) *double_val; + break; + default: + fprintf (stderr, "get_binary_item: bad type = %d\n", type); + exit (-1); + } + } + + + /****************************************************************************** + Extract the value of an item from an ascii word, and place the result + into an integer, an unsigned integer and a double. + + Entry: + word - word to extract value from + type - data type supposedly in the word + + Exit: + int_val - integer value + uint_val - unsigned integer value + double_val - double-precision floating point value + ******************************************************************************/ + + void get_ascii_item( + char *word, + int type, + int *int_val, + unsigned int *uint_val, + double *double_val + ) + { + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + case PLY_UCHAR: + case PLY_UINT_8: + case PLY_SHORT: + case PLY_INT_16: + case PLY_USHORT: + case PLY_UINT_16: + case PLY_INT: + case PLY_INT_32: + *int_val = atoi (word); + *uint_val = (unsigned int) *int_val; + *double_val = (double) *int_val; + break; + + case PLY_UINT: + case PLY_UINT_32: + *uint_val = strtol (word, (char **) NULL, 10); + *int_val = (int) *uint_val; + *double_val = (double) *uint_val; + break; + + case PLY_FLOAT: + case PLY_FLOAT_32: + case PLY_DOUBLE: + case PLY_FLOAT_64: + *double_val = atof (word); + *int_val = (int) *double_val; + *uint_val = (unsigned int) *double_val; + break; + + default: + fprintf (stderr, "get_ascii_item: bad type = %d\n", type); + exit (-1); + } + } + + + /****************************************************************************** + Store a value into a place being pointed to, guided by a data type. + + Entry: + item - place to store value + type - data type + int_val - integer version of value + uint_val - unsigned integer version of value + double_val - double version of value + + Exit: + item - pointer to stored value + ******************************************************************************/ + + void store_item ( + char *item, + int type, + int int_val, + unsigned int uint_val, + double double_val + ) + { + unsigned char *puchar; + short int *pshort; + unsigned short int *pushort; + int *pint; + unsigned int *puint; + float *pfloat; + double *pdouble; + + + switch (type) { + case PLY_CHAR: + case PLY_INT_8: + *item = char(int_val); + break; + case PLY_UCHAR: + case PLY_UINT_8: + puchar = (unsigned char *) item; + *puchar = (unsigned char)(uint_val); + break; + case PLY_SHORT: + case PLY_INT_16: + pshort = (short *) item; + *pshort = short(int_val); + break; + case PLY_USHORT: + case PLY_UINT_16: + pushort = (unsigned short *) item; + *pushort = (unsigned short)(uint_val); + break; + case PLY_INT: + case PLY_INT_32: + pint = (int *) item; + *pint = int_val; + break; + case PLY_UINT: + case PLY_UINT_32: + puint = (unsigned int *) item; + *puint = uint_val; + break; + case PLY_FLOAT: + case PLY_FLOAT_32: + pfloat = (float *) item; + *pfloat = (float)double_val; + break; + case PLY_DOUBLE: + case PLY_FLOAT_64: + pdouble = (double *) item; + *pdouble = double_val; + break; + default: + fprintf (stderr, "store_item: bad type = %d\n", type); + exit (-1); + } + } + + + /****************************************************************************** + Add an element to a PLY file descriptor. + + Entry: + plyfile - PLY file descriptor + words - list of words describing the element + nwords - number of words in the list + ******************************************************************************/ + + void add_element (PlyFile *plyfile, char **words) + { + PlyElement *elem; + + /* create the new element */ + elem = (PlyElement *) myalloc (sizeof (PlyElement)); + elem->name = _strdup (words[1]); + elem->num = atoi (words[2]); + elem->nprops = 0; + + /* make room for new element in the object's list of elements */ + if (plyfile->nelems == 0) + plyfile->elems = (PlyElement **) myalloc (sizeof (PlyElement *)); + else + plyfile->elems = (PlyElement **) realloc (plyfile->elems, + sizeof (PlyElement *) * (plyfile->nelems + 1)); + + /* add the new element to the object's list */ + plyfile->elems[plyfile->nelems] = elem; + plyfile->nelems++; + } + + + /****************************************************************************** + Return the type of a property, given the name of the property. + + Entry: + name - name of property type + + Exit: + returns integer code for property, or 0 if not found + ******************************************************************************/ + + int get_prop_type(char *type_name) + { + int i; + + for (i = PLY_START_TYPE + 1; i < PLY_END_TYPE; i++) + if (equal_strings (type_name, type_names[i])) + return (i); + + /* if we get here, we didn't find the type */ + return (0); + } + + + /****************************************************************************** + Add a property to a PLY file descriptor. + + Entry: + plyfile - PLY file descriptor + words - list of words describing the property + nwords - number of words in the list + ******************************************************************************/ + + void add_property (PlyFile *plyfile, char **words) + { + PlyProperty *prop; + PlyElement *elem; + + /* create the new property */ + + prop = (PlyProperty *) myalloc (sizeof (PlyProperty)); + + if (equal_strings (words[1], "list")) { /* is a list */ + prop->count_external = get_prop_type (words[2]); + prop->external_type = get_prop_type (words[3]); + prop->name = _strdup (words[4]); + prop->is_list = 1; + } + else { /* not a list */ + prop->external_type = get_prop_type (words[1]); + prop->name = _strdup (words[2]); + prop->is_list = 0; + } + + /* add this property to the list of properties of the current element */ + + elem = plyfile->elems[plyfile->nelems - 1]; + + if (elem->nprops == 0) + elem->props = (PlyProperty **) myalloc (sizeof (PlyProperty *)); + else + elem->props = (PlyProperty **) realloc (elem->props, + sizeof (PlyProperty *) * (elem->nprops + 1)); + + elem->props[elem->nprops] = prop; + elem->nprops++; + } + + + /****************************************************************************** + Add a comment to a PLY file descriptor. + + Entry: + plyfile - PLY file descriptor + line - line containing comment + ******************************************************************************/ + + void add_comment (PlyFile *plyfile, char *line) + { + int i; + + /* skip over "comment" and leading spaces and tabs */ + i = 7; + while (line[i] == ' ' || line[i] == '\t') + i++; + + ply_put_comment (plyfile, &line[i]); + } + + + /****************************************************************************** + Add a some object information to a PLY file descriptor. + + Entry: + plyfile - PLY file descriptor + line - line containing text info + ******************************************************************************/ + + void add_obj_info (PlyFile *plyfile, char *line) + { + int i; + + /* skip over "obj_info" and leading spaces and tabs */ + i = 8; + while (line[i] == ' ' || line[i] == '\t') + i++; + + ply_put_obj_info (plyfile, &line[i]); + } + + + /****************************************************************************** + Copy a property. + ******************************************************************************/ + + void copy_property(PlyProperty *dest, PlyProperty *src) + { + dest->name = _strdup (src->name); + dest->external_type = src->external_type; + dest->internal_type = src->internal_type; + dest->offset = src->offset; + + dest->is_list = src->is_list; + dest->count_external = src->count_external; + dest->count_internal = src->count_internal; + dest->count_offset = src->count_offset; + } + + + /****************************************************************************** + Allocate some memory. + + Entry: + size - amount of memory requested (in bytes) + lnum - line number from which memory was requested + fname - file name from which memory was requested + ******************************************************************************/ + + char *my_alloc(int size, int lnum, char *fname) + { + char *ptr; + + ptr = (char *) malloc (size); + + if (ptr == 0) { + fprintf(stderr, "Memory allocation bombed on line %d in %s\n", lnum, fname); + } + + return (ptr); + } + diff --git a/src/plugins_experimental/filter_screened_poisson/Src/PlyFile.h b/src/plugins_experimental/filter_screened_poisson/Src/PlyFile.h new file mode 100755 index 000000000..1c72360b8 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/PlyFile.h @@ -0,0 +1,221 @@ +/* + + Header for PLY polygon files. + + - Greg Turk, March 1994 + + A PLY file contains a single polygonal _object_. + + An object is composed of lists of _elements_. Typical elements are + vertices, faces, edges and materials. + + Each type of element for a given object has one or more _properties_ + associated with the element type. For instance, a vertex element may + have as properties three floating-point values x,y,z and three unsigned + chars for red, green and blue. + + --------------------------------------------------------------- + + Copyright (c) 1994 The Board of Trustees of The Leland Stanford + Junior University. All rights reserved. + + Permission to use, copy, modify and distribute this software and its + documentation for any purpose is hereby granted without fee, provided + that the above copyright notice and this permission notice appear in + all copies of this software and that you do not sell the software. + + THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND, + EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY + WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. + +*/ + +#ifndef PLY_FILE_INCLUDED +#define PLY_FILE_INCLUDED + + +#ifndef WIN32 +#define _strdup strdup +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +#include +#include +#include +#include + +#define PLY_ASCII 1 /* ascii PLY file */ +#define PLY_BINARY_BE 2 /* binary PLY file, big endian */ +#define PLY_BINARY_LE 3 /* binary PLY file, little endian */ +#define PLY_BINARY_NATIVE 4 /* binary PLY file, same endianness as current architecture */ + +#define PLY_OKAY 0 /* ply routine worked okay */ +#define PLY_ERROR -1 /* error in ply routine */ + + /* scalar data types supported by PLY format */ + +#define PLY_START_TYPE 0 +#define PLY_CHAR 1 +#define PLY_SHORT 2 +#define PLY_INT 3 +#define PLY_UCHAR 4 +#define PLY_USHORT 5 +#define PLY_UINT 6 +#define PLY_FLOAT 7 +#define PLY_DOUBLE 8 +#define PLY_INT_8 9 +#define PLY_UINT_8 10 +#define PLY_INT_16 11 +#define PLY_UINT_16 12 +#define PLY_INT_32 13 +#define PLY_UINT_32 14 +#define PLY_FLOAT_32 15 +#define PLY_FLOAT_64 16 + +#define PLY_END_TYPE 17 + +#define PLY_SCALAR 0 +#define PLY_LIST 1 + +#define PLY_STRIP_COMMENT_HEADER 0 + +typedef struct PlyProperty { /* description of a property */ + + char *name; /* property name */ + int external_type; /* file's data type */ + int internal_type; /* program's data type */ + int offset; /* offset bytes of prop in a struct */ + + int is_list; /* 1 = list, 0 = scalar */ + int count_external; /* file's count type */ + int count_internal; /* program's count type */ + int count_offset; /* offset byte for list count */ + +} PlyProperty; + +typedef struct PlyElement { /* description of an element */ + char *name; /* element name */ + int num; /* number of elements in this object */ + int size; /* size of element (bytes) or -1 if variable */ + int nprops; /* number of properties for this element */ + PlyProperty **props; /* list of properties in the file */ + char *store_prop; /* flags: property wanted by user? */ + int other_offset; /* offset to un-asked-for props, or -1 if none*/ + int other_size; /* size of other_props structure */ +} PlyElement; + +typedef struct PlyOtherProp { /* describes other properties in an element */ + char *name; /* element name */ + int size; /* size of other_props */ + int nprops; /* number of properties in other_props */ + PlyProperty **props; /* list of properties in other_props */ +} PlyOtherProp; + +typedef struct OtherData { /* for storing other_props for an other element */ + void *other_props; +} OtherData; + +typedef struct OtherElem { /* data for one "other" element */ + char *elem_name; /* names of other elements */ + int elem_count; /* count of instances of each element */ + OtherData **other_data; /* actual property data for the elements */ + PlyOtherProp *other_props; /* description of the property data */ +} OtherElem; + +typedef struct PlyOtherElems { /* "other" elements, not interpreted by user */ + int num_elems; /* number of other elements */ + OtherElem *other_list; /* list of data for other elements */ +} PlyOtherElems; + +typedef struct PlyFile { /* description of PLY file */ + FILE *fp; /* file pointer */ + int file_type; /* ascii or binary */ + float version; /* version number of file */ + int nelems; /* number of elements of object */ + PlyElement **elems; /* list of elements */ + int num_comments; /* number of comments */ + char **comments; /* list of comments */ + int num_obj_info; /* number of items of object information */ + char **obj_info; /* list of object info items */ + PlyElement *which_elem; /* which element we're currently writing */ + PlyOtherElems *other_elems; /* "other" elements from a PLY file */ +} PlyFile; + + /* memory allocation */ +extern char *my_alloc(); +#define myalloc(mem_size) my_alloc((mem_size), __LINE__, __FILE__) + +#ifndef ALLOCN +#define REALLOCN(PTR,TYPE,OLD_N,NEW_N) \ +{ \ + if ((OLD_N) == 0) \ +{ ALLOCN((PTR),TYPE,(NEW_N));} \ + else \ +{ \ + (PTR) = (TYPE *)realloc((PTR),(NEW_N)*sizeof(TYPE)); \ + if (((PTR) == NULL) && ((NEW_N) != 0)) \ +{ \ + fprintf(stderr, "Memory reallocation failed on line %d in %s\n", \ + __LINE__, __FILE__); \ + fprintf(stderr, " tried to reallocate %d->%d\n", \ + (OLD_N), (NEW_N)); \ + exit(-1); \ +} \ + if ((NEW_N)>(OLD_N)) \ + memset((char *)(PTR)+(OLD_N)*sizeof(TYPE), 0, \ + ((NEW_N)-(OLD_N))*sizeof(TYPE)); \ +} \ +} + +#define ALLOCN(PTR,TYPE,N) \ +{ (PTR) = (TYPE *) calloc(((unsigned)(N)),sizeof(TYPE));\ + if ((PTR) == NULL) { \ + fprintf(stderr, "Memory allocation failed on line %d in %s\n", \ + __LINE__, __FILE__); \ + exit(-1); \ + } \ +} + + +#define FREE(PTR) { free((PTR)); (PTR) = NULL; } +#endif + + +/*** delcaration of routines ***/ + +extern PlyFile *ply_write(FILE *, int, const char **, int); +extern PlyFile *ply_open_for_writing(char *, int, const char **, int, float *); +extern void ply_describe_element(PlyFile *, char *, int, int, PlyProperty *); +extern void ply_describe_property(PlyFile *, char *, PlyProperty *); +extern void ply_element_count(PlyFile *, char *, int); +extern void ply_header_complete(PlyFile *); +extern void ply_put_element_setup(PlyFile *, char *); +extern void ply_put_element(PlyFile *, void *); +extern void ply_put_comment(PlyFile *, char *); +extern void ply_put_obj_info(PlyFile *, char *); +extern PlyFile *ply_read(FILE *, int *, char ***); +extern PlyFile *ply_open_for_reading( char *, int *, char ***, int *, float *); +extern PlyProperty **ply_get_element_description(PlyFile *, char *, int*, int*); +extern void ply_get_element_setup( PlyFile *, char *, int, PlyProperty *); +extern int ply_get_property(PlyFile *, char *, PlyProperty *); +extern PlyOtherProp *ply_get_other_properties(PlyFile *, char *, int); +extern void ply_get_element(PlyFile *, void *); +extern char **ply_get_comments(PlyFile *, int *); +extern char **ply_get_obj_info(PlyFile *, int *); +extern void ply_close(PlyFile *); +extern void ply_get_info(PlyFile *, float *, int *); +extern PlyOtherElems *ply_get_other_element (PlyFile *, char *, int); +extern void ply_describe_other_elements ( PlyFile *, PlyOtherElems *); +extern void ply_put_other_elements (PlyFile *); +extern void ply_free_other_elements (PlyOtherElems *); +extern void ply_describe_other_properties(PlyFile *, PlyOtherProp *, int); + +extern int equal_strings(const char *, const char *); + +#ifdef __cplusplus +} +#endif +#endif // PLY_FILE_INCLUDED diff --git a/src/plugins_experimental/filter_screened_poisson/Src/PointStream.inl b/src/plugins_experimental/filter_screened_poisson/Src/PointStream.inl new file mode 100755 index 000000000..c359c5993 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/PointStream.inl @@ -0,0 +1,183 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ +#include "Ply.h" + +template< class Real > +ASCIIPointStream< Real >::ASCIIPointStream( const char* fileName ) +{ + _fp = fopen( fileName , "r" ); + if( !_fp ) fprintf( stderr , "Failed to open file for reading: %s\n" , fileName ) , exit( 0 ); +} +template< class Real > +ASCIIPointStream< Real >::~ASCIIPointStream( void ) +{ + fclose( _fp ); + _fp = NULL; +} +template< class Real > +void ASCIIPointStream< Real >::reset( void ) { fseek( _fp , SEEK_SET , 0 ); } +template< class Real > +bool ASCIIPointStream< Real >::nextPoint( Point3D< Real >& p , Point3D< Real >& n ) +{ + float c[2*DIMENSION]; + if( fscanf( _fp , " %f %f %f %f %f %f " , &c[0] , &c[1] , &c[2] , &c[3] , &c[4] , &c[5] )!=2*DIMENSION ) return false; + p[0] = c[0] , p[1] = c[1] , p[2] = c[2]; + n[0] = c[3] , n[1] = c[4] , n[2] = c[5]; + return true; +} +template< class Real > +BinaryPointStream< Real >::BinaryPointStream( const char* fileName ) +{ + _pointsInBuffer = _currentPointIndex = 0; + _fp = fopen( fileName , "rb" ); + if( !_fp ) fprintf( stderr , "Failed to open file for reading: %s\n" , fileName ) , exit( 0 ); +} +template< class Real > +BinaryPointStream< Real >::~BinaryPointStream( void ) +{ + fclose( _fp ); + _fp = NULL; +} +template< class Real > +void BinaryPointStream< Real >::reset( void ) +{ + fseek( _fp , SEEK_SET , 0 ); + _pointsInBuffer = _currentPointIndex = 0; +} +template< class Real > +bool BinaryPointStream< Real >::nextPoint( Point3D< Real >& p , Point3D< Real >& n ) +{ + if( _currentPointIndex<_pointsInBuffer ) + { + p[0] = _pointBuffer[ _currentPointIndex*6+0 ]; + p[1] = _pointBuffer[ _currentPointIndex*6+1 ]; + p[2] = _pointBuffer[ _currentPointIndex*6+2 ]; + n[0] = _pointBuffer[ _currentPointIndex*6+3 ]; + n[1] = _pointBuffer[ _currentPointIndex*6+4 ]; + n[2] = _pointBuffer[ _currentPointIndex*6+5 ]; + _currentPointIndex++; + return true; + } + else + { + _currentPointIndex = 0; + _pointsInBuffer = int( fread( _pointBuffer , sizeof( Real ) * 6 , POINT_BUFFER_SIZE , _fp ) ); + if( !_pointsInBuffer ) return false; + else return nextPoint( p , n ); + } +} + +template< class Real > +PLYPointStream< Real >::PLYPointStream( const char* fileName ) +{ + _fileName = new char[ strlen( fileName )+1 ]; + strcpy( _fileName , fileName ); + _ply = NULL; + reset(); +} +template< class Real > +void PLYPointStream< Real >::reset( void ) +{ + int fileType; + float version; + PlyProperty** plist; + if( _ply ) _free(); + _ply = ply_open_for_reading( _fileName, &_nr_elems, &_elist, &fileType, &version ); + if( !_ply ) + { + fprintf( stderr, "[ERROR] Failed to open ply file for reading: %s\n" , _fileName ); + exit( 0 ); + } + bool foundVertices = false; + for( int i=0 ; i<_nr_elems ; i++ ) + { + int num_elems; + int nr_props; + char* elem_name = _elist[i]; + plist = ply_get_element_description( _ply , elem_name , &num_elems , &nr_props ); + if( !plist ) + { + fprintf( stderr , "[ERROR] Failed to get element description: %s\n" , elem_name ); + exit( 0 ); + } + + if( equal_strings( "vertex" , elem_name ) ) + { + foundVertices = true; + _pCount = num_elems , _pIdx = 0; + for( int i=0 ; i::Components ; i++ ) + if( !ply_get_property( _ply , elem_name , &(PlyOrientedVertex< Real >::Properties[i]) ) ) + { + fprintf( stderr , "[ERROR] Failed to find property in ply file: %s\n" , PlyOrientedVertex< Real >::Properties[i].name ); + exit( 0 ); + } + } + for( int j=0 ; jname ); + free( plist[j] ); + } + free( plist ); + if( foundVertices ) break; + } + if( !foundVertices ) + { + fprintf( stderr , "[ERROR] Could not find vertices in ply file\n" ); + exit( 0 ); + } +} +template< class Real > +void PLYPointStream< Real >::_free( void ) +{ + if( _ply ) ply_close( _ply ) , _ply = NULL; + if( _elist ) + { + for( int i=0 ; i<_nr_elems ; i++ ) free( _elist[i] ); + free( _elist ); + } +} +template< class Real > +PLYPointStream< Real >::~PLYPointStream( void ) +{ + _free(); + if( _fileName ) delete[] _fileName , _fileName = NULL; +} +template< class Real > +bool PLYPointStream< Real >::nextPoint( Point3D< Real >& p , Point3D< Real >& n ) +{ + if( _pIdx<_pCount ) + { + PlyOrientedVertex< Real > op; + ply_get_element( _ply, (void *)&op ); + p = op.point; + n = op.normal; + _pIdx++; + return true; + } + else return false; +} diff --git a/src/plugins_experimental/filter_screened_poisson/Src/SparseMatrix.h b/src/plugins_experimental/filter_screened_poisson/Src/SparseMatrix.h new file mode 100755 index 000000000..85154ded6 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/SparseMatrix.h @@ -0,0 +1,218 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#ifndef __SPARSEMATRIX_HPP +#define __SPARSEMATRIX_HPP + +#define ZERO_TESTING_JACOBI 1 + + +#include "Vector.h" +#include "Array.h" + +template +struct MatrixEntry +{ + MatrixEntry( void ) { N =-1; Value = 0; } + MatrixEntry( int i ) { N = i; Value = 0; } + MatrixEntry( int i , T v ) { N = i; Value = v; } + int N; + T Value; +}; + +template class SparseMatrix +{ +private: + bool _contiguous; + int _maxEntriesPerRow; + void _init( void ); +public: + int rows; + Pointer( int ) rowSizes; + Pointer( Pointer( MatrixEntry< T > ) ) m_ppElements; + Pointer( MatrixEntry< T > ) operator[] ( int idx ) { return m_ppElements[idx]; } + ConstPointer( MatrixEntry< T > ) operator[] ( int idx ) const { return m_ppElements[idx]; } + + SparseMatrix( void ); + SparseMatrix( int rows ); + SparseMatrix( int rows , int maxEntriesPerRow ); + void Resize( int rows ); + void Resize( int rows , int maxEntriesPerRow ); + void SetRowSize( int row , int count ); + int Entries( void ) const; + + SparseMatrix( const SparseMatrix& M ); + ~SparseMatrix(); + + void SetZero(); + void SetIdentity(); + + SparseMatrix& operator = (const SparseMatrix& M); + + SparseMatrix operator * (const T& V) const; + SparseMatrix& operator *= (const T& V); + + + SparseMatrix operator * (const SparseMatrix& M) const; + SparseMatrix Multiply( const SparseMatrix& M ) const; + SparseMatrix MultiplyTranspose( const SparseMatrix& Mt ) const; + + template + Vector operator * (const Vector& V) const; + template + Vector Multiply( const Vector& V ) const; + template + void Multiply( const Vector& In , Vector& Out , int threads=1 ) const; + + + SparseMatrix Transpose() const; + + static int Solve (const SparseMatrix& M,const Vector& b, int iters,Vector& solution,const T eps=1e-8); + + template + static int SolveSymmetric( const SparseMatrix& M , const Vector& b , int iters , Vector& solution , const T2 eps=1e-8 , int reset=1 , int threads=1 ); + + bool write( FILE* fp ) const; + bool write( const char* fileName ) const; + bool read( FILE* fp ); + bool read( const char* fileName ); + + template< class T2 > + static int SolveJacobi( const SparseMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , Vector& Mx , T2 sor , int threads=1 , int offset=0 ); + template< class T2 > + static int SolveGS( const SparseMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , bool forward , int offset=0 ); + template< class T2 > + static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , bool forward , int threads=1 , int offset=0 ); + + template< class T2 > + static int SolveJacobi( const SparseMatrix& M , const Vector& b , Vector& x , Vector& Mx , T2 sor , int threads=1 , int offset=0 ); + template< class T2 > + static int SolveGS( const SparseMatrix& M , const Vector& b , Vector& x , bool forward , int offset=0 ); + template< class T2 > + static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix& M , const Vector& b , Vector& x , bool forward , int threads=1 , int offset=0 ); + + template< class T2 > + void getDiagonal( Vector< T2 >& diagonal , int threads=1 , int offset=0 ) const; +}; + + +template< class T2 > +struct MapReduceVector +{ +private: + int _dim; +public: + std::vector< T2* > out; + MapReduceVector( void ) { _dim = 0; } + ~MapReduceVector( void ) + { + if( _dim ) for( int t=0 ; t +class SparseSymmetricMatrix : public SparseMatrix< T > +{ +public: + + template< class T2 > + Vector< T2 > operator * ( const Vector& V ) const; + + template< class T2 > + Vector< T2 > Multiply( const Vector& V ) const; + + template< class T2 > + void Multiply( const Vector& In, Vector& Out , bool addDCTerm=false ) const; + + template< class T2 > + void Multiply( const Vector& In, Vector& Out , MapReduceVector< T2 >& OutScratch , bool addDCTerm=false ) const; + + template< class T2 > + void Multiply( const Vector& In, Vector& Out , std::vector< T2* >& OutScratch , const std::vector< int >& bounds ) const; + + template< class T2 > + static int SolveCG( const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , T2 eps=1e-8 , int reset=1 , int threads=0 , bool addDCTerm=false , bool solveNormal=false ); + + template< class T2 > + static int SolveCG( const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , MapReduceVector& scratch , T2 eps=1e-8 , int reset=1 , bool addDCTerm=false , bool solveNormal=false ); +#ifdef WIN32 + template< class T2 > + static int SolveCGAtomic( const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , T2 eps=1e-8 , int reset=1 , int threads=0 , bool solveNormal=false ); +#endif // WIN32 + template< class T2 > + static int SolveJacobi( const SparseSymmetricMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , MapReduceVector& scratch , Vector& Mx , T2 sor , int reset ); + template< class T2 > + static int SolveJacobi( const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , MapReduceVector& scratch , T2 sor=T2(1.) , int reset=1 ); + template< class T2 > + static int SolveJacobi( const SparseSymmetricMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , Vector& Mx , T2 sor , int reset ); + template< class T2 > + static int SolveJacobi( const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , T2 sor=T2(1.) , int reset=1 ); + + enum + { + ORDERING_UPPER_TRIANGULAR , + ORDERING_LOWER_TRIANGULAR , + ORDERING_NONE + }; + template< class T2 > + static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseSymmetricMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , MapReduceVector& scratch , Vector& Mx , Vector& dx , bool forward , int reset ); + template< class T2 > + static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , MapReduceVector& scratch , bool forward , int reset=1 ); + + template< class T2 > + static int SolveGS( const SparseSymmetricMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , MapReduceVector& scratch , Vector& Mx , Vector& dx , bool forward , int reset , int ordering ); + template< class T2 > + static int SolveGS( const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , MapReduceVector& scratch , bool forward , int reset=1 , int ordering=ORDERING_NONE ); + template< class T2 > + static int SolveGS( const SparseSymmetricMatrix& M , const Vector& diagonal , const Vector& b , Vector& x , Vector& Mx , Vector& dx , bool forward , int reset , int ordering ); + template< class T2 > + static int SolveGS( const SparseSymmetricMatrix& M , const Vector& b , int iters , Vector& x , bool forward , int reset=1 , int ordering=ORDERING_NONE ); + + template< class T2 > + void getDiagonal( Vector< T2 >& diagonal , int threads=1 ) const; +}; + +#include "SparseMatrix.inl" + +#endif + diff --git a/src/plugins_experimental/filter_screened_poisson/Src/SurfaceTrimmer.cpp b/src/plugins_experimental/filter_screened_poisson/Src/SurfaceTrimmer.cpp new file mode 100755 index 000000000..3821e6803 --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/SurfaceTrimmer.cpp @@ -0,0 +1,426 @@ +/* +Copyright (c) 2013, Michael Kazhdan +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#include +#include +#include +#include +#include +#include "CmdLineParser.h" +#include "Geometry.h" +#include "Ply.h" +#include "MAT.h" +#include "Time.h" + +#define FOR_RELEASE 1 + +cmdLineString In( "in" ) , Out( "out" ); +cmdLineInt Smooth( "smooth" , 5 ); +cmdLineFloat Trim( "trim" ) , IslandAreaRatio( "aRatio" , 0.001f ); +cmdLineFloatArray< 2 > ColorRange( "color" ); +cmdLineReadable PolygonMesh( "polygonMesh" ); + +cmdLineReadable* params[] = +{ + &In , &Out , &Trim , &PolygonMesh , &ColorRange , &Smooth , &IslandAreaRatio +}; + +void ShowUsage( char* ex ) +{ + printf( "Usage: %s\n" , ex ); + printf( "\t --%s \n" , In.name ); + printf( "\t[--%s ]\n" , Out.name ); + printf( "\t[--%s =%d]\n" , Smooth.name , Smooth.value ); + printf( "\t[--%s ]\n" , Trim.name ); + printf( "\t[--%s =%f]\n" , IslandAreaRatio.name , IslandAreaRatio.value ); + printf( "\t[--%s]\n" , PolygonMesh.name ); +#if !FOR_RELEASE + printf( "\t[--%s ]\n" , ColorRange.name ); +#endif // !FOR_RELEASE +} + +long long EdgeKey( int key1 , int key2 ) +{ + if( key1 +PlyValueVertex< Real > InterpolateVertices( const PlyValueVertex< Real >& v1 , const PlyValueVertex< Real >& v2 , const float& value ) +{ + if( v1.value==v2.value ) return (v1+v2)/Real(2.); + PlyValueVertex< Real > v; + + Real dx = (v1.value-value)/(v1.value-v2.value); + for( int i=0 ; i<3 ; i++ ) v.point.coords[i]=v1.point.coords[i]*(1.f-dx)+v2.point.coords[i]*dx; + v.value=v1.value*(1.f-dx)+v2.value*dx; + return v; +} + +template< class Real > +void ColorVertices( const std::vector< PlyValueVertex< Real > >& inVertices , std::vector< PlyColorVertex< Real > >& outVertices , float min , float max ) +{ + outVertices.resize( inVertices.size() ); + for( size_t i=0 ; i( 0.f , std::min< float >( 1.f , temp ) ); + temp *= 255; + outVertices[i].color[0] = outVertices[i].color[1] = outVertices[i].color[2] = (int)temp; + } +} +template< class Real > +void SmoothValues( std::vector< PlyValueVertex< Real > >& vertices , const std::vector< std::vector< int > >& polygons ) +{ + std::vector< int > count( vertices.size() ); + std::vector< Real > sums( vertices.size() , 0 ); + for( size_t i=0 ; i +void SmoothValues( std::vector< PlyValueVertex< Real > >& vertices , const std::vector< std::vector< int > >& polygons , Real min , Real max ) +{ + std::vector< int > count( vertices.size() ); + std::vector< Real > sums( vertices.size() , 0 ); + for( int i=0 ; imin && vertices[v1].valuemin && vertices[v2].value +void SplitPolygon + ( + const std::vector< int >& polygon , + std::vector< PlyValueVertex< Real > >& vertices , + std::vector< std::vector< int > >* ltPolygons , std::vector< std::vector< int > >* gtPolygons , + std::vector< bool >* ltFlags , std::vector< bool >* gtFlags , + hash_map< long long , int >& vertexTable , + Real trimValue + ) +{ + int sz = int( polygon.size() ); + std::vector< bool > gt( sz ); + int gtCount = 0; + for( int j=0 ; jtrimValue ); + if( gt[j] ) gtCount++; + } + if ( gtCount==sz ){ if( gtPolygons ) gtPolygons->push_back( polygon ) ; if( gtFlags ) gtFlags->push_back( false ); } + else if( gtCount==0 ){ if( ltPolygons ) ltPolygons->push_back( polygon ) ; if( ltFlags ) ltFlags->push_back( false ); } + else + { + int start; + for( start=0 ; start poly; + + // Add the initial vertex + { + int j1 = (start+int(sz)-1)%sz , j2 = start; + int v1 = polygon[j1] , v2 = polygon[j2]; + int vIdx; + hash_map< long long , int >::iterator iter = vertexTable.find( EdgeKey( v1 , v2 ) ); + if( iter==vertexTable.end() ) + { + vertexTable[ EdgeKey( v1 , v2 ) ] = vIdx = int( vertices.size() ); + vertices.push_back( InterpolateVertices( vertices[v1] , vertices[v2] , trimValue ) ); + } + else vIdx = iter->second; + poly.push_back( vIdx ); + } + + for( int _j=0 ; _j<=sz ; _j++ ) + { + int j1 = (_j+start+sz-1)%sz , j2 = (_j+start)%sz; + int v1 = polygon[j1] , v2 = polygon[j2]; + if( gt[j2]==gtFlag ) poly.push_back( v2 ); + else + { + int vIdx; + hash_map< long long , int >::iterator iter = vertexTable.find( EdgeKey( v1 , v2 ) ); + if( iter==vertexTable.end() ) + { + vertexTable[ EdgeKey( v1 , v2 ) ] = vIdx = int( vertices.size() ); + vertices.push_back( InterpolateVertices( vertices[v1] , vertices[v2] , trimValue ) ); + } + else vIdx = iter->second; + poly.push_back( vIdx ); + if( gtFlag ){ if( gtPolygons ) gtPolygons->push_back( poly ) ; if( ltFlags ) ltFlags->push_back( true ); } + else { if( ltPolygons ) ltPolygons->push_back( poly ) ; if( gtFlags ) gtFlags->push_back( true ); } + poly.clear() , poly.push_back( vIdx ) , poly.push_back( v2 ); + gtFlag = !gtFlag; + } + } + } +} +template< class Real > +void Triangulate( const std::vector< PlyValueVertex< Real > >& vertices , const std::vector< std::vector< int > >& polygons , std::vector< std::vector< int > >& triangles ) +{ + triangles.clear(); + for( size_t i=0 ; i3 ) + { + MinimalAreaTriangulation< Real > mat; + std::vector< Point3D< Real > > _vertices( polygons[i].size() ); + std::vector< TriangleIndex > _triangles; + for( int j=0 ; j +void RemoveHangingVertices( std::vector< Vertex >& vertices , std::vector< std::vector< int > >& polygons ) +{ + hash_map< int , int > vMap; + std::vector< bool > vertexFlags( vertices.size() , false ); + for( size_t i=0 ; i _vertices( vCount ); + for( int i=0 ; i >& polygons , std::vector< std::vector< int > >& components ) +{ + std::vector< int > polygonRoots( polygons.size() ); + for( size_t i=0 ; i edgeTable; + for( size_t i=0 ; i::iterator iter = edgeTable.find( eKey ); + if( iter==edgeTable.end() ) edgeTable[ eKey ] = int(i); + else + { + int p = iter->second; + while( polygonRoots[p]!=p ) + { + int temp = polygonRoots[p]; + polygonRoots[p] = int(i); + p = temp; + } + polygonRoots[p] = int(i); + } + } + } + for( size_t i=0 ; i vMap; + for( int i= 0 ; i +inline Point3D< Real > CrossProduct( Point3D< Real > p1 , Point3D< Real > p2 ){ return Point3D< Real >( p1[1]*p2[2]-p1[2]*p2[1] , p1[2]*p2[0]-p1[0]*p2[2] , p1[0]*p1[1]-p1[1]*p2[0] ); } +template< class Real > +double TriangleArea( Point3D< Real > v1 , Point3D< Real > v2 , Point3D< Real > v3 ) +{ + Point3D< Real > n = CrossProduct( v2-v1 , v3-v1 ); + return sqrt( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] ) / 2.; +} +template< class Real > +double PolygonArea( const std::vector< PlyValueVertex< Real > >& vertices , const std::vector< int >& polygon ) +{ + if( polygon.size()<3 ) return 0.; + else if( polygon.size()==3 ) return TriangleArea( vertices[polygon[0]].point , vertices[polygon[1]].point , vertices[polygon[2]].point ); + else + { + Point3D< Real > center; + for( size_t i=0 ; i > vertices; + std::vector< std::vector< int > > polygons; + + int ft , commentNum = paramNum+2; + char** comments; + bool readFlags[ PlyValueVertex< float >::Components ]; + PlyReadPolygons( In.value , vertices , polygons , PlyValueVertex< float >::Properties , PlyValueVertex< float >::Components , ft , &comments , &commentNum , readFlags ); + if( !readFlags[3] ){ fprintf( stderr , "[ERROR] vertices do not have value flag\n" ) ; return EXIT_FAILURE; } +#if 0 + if( Trim.set ) for( int i=0 ; i( min , vertices[i].value ) , max = std::max< float >( max , vertices[i].value ); + printf( "Value Range: [%f,%f]\n" , min , max ); + + + if( Trim.set ) + { + hash_map< long long , int > vertexTable; + std::vector< std::vector< int > > ltPolygons , gtPolygons; + std::vector< bool > ltFlags , gtFlags; + + for( int i=0 ; i0 ) + { + std::vector< std::vector< int > > _ltPolygons , _gtPolygons; + std::vector< std::vector< int > > ltComponents , gtComponents; + SetConnectedComponents( ltPolygons , ltComponents ); + SetConnectedComponents( gtPolygons , gtComponents ); + std::vector< double > ltAreas( ltComponents.size() , 0. ) , gtAreas( gtComponents.size() , 0. ); + std::vector< bool > ltComponentFlags( ltComponents.size() , false ) , gtComponentFlags( gtComponents.size() , false ); + double area = 0.; + for( size_t i=0 ; i > polys = ltPolygons; + Triangulate( vertices , ltPolygons , polys ) , ltPolygons = polys; + } + { + std::vector< std::vector< int > > polys = gtPolygons; + Triangulate( vertices , gtPolygons , polys ) , gtPolygons = polys; + } + } + + RemoveHangingVertices( vertices , gtPolygons ); + sprintf( comments[commentNum++] , "#Trimmed In: %9.1f (s)" , Time()-t ); + if( Out.set ) PlyWritePolygons( Out.value , vertices , gtPolygons , PlyValueVertex< float >::Properties , PlyValueVertex< float >::Components , ft , comments , commentNum ); + } + else + { + if( ColorRange.set ) min = ColorRange.values[0] , max = ColorRange.values[1]; + std::vector< PlyColorVertex< float > > outVertices; + ColorVertices( vertices , outVertices , min , max ); + if( Out.set ) PlyWritePolygons( Out.value , outVertices , polygons , PlyColorVertex< float >::Properties , PlyColorVertex< float >::Components , ft , comments , commentNum ); + } + + return EXIT_SUCCESS; +} + diff --git a/src/plugins_experimental/filter_screened_poisson/Src/Vector.inl b/src/plugins_experimental/filter_screened_poisson/Src/Vector.inl new file mode 100755 index 000000000..432f4341f --- /dev/null +++ b/src/plugins_experimental/filter_screened_poisson/Src/Vector.inl @@ -0,0 +1,302 @@ +/* +Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +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. + +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. + +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 +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 +TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR +BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. +*/ + +#ifndef __VECTORIMPL_HPP +#define __VECTORIMPL_HPP + +//////////// +// Vector // +//////////// +template +Vector::Vector( void ) +{ + m_N = 0; + m_pV = NullPointer< T >(); +} +template< class T > +Vector< T >::Vector( const Vector& V ) +{ + m_N = 0; + m_pV = NullPointer< T >(); + Resize( V.m_N ); + memcpy( m_pV , V.m_pV , m_N*sizeof(T) ); +} +template +Vector::Vector( size_t N ) +{ + m_N=0; + m_pV = NullPointer< T >(); + Resize(N); +} +template +void Vector::Resize( size_t N ) +{ + if( m_N!=N ) + { + if( m_N ) DeletePointer( m_pV ); + m_N = N; + m_pV = NewPointer< T >( N ); + } + if( N ) memset( m_pV , 0 , N*sizeof(T) ); +} + +template +Vector::Vector( size_t N, ConstPointer( T ) pV ) +{ + Resize(N); + memcpy( m_pV, pV, N*sizeof(T) ); +} +template +Vector::~Vector(){Resize(0);} +template +Vector& Vector::operator = (const Vector& V) +{ + Resize(V.m_N); + memcpy( m_pV, V.m_pV, m_N*sizeof(T) ); + return *this; +} +template +size_t Vector::Dimensions() const{return m_N;} +template +void Vector::SetZero(void){for (size_t i=0; i +const T& Vector::operator () (size_t i) const +{ + return m_pV[i]; +} +template +T& Vector::operator () (size_t i) +{ + return m_pV[i]; +} +template +const T& Vector::operator [] (size_t i) const +{ + return m_pV[i]; +} +template +T& Vector::operator [] (size_t i) +{ + return m_pV[i]; +} +template +Vector Vector::operator * (const T& A) const +{ + Vector V(*this); + for (size_t i=0; i +Vector& Vector::operator *= (const T& A) +{ + for (size_t i=0; i +Vector Vector::operator / (const T& A) const +{ + Vector V(*this); + for (size_t i=0; i +Vector& Vector::operator /= (const T& A) +{ + for (size_t i=0; i +Vector Vector::operator + (const T& A) const +{ + Vector V(*this); + for (size_t i=0; i +Vector& Vector::operator += (const T& A) +{ + for (size_t i=0; i +Vector Vector::operator - (const T& A) const +{ + Vector V(*this); + for (size_t i=0; i +Vector& Vector::operator -= (const T& A) +{ + for (size_t i=0; i +Vector Vector::operator + (const Vector& V0) const +{ + Vector V(m_N); + for (size_t i=0; i +Vector& Vector::operator += (const Vector& V) +{ + for ( size_t i=0 ; i +Vector Vector::operator - (const Vector& V0) const +{ + Vector V(m_N); + for (size_t i=0; i +Vector& Vector::operator -= (const Vector& V) +{ + for (size_t i=0; i +Vector Vector::operator - (void) const +{ + Vector V(m_N); + for (size_t i=0; i +Vector< T >& Vector< T >::Add( const Vector< T >* V , int count ) +{ + for( int c=0 ; c +Vector< T >& Vector< T >::AddScaled( const Vector& V , const T& scale ) +{ + for (size_t i=0; i +Vector& Vector::SubtractScaled(const Vector& V,const T& scale) +{ + for (size_t i=0; i +void Vector::Add( const Vector& V1 , const T& scale1 , const Vector& V2 , const T& scale2 , Vector& Out ) +{ + for( size_t i=0 ; i +void Vector::Add(const Vector& V1,const T& scale1,const Vector& V2,Vector& Out) +{ + for( size_t i=0 ; i +T Vector::Norm( size_t Ln ) const +{ + T N = T(); + for (size_t i = 0; i +void Vector::Normalize() +{ + T N = 1.0f/Norm(2); + for (size_t i = 0; i +T Vector< T >::Average( void ) const +{ + T N = T(); + for( size_t i=0 ; i +T Vector::Length() const +{ + T N = T(); + for (size_t i = 0; i +T Vector::Dot( const Vector& V ) const +{ + T V0 = T(); + for( size_t i=0 ; i +bool Vector< T >::read( const char* fileName ) +{ + FILE* fp = fopen( fileName , "rb" ); + if( !fp ) return false; + bool ret = read( fp ); + fclose( fp ); + return ret; +} +template< class T > +bool Vector< T >::write( const char* fileName ) const +{ + FILE* fp = fopen( fileName , "wb" ); + if( !fp ) return false; + bool ret = write( fp ); + fclose( fp ); + return ret; +} +template< class T > +bool Vector< T >::read( FILE* fp ) +{ + int d; + if( fread( &d , sizeof(int) , 1 , fp )!=1 ) return false; + Resize( d ); + if( fread( &(*this)[0] , sizeof( T ) , d , fp )!=d ) return false; + return true; +} +template< class T > +bool Vector< T >::write( FILE* fp ) const +{ + if( fwrite( &m_N , sizeof( int ) , 1 , fp )!=1 ) return false; + if( fwrite( &(*this)[0] , sizeof( T ) , m_N , fp )!=m_N ) return false; + return true; +} + + + +#endif