diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..70ee0e9
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,68 @@
+# /********************************************************************************************
+# * File: Makefile
+# * Author: $LastChangedBy: matthew $
+# * Revision: $Revision: 233 $
+# * Last Updated: $LastChangedDate: 2006-11-10 15:03:28 -0500 (Fri, 10 Nov 2006) $
+# ********************************************************************************************/
+
+PR_TARGET=PoissonRecon
+ST_TARGET=SurfaceTrimmer
+PR_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cpp PoissonRecon.cpp
+ST_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cpp SurfaceTrimmer.cpp
+
+CFLAGS += -fopenmp -Wno-deprecated
+LFLAGS += -lgomp
+
+CFLAGS_DEBUG = -DDEBUG -g3
+LFLAGS_DEBUG =
+
+CFLAGS_RELEASE = -O3 -DRELEASE -funroll-loops -ffast-math
+LFLAGS_RELEASE = -O3
+
+SRC = Src/
+BIN = Bin/Linux/
+INCLUDE = /usr/include/
+
+CC=gcc
+CXX=g++
+MD=mkdir
+
+PR_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(PR_SOURCE))))
+ST_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(ST_SOURCE))))
+
+
+all: CFLAGS += $(CFLAGS_DEBUG)
+all: LFLAGS += $(LFLAGS_DEBUG)
+all: $(BIN)
+all: $(BIN)$(PR_TARGET)
+all: $(BIN)$(ST_TARGET)
+
+release: CFLAGS += $(CFLAGS_RELEASE)
+release: LFLAGS += $(LFLAGS_RELEASE)
+release: $(BIN)
+release: $(BIN)$(PR_TARGET)
+release: $(BIN)$(ST_TARGET)
+
+clean:
+ rm -f $(BIN)$(PR_TARGET)
+ rm -f $(BIN)$(ST_TARGET)
+ rm -f $(PR_OBJECTS)
+ rm -f $(ST_OBJECTS)
+
+$(BIN):
+ $(MD) -p $(BIN)
+
+$(BIN)$(PR_TARGET): $(PR_OBJECTS)
+ $(CXX) -o $@ $(PR_OBJECTS) $(LFLAGS)
+
+$(BIN)$(ST_TARGET): $(ST_OBJECTS)
+ $(CXX) -o $@ $(ST_OBJECTS) $(LFLAGS)
+
+$(BIN)%.o: $(SRC)%.c
+ mkdir -p $(BIN)
+ $(CC) -c -o $@ $(CFLAGS) -I$(INCLUDE) $<
+
+$(BIN)%.o: $(SRC)%.cpp
+ mkdir -p $(BIN)
+ $(CXX) -c -o $@ $(CFLAGS) -I$(INCLUDE) $<
+
diff --git a/PoissonRecon.sln b/PoissonRecon.sln
new file mode 100644
index 0000000..2fc4c79
--- /dev/null
+++ b/PoissonRecon.sln
@@ -0,0 +1,45 @@
+
+Microsoft Visual Studio Solution File, Format Version 12.00
+# Visual Studio 2013
+VisualStudioVersion = 12.0.30501.0
+MinimumVisualStudioVersion = 10.0.40219.1
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "PoissonRecon", "PoissonRecon.vcxproj", "{46F87D0E-C53A-4F95-AB48-A5DBA8014340}"
+EndProject
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "SurfaceTrimmer", "SurfaceTrimmer.vcxproj", "{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}"
+EndProject
+Global
+ GlobalSection(SolutionConfigurationPlatforms) = preSolution
+ Debug|Win32 = Debug|Win32
+ Debug|x64 = Debug|x64
+ Release|Win32 = Release|Win32
+ Release|x64 = Release|x64
+ EndGlobalSection
+ GlobalSection(ProjectConfigurationPlatforms) = postSolution
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|Win32.ActiveCfg = Debug|Win32
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|Win32.Build.0 = Debug|Win32
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|x64.ActiveCfg = Debug|x64
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|x64.Build.0 = Debug|x64
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|Win32.ActiveCfg = Release|Win32
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|Win32.Build.0 = Release|Win32
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|x64.ActiveCfg = Release|x64
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|x64.Build.0 = Release|x64
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|Win32.ActiveCfg = Debug|Win32
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|Win32.Build.0 = Debug|Win32
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|x64.ActiveCfg = Debug|x64
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|x64.Build.0 = Debug|x64
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|Win32.ActiveCfg = Release|Win32
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|Win32.Build.0 = Release|Win32
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|x64.ActiveCfg = Release|x64
+ {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|x64.Build.0 = Release|x64
+ {7DB9957D-3644-4441-85A1-1A0928E087F0}.Debug|Win32.ActiveCfg = Debug|Win32
+ {7DB9957D-3644-4441-85A1-1A0928E087F0}.Debug|Win32.Build.0 = Debug|Win32
+ {7DB9957D-3644-4441-85A1-1A0928E087F0}.Debug|x64.ActiveCfg = Debug|x64
+ {7DB9957D-3644-4441-85A1-1A0928E087F0}.Debug|x64.Build.0 = Debug|x64
+ {7DB9957D-3644-4441-85A1-1A0928E087F0}.Release|Win32.ActiveCfg = Release|Win32
+ {7DB9957D-3644-4441-85A1-1A0928E087F0}.Release|Win32.Build.0 = Release|Win32
+ {7DB9957D-3644-4441-85A1-1A0928E087F0}.Release|x64.ActiveCfg = Release|x64
+ EndGlobalSection
+ GlobalSection(SolutionProperties) = preSolution
+ HideSolutionNode = FALSE
+ EndGlobalSection
+EndGlobal
diff --git a/PoissonRecon.vcxproj b/PoissonRecon.vcxproj
new file mode 100644
index 0000000..edbc0c2
--- /dev/null
+++ b/PoissonRecon.vcxproj
@@ -0,0 +1,232 @@
+
+
+
+
+ Debug
+ Win32
+
+
+ Debug
+ x64
+
+
+ Release
+ Win32
+
+
+ Release
+ x64
+
+
+
+ PoissonRecon
+ {46F87D0E-C53A-4F95-AB48-A5DBA8014340}
+ PoissonRecon
+ Win32Proj
+
+
+ Application
+ MultiByte
+ true
+ v120
+
+
+ Application
+ MultiByte
+ v120
+
+
+ Application
+ MultiByte
+ true
+ v120
+
+
+ Application
+ MultiByte
+ v120
+
+
+
+
+
+
+
+
+
+
+
+
+ <_ProjectFileVersion>10.0.30319.1
+ $(SolutionDir)\Bin\$(Platform)\$(Configuration)\
+ $(SolutionDir)\Obj\$(Platform)\$(Configuration)\
+ true
+ $(SolutionDir)\Bin\$(Platform)\$(Configuration)\
+ $(SolutionDir)\Obj\$(Platform)\$(Configuration)\
+ true
+ $(SolutionDir)Bin\$(Platform)\$(Configuration)\
+ $(SolutionDir)\Obj\$(TargetName)\$(Platform)\$(Configuration)\
+ false
+ $(SolutionDir)\Bin\$(Platform)\$(Configuration)\
+ $(SolutionDir)\Obj\$(TargetName)\$(Platform)\$(Configuration)\
+ false
+ .exe
+
+
+
+ Disabled
+ WIN32;_DEBUG;_CONSOLE;_CRT_SECURE_NO_DEPRECATE;%(PreprocessorDefinitions)
+ true
+ EnableFastChecks
+ MultiThreadedDebugDLL
+
+
+ Level3
+ EditAndContinue
+
+ true
+ Console
+ false
+
+
+ MachineX86
+
+
+
+
+ X64
+
+
+ Disabled
+ WIN32;_DEBUG;_CONSOLE;_CRT_SECURE_NO_DEPRECATE;%(PreprocessorDefinitions)
+ true
+ EnableFastChecks
+ MultiThreadedDebugDLL
+
+
+ Level3
+ ProgramDatabase
+
+ true
+ Console
+ false
+
+
+ MachineX64
+
+
+
+
+ WIN32;NDEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;_CRT_SECURE_NO_DEPRECATE;%(PreprocessorDefinitions)
+ MultiThreadedDLL
+
+
+ Level3
+ ProgramDatabase
+ true
+ %(AdditionalIncludeDirectories)
+
+ true
+ Console
+ true
+ true
+ true
+ false
+
+
+ MachineX86
+ psapi.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)
+
+
+
+
+ X64
+
+
+ %(AdditionalIncludeDirectories)
+ WIN32;NDEBUG;_CONSOLE;_CRT_SECURE_NO_DEPRECATE;%(PreprocessorDefinitions)
+ MultiThreadedDLL
+
+
+ Level3
+ ProgramDatabase
+ Precise
+ true
+ false
+ AdvancedVectorExtensions2
+
+ true
+ Console
+ true
+ true
+ false
+
+
+ MachineX64
+ psapi.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)
+
+
+ $(OutDir)$(TargetName)$(TargetExt)
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/PoissonRecon.vcxproj.filters b/PoissonRecon.vcxproj.filters
new file mode 100644
index 0000000..852bbac
--- /dev/null
+++ b/PoissonRecon.vcxproj.filters
@@ -0,0 +1,143 @@
+
+
+
+
+ {4FC737F1-C7A5-4376-A066-2A32D752A2FF}
+ cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx
+
+
+ {93995380-89BD-4b04-88EB-625FBE52EBFB}
+ h;hpp;hxx;hm;inl;inc;xsd
+
+
+ {67DA6AB6-F800-4c08-8B7A-83BB121AAD01}
+ inc;inl
+
+
+
+
+ Source Files
+
+
+ Source Files
+
+
+ Source Files
+
+
+ Source Files
+
+
+ Source Files
+
+
+ Source Files
+
+
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+ Header Files
+
+
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+ Include Files
+
+
+
\ No newline at end of file
diff --git a/README.md b/README.md
index 0e6e59a..0578f41 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,290 @@
-# PoissonRecon
-Poisson Surface Reconstruction
+Screened Poisson Surface Reconstruction (Version 6.13a)
+links
+executables
+usage
+changes
+
+LINKS SGP 2006 ,
+ToG 2013
+Win32 , Win64 )V6.13 ,
+V6.12 ,
+V6.11 ,
+V6.1 ,
+V6 ,
+V5.71 ,
+V5.6 ,
+V5.5a ,
+V5.1 ,
+V5 ,
+V4.51 ,
+V4.5 ,
+V4 ,
+V3 ,
+V2 ,
+V1 )
+
+PoissonRecon:
+
+--in <input points >
+ This string is the name of the file from which the point set will be read..ply , the file should be in
+PLY format, giving the list of oriented
+vertices with the x-, y-, and z-coordinates of the positions encoded by the properties x , y , and
+z and the x-, y-, and z-coordinates of the normals encoded by the properties nx , ny , and
+nz ..bnpts , the file should be a binary file, consisting of blocks of 6 32-bit
+floats: x-, y-, and z-coordinates of the point's position, followed by the x-, y-, and z-coordinates
+of the point's normal. (No information about the number of oriented point samples should be specified.) [--out <output triangle mesh >]
+ This string is the name of the file to which the triangle mesh will be written.
+The file is written in PLY format.
+
+ [--voxel <output voxel grid >]
+ This string is the name of the file to which the sampled implicit function will be written.
+The filw is wrtten out in binary, with the first 4 bytes corresponding to the (integer) sampling resolution, 2^d ,
+and the next 4 x 2^d x 2^d x 2^d bytes corresponding to the (single precision) floating point values
+of the implicit function.
+
+ [--depth <reconstruction depth >]
+ This integer is the maximum depth of the tree that will be used for surface reconstruction.
+Running at depth d corresponds to solving on a voxel grid whose resolution is no larger than
+2^d x 2^d x 2^d . Note that since the reconstructor adapts the octree to the
+sampling density, the specified reconstruction depth is only an upper bound. [--fullDepth <adaptive octree depth >]
+ This integer specifies the depth beyond depth the octree will be adapted.
+At coarser depths, the octree will be complete, containing all
+2^d x 2^d x 2^d nodes. [--voxelDepth <voxel sampling depth >]
+ This integer is the depth of the regular grid over which the implicit function is to be sampled.
+Running at depth d corresponds to sampling on a voxel grid whose resolution is 2^d x 2^d x 2^d .--depth parameter.
+
+ [--cgDepth <conjugate gradients solver depth >]
+ This integer is the depth up to which a conjugate-gradients solver will be used to solve the linear system. Beyond this depth Gauss-Seidel relaxation will be used. [--scale <scale factor >]
+ This floating point value specifies the ratio between the diameter of the cube used for reconstruction
+and the diameter of the samples' bounding cube. [--samplesPerNode <minimum number of samples >]
+ This floating point value specifies the minimum number of sample points that should fall within an
+octree node as the octree construction is adapted to sampling density. For noise-free samples, small values
+in the range [1.0 - 5.0] can be used. For more noisy samples, larger values in the range [15.0 - 20.0] may
+be needed to provide a smoother, noise-reduced, reconstruction. [--pointWeight <interpolation weight >]
+ This floating point value specifies the importants that interpolation of the point samples
+is given in the formulation of the screened Poisson equation. [--iters <GS iters >]
+ This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy. [--threads <number of processing threads >]
+ This integer specifies the number of threads across which the reconstruction
+algorithm should be parallelized. [--confidence ]
+ Enabling this flag tells the reconstructor to use the size of the normals as confidence information. When the flag
+is not enabled, all normals are normalized to have unit-length prior to reconstruction.
+
+ [--nWeights ]
+ Enabling this flag tells the reconstructor to use the size of the normals to modulate the interpolation weights. When the flag
+is not enabled, all points are given the same weight.
+
+ [--polygonMesh ]
+ Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes).
+
+ [--density ]
+ Enabling this flag tells the reconstructor to output the estimated depth values of the iso-surface vertices.
+
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of
+individual components of the surface reconstructor.
+
+
+SurfaceTrimmer:
+
+--in <input triangle mesh >
+ This string is the name of the file from which the triangle mesh will be read.
+The file is read in PLY format and it is assumed that the vertices have a value field which stores the signal's value. (When run with --density flag, the reconstructor will output this field with the mesh vertices.)
+
+ --trim <trimming value >
+ This floating point values specifies the value for mesh trimming. The subset of the mesh with signal value less than the trim value is discarded.
+
+ [--out <output triangle mesh >]
+ This string is the name of the file to which the triangle mesh will be written.
+The file is written in PLY format.
+
+ [--smooth <smoothing iterations >]
+ This integer values the number of umbrella smoothing operations to perform on the signal before trimming. [--aRatio <island area ratio >]
+ This floating point value specifies the area ratio that defines a disconnected component as an "island". Connected components whose area, relative to the total area of the mesh, are smaller than this value will be merged into the output surface to close small holes, and will be discarded from the output surface to remove small disconnected components. [--polygonMesh ]
+ Enabling this flag tells the trimmer to output a polygon mesh (rather than triangulating the trimming results).
+
+
+
+ USAGE
+
+ Bunny range scans . The orientation of the sample points was estimated
+using the connectivity information within individual range scans.% PoissonRecon --in bunny.points.ply --out bunny.unscreened.ply --depth 10 --pointWeight 0--pointWeight 0 argument to disable the screening.% PoissonRecon --in bunny.points.ply --out bunny.screened.ply --depth 10% PoissonRecon --in bunny.points.ply --out bunny.screened.ply --depth 10 --density% SurfaceTrimmer --in bunny.screened.ply --out bunny.screened.trimmed.ply --trim 7 --aRatio 0% SurfaceTrimmer --in bunny.screened.ply --out bunny.screened.trimmed.ply --trim 7 Horse % PoissonRecon --in horse.npts --out horse.ply --depth 10
+
+To convert the binary PLY format to
+Hugues Hoppe's ASCII
+mesh format, a Perl script is provided.% ply2mesh.pl bunny.ply > bunny.mCHANGES Version 3 :
+
+ The implementation of the --samplesPerNode parameter has been modified so that a value of "1" more closely corresponds to a distribution with one sample per leaf node.
+ The code has been modified to support compilation under MSVC 2010 and the associated solution and project files are now provided. (Due to a bug in the Visual Studios compiler, this required modifying the implementation of some of the bit-shifting operators.)
+
+Version 4 :
+
+ The code supports screened reconstruction, with interpolation weight specified through the --pointWeight parameter.
+ The code has been implemented to support parallel processing, with the number of threads used for parallelization specified by the --threads parameter.
+ The input point set can now also be in PLY format, and the file-type is determined by the extension, so that the --binary flag is now obsolete.
+ At depths coarser than the one specified by the value --minDepth the octree is no longer adaptive but rather complete, simplifying the prolongation operator.
+
+Version 4.5 :
+
+ The algorithmic complexity of the solver was reduced from log-linear to linear.
+
+Version 4.51 :
+
+ Smart pointers were added to ensure that memory accesses were in bounds.
+
+Version 5 :
+
+ The --density flag was added to the reconstructor to output the estimated depth of the iso-vertices.
+ The SurfaceTrimmer executable was added to support trimming off the subset of the reconstructed surface that are far away from the input samples, thereby allowing for the generation of non-water-tight surface.
+
+
+Version 5.1 :
+
+ Minor bug-fix to address incorrect neighborhood estimation in the octree finalization.
+
+
+Version 5.5a :
+
+ Modified to support depths greater than 14. (Should work up to 18 or 19 now.)
+ Improved speed and memory performance by removing the construction of integral and value tables.
+ Fixed a bug in Version 5.5 that used memory and took more time without doing anything useful.
+
+
+Version 5.6 :
+
+ Added the --normalWeight flag to support setting a point's interpolation weight in proportion to the magnitude of its normal.
+
+
+Version 5.7 :
+
+ Modified the setting of the constraints, replacing the map/reduce implementation with OpenMP atomics to reduce memory usage.
+ Fixed bugs that caused numerical overflow when processing large point clouds on multi-core machines.
+ Improved efficiency of the iso-surface extraction phse.
+
+
+Version 5.71 :
+
+ Added the function GetSolutionValue to support the evaluation of the implicit function at a specific point.
+
+
+Version 6 :
+
+ Modified the solver to use Gauss-Seidel relaxation instead of conjugate-gradients at finer resolution.
+ Re-ordered the implementation of the solver so that only a windowed subset of the matrix is in memory at any time, thereby reducing the memory usage during the solver phase.
+ Separated the storage of the data associated with the octree nodes from the topology.
+
+
+Version 6.1 :
+
+ Re-ordered the implementation of the iso-surface extraction so that only a windowed subset of the octree is in memory at any time, thereby reducing the memory usage during the extracted phase.
+
+
+Version 6.11 :
+
+ Fixed a bug that created a crash in the evaluation phase when --pointWeight is set zero.
+
+
+Version 6.12 :
+
+ Removed the OpenMP firstprivate directive as it seemed to cause trouble under Linux compilations.
+
+
+Version 6.13 :
+
+ Added a MemoryPointStream class in PointStream.inl to support in-memory point clouds.
+ Modified the signature of Octree::SetTree in MultiGridOctreeData.h to take in a pointer to an object of type PointStream rather than a file-name.
+
+
+Version 6.13a :
+
+ Modified the signature of Octree::SetIsoSurface to rerun a void . [cloudcompare ]
+ Added a definition of SetIsoVertexValue supporting double precision vertices. [cloudcompare ]
+ Removed Time.[h/cpp] from the repository. [cloudcompare /asmaloney ]
+ Fixed assignment bug in Octree::SetSliceIsoVertices . [asmaloney ]
+ Fixed initialization bug in SortedTreeNodes::SliceTableData and SortedTreeNodes::XSliceTableData . [asmaloney ]
+ Included stdlib.h in Geometry.h . [asmaloney ]
+ Fixed default value bug in declaration of Octree::SetTree . [asmaloney ]
+
+
+HOME
diff --git a/Src/Allocator.h b/Src/Allocator.h
new file mode 100644
index 0000000..e2723b9
--- /dev/null
+++ b/Src/Allocator.h
@@ -0,0 +1,163 @@
+/*
+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 ALLOCATOR_INCLUDED
+#define ALLOCATOR_INCLUDED
+#include
+
+class AllocatorState{
+public:
+ int index,remains;
+};
+/** This templated class assists in memory allocation and is well suited for instances
+ * when it is known that the sequence of memory allocations is performed in a stack-based
+ * manner, so that memory allocated last is released first. It also preallocates memory
+ * in chunks so that multiple requests for small chunks of memory do not require separate
+ * system calls to the memory manager.
+ * The allocator is templated off of the class of objects that we would like it to allocate,
+ * ensuring that appropriate constructors and destructors are called as necessary.
+ */
+template
+class Allocator
+{
+ int blockSize;
+ int index , remains;
+ std::vector< T* > memory;
+public:
+ Allocator( void ){ blockSize = index = remains = 0; }
+ ~Allocator( void ){ reset(); }
+
+ /** This method is the allocators destructor. It frees up any of the memory that
+ * it has allocated. */
+ void reset( void )
+ {
+ for(size_t i=0;iblockSize = blockSize;
+ index=-1;
+ remains=0;
+ }
+
+ /** This method returns a pointer to an array of elements objects. If there is left over pre-allocated
+ * memory, this method simply returns a pointer to the next free piece of memory, otherwise it pre-allocates
+ * more memory. Note that if the number of objects requested is larger than the value blockSize with which
+ * the allocator was initialized, the request for memory will fail.
+ */
+ T* newElements( int elements=1 )
+ {
+ T* mem;
+ if( !elements ) return NULL;
+ if( elements>blockSize ) fprintf( stderr , "[ERROR] Allocator: elements bigger than block-size: %d>%d\n" , elements , blockSize ) , exit( 0 );
+ if( remains
+
+#define ARRAY_DEBUG 0
+#ifdef _WIN64
+#define ASSERT( x ) { if( !( x ) ) __debugbreak(); }
+#else // !_WIN64
+#ifdef _WIN32
+#define ASSERT( x ) { if( !( x ) ) _asm{ int 0x03 } }
+#else // !_WIN32
+#define ASSERT( x ) { if( !( x ) ) exit(0); }
+#endif // _WIN32
+#endif // _WIN64
+
+// Code from http://stackoverflow.com
+void* aligned_malloc( size_t size , size_t align )
+{
+ // Align enough for the data, the alignment padding, and room to store a pointer to the actual start of the memory
+ void* mem = malloc( size + align + sizeof( void* ) );
+ // The position at which we could potentially start addressing
+ char* amem = ( (char*)mem ) + sizeof( void* );
+ // Add align-1 to the start of the address and then zero out at most of the first align-1 bits.
+ amem = ( char* )( ( (size_t)( ( (char*)amem ) + (align-1) ) ) & ~( align-1 ) );
+ // Pre-write the actual address
+ ( ( void** ) amem )[-1] = mem;
+ return amem;
+}
+void aligned_free( void* mem ) { free( ( ( void** )mem )[-1] ); }
+
+#if ARRAY_DEBUG
+#pragma message ( "[WARNING] Array debugging is enabled" )
+#include "Array.inl"
+#define Pointer( ... ) Array< __VA_ARGS__ >
+#define ConstPointer( ... ) ConstArray< __VA_ARGS__ >
+template< class C > void FreePointer( Array< C >& a ){ a.Free( ); }
+template< class C > void AlignedFreePointer( Array< C >& a ){ a.Free( ); }
+template< class C > void VFreePointer( Array< C >& a ){ a.Free( ); }
+template< class C > void DeletePointer( Array< C >& a ){ a.Delete( ); }
+
+template< class C > Array< C > NewPointer( size_t size , const char* name=NULL ){ return Array< C >::New ( size , name ); }
+template< class C > Array< C > AllocPointer( size_t size , const char* name=NULL ){ return Array< C >::Alloc ( size , false , name ); }
+template< class C > Array< C > AlignedAllocPointer( size_t size , size_t alignment , const char* name=NULL ){ return Array< C >::AlignedAlloc( size , alignment , false , name ); }
+template< class C > Array< C > ReAllocPointer( Array< C >& a , size_t size , const char* name=NULL ){ return Array< C >::ReAlloc ( a , size , false , name ); }
+
+template< class C > Array< C > NullPointer( void ){ return Array< C >( ); }
+
+template< class C > C* PointerAddress( Array< C >& a ) { return a.pointer(); }
+template< class C > const C* PointerAddress( ConstArray< C >& a ) { return a.pointer(); }
+template< class C > Array< C > GetPointer( C& c ) { return Array< C >::FromPointer( &c , 1 ); }
+template< class C > ConstArray< C > GetPointer( const C& c ) { return ConstArray< C >::FromPointer( &c , 1 ); }
+template< class C > Array< C > GetPointer( std::vector< C >& v ){ return Array< C >::FromPointer( &v[0] , v.size() ); }
+template< class C > ConstArray< C > GetPointer( const std::vector< C >& v ){ return ConstArray< C >::FromPointer( &v[0] , v.size() ); }
+
+#else // !ARRAY_DEBUG
+#define Pointer( ... ) __VA_ARGS__*
+#define ConstPointer( ... ) const __VA_ARGS__*
+
+#define FreePointer( ... ) { if( __VA_ARGS__ ) free( __VA_ARGS__ ) , __VA_ARGS__ = NULL; }
+#define AlignedFreePointer( ... ) { if( __VA_ARGS__ ) aligned_free( __VA_ARGS__ ) , __VA_ARGS__ = NULL; }
+#define DeletePointer( ... ) { if( __VA_ARGS__ ) delete[] __VA_ARGS__ , __VA_ARGS__ = NULL; }
+
+template< class C > C* NewPointer( size_t size , const char* name=NULL ){ return new C[size]; }
+template< class C > C* AllocPointer( size_t size , const char* name=NULL ){ return (C*) malloc( sizeof(C) * size ); }
+template< class C > C* AlignedAllocPointer( size_t size , size_t alignment , const char* name=NULL ){ return (C*)aligned_malloc( sizeof(C) * size , alignment ); }
+template< class C > C* ReAllocPointer( C* c , size_t size , const char* name=NULL ){ return (C*) realloc( c , sizeof(C) * size ); }
+
+template< class C > C* NullPointer( void ){ return NULL; }
+
+template< class C > C* PointerAddress( C* c ){ return c; }
+template< class C > const C* PointerAddress( const C* c ){ return c; }
+template< class C > C* GetPointer( C& c ){ return &c; }
+template< class C > const C* GetPointer( const C& c ){ return &c; }
+template< class C > C* GetPointer( std::vector< C >& v ){ return &v[0]; }
+template< class C > const C* GetPointer( const std::vector< C >& v ){ return &v[0]; }
+#endif // ARRAY_DEBUG
+#endif // ARRAY_INCLUDED
diff --git a/Src/Array.inl b/Src/Array.inl
new file mode 100644
index 0000000..511b9a0
--- /dev/null
+++ b/Src/Array.inl
@@ -0,0 +1,658 @@
+/*
+Copyright (c) 2011, Michael Kazhdan and Ming Chuang
+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.
+*/
+#define FULL_ARRAY_DEBUG 0 // Note that this is not thread-safe
+
+#include
+#include
+#include
+#ifdef _WIN32
+#include
+#endif // _WIN32
+#include
+
+inline bool isfinitef( float fp ){ float f=fp; return ((*(unsigned *)&f)&0x7f800000)!=0x7f800000; }
+
+
+template< class C > bool IsValid( const C& c );
+#if _DEBUG
+template< > inline bool IsValid< float >( const float& f ) { return isfinitef( f ) && ( f==0.f || abs(f)>1e-31f ); }
+#else // !_DEBUG
+template< > inline bool IsValid< float >( const float& f ) { return isfinitef( f ); }
+#endif // _DEBUG
+template< > inline bool IsValid< __m128 >( const __m128& m )
+{
+ const __m128* addr = &m;
+ if( size_t(addr) & 15 ) return false;
+ else return true;
+}
+template< class C > inline bool IsValid( const C& c ){ return true; }
+
+
+#if FULL_ARRAY_DEBUG
+class DebugMemoryInfo
+{
+public:
+ const void* address;
+ char name[512];
+};
+static std::vector< DebugMemoryInfo > memoryInfo;
+#endif // FULL_ARRAY_DEBUG
+
+template< class C >
+class Array
+{
+ void _assertBounds( long long idx ) const
+ {
+ if( idx=max )
+ {
+ fprintf( stderr , "Array index out-of-bounds: %lld <= %lld < %lld\n" , min , idx , max );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+protected:
+ C *data , *_data;
+ long long min , max;
+#if FULL_ARRAY_DEBUG
+ static void _AddMemoryInfo( const void* ptr , const char* name )
+ {
+ size_t sz = memoryInfo.size();
+ memoryInfo.resize( sz + 1 );
+ memoryInfo[sz].address = ptr;
+ if( name ) strcpy( memoryInfo[sz].name , name );
+ else memoryInfo[sz].name[0] = 0;
+ }
+ static void _RemoveMemoryInfo( const void* ptr )
+ {
+ {
+ size_t idx;
+ for( idx=0 ; idx
+ Array( Array< D >& a )
+ {
+ _data = NULL;
+ if( !a )
+ {
+ data = NULL;
+ min = max = 0;
+ }
+ else
+ {
+ // [WARNING] Chaning szC and szD to size_t causes some really strange behavior.
+ long long szC = sizeof( C );
+ long long szD = sizeof( D );
+ data = (C*)&a[0];
+ min = ( a.minimum() * szD ) / szC;
+ max = ( a.maximum() * szD ) / szC;
+ if( min*szC!=a.minimum()*szD || max*szC!=a.maximum()*szD )
+ {
+ fprintf( stderr , "Could not convert array [ %lld , %lld ] * %lld => [ %lld , %lld ] * %lld\n" , a.minimum() , a.maximum() , szD , min , max , szC );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+ }
+ static Array FromPointer( C* data , long long max )
+ {
+ Array a;
+ a._data = NULL;
+ a.data = data;
+ a.min = 0;
+ a.max = max;
+ return a;
+ }
+ static Array FromPointer( C* data , long long min , long long max )
+ {
+ Array a;
+ a._data = NULL;
+ a.data = data;
+ a.min = min;
+ a.max = max;
+ return a;
+ }
+ inline bool operator == ( const Array< C >& a ) const { return data==a.data; }
+ inline bool operator != ( const Array< C >& a ) const { return data!=a.data; }
+ inline bool operator == ( const C* c ) const { return data==c; }
+ inline bool operator != ( const C* c ) const { return data!=c; }
+ inline C* operator -> ( void )
+ {
+ _assertBounds( 0 );
+ return data;
+ }
+ inline const C* operator -> ( ) const
+ {
+ _assertBounds( 0 );
+ return data;
+ }
+ inline C& operator[]( long long idx )
+ {
+ _assertBounds( idx );
+ return data[idx];
+ }
+ inline const C& operator[]( long long idx ) const
+ {
+ _assertBounds( idx );
+ return data[idx];
+ }
+ inline Array operator + ( int idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array operator + ( long long idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array operator + ( unsigned int idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array operator + ( unsigned long long idx ) const
+ {
+ Array a;
+ a._data = _data;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline Array& operator += ( int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator += ( long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator += ( unsigned int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator += ( unsigned long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline Array& operator ++ ( void ) { return (*this) += 1; }
+ Array operator - ( int idx ) const { return (*this) + (-idx); }
+ Array operator - ( long long idx ) const { return (*this) + (-idx); }
+ Array operator - ( unsigned int idx ) const { return (*this) + (-idx); }
+ Array operator - ( unsigned long long idx ) const { return (*this) + (-idx); }
+ Array& operator -= ( int idx ) { return (*this) += (-idx); }
+ Array& operator -= ( long long idx ) { return (*this) += (-idx); }
+ Array& operator -= ( unsigned int idx ) { return (*this) += (-idx); }
+ Array& operator -= ( unsigned long long idx ) { return (*this) += (-idx); }
+ Array& operator -- ( void ) { return (*this) -= 1; }
+ long long operator - ( const Array& a ) const { return ( long long )( data - a.data ); }
+
+ void Free( void )
+ {
+ if( _data )
+ {
+ free( _data );
+#if FULL_ARRAY_DEBUG
+ _RemoveMemoryInfo( _data );
+#endif // FULL_ARRAY_DEBUG
+ }
+ (*this) = Array( );
+ }
+ void Delete( void )
+ {
+ if( _data )
+ {
+ delete[] _data;
+#if FULL_ARRAY_DEBUG
+ _RemoveMemoryInfo( _data );
+#endif // FULL_ARRAY_DEBUG
+ }
+ (*this) = Array( );
+ }
+ C* pointer( void ){ return data; }
+ const C* pointer( void ) const { return data; }
+ bool operator !( void ) const { return data==NULL; }
+ operator bool( ) const { return data!=NULL; }
+};
+
+template< class C >
+class ConstArray
+{
+ void _assertBounds( long long idx ) const
+ {
+ if( idx=max )
+ {
+ fprintf( stderr , "ConstArray index out-of-bounds: %lld <= %lld < %lld\n" , min , idx , max );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+protected:
+ const C *data;
+ long long min , max;
+public:
+ long long minimum( void ) const { return min; }
+ long long maximum( void ) const { return max; }
+
+ inline ConstArray( void )
+ {
+ data = NULL;
+ min = max = 0;
+ }
+ inline ConstArray( const Array< C >& a )
+ {
+ // [WARNING] Changing szC and szD to size_t causes some really strange behavior.
+ data = ( const C* )a.pointer( );
+ min = a.minimum();
+ max = a.maximum();
+ }
+ template< class D >
+ inline ConstArray( const Array< D >& a )
+ {
+ // [WARNING] Changing szC and szD to size_t causes some really strange behavior.
+ long long szC = ( long long ) sizeof( C );
+ long long szD = ( long long ) sizeof( D );
+ data = ( const C* )a.pointer( );
+ min = ( a.minimum() * szD ) / szC;
+ max = ( a.maximum() * szD ) / szC;
+ if( min*szC!=a.minimum()*szD || max*szC!=a.maximum()*szD )
+ {
+// fprintf( stderr , "Could not convert const array [ %lld , %lld ] * %lld => [ %lld , %lld ] * %lld\n" , a.minimum() , a.maximum() , szD , min , max , szC );
+ fprintf( stderr , "Could not convert const array [ %lld , %lld ] * %lld => [ %lld , %lld ] * %lld\n %lld %lld %lld\n" , a.minimum() , a.maximum() , szD , min , max , szC , a.minimum() , a.minimum()*szD , (a.minimum()*szD)/szC );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+ template< class D >
+ inline ConstArray( const ConstArray< D >& a )
+ {
+ // [WARNING] Chaning szC and szD to size_t causes some really strange behavior.
+ long long szC = sizeof( C );
+ long long szD = sizeof( D );
+ data = ( const C*)a.pointer( );
+ min = ( a.minimum() * szD ) / szC;
+ max = ( a.maximum() * szD ) / szC;
+ if( min*szC!=a.minimum()*szD || max*szC!=a.maximum()*szD )
+ {
+ fprintf( stderr , "Could not convert array [ %lld , %lld ] * %lld => [ %lld , %lld ] * %lld\n" , a.minimum() , a.maximum() , szD , min , max , szC );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ }
+ static ConstArray FromPointer( const C* data , long long max )
+ {
+ ConstArray a;
+ a.data = data;
+ a.min = 0;
+ a.max = max;
+ return a;
+ }
+ static ConstArray FromPointer( const C* data , long long min , long long max )
+ {
+ ConstArray a;
+ a.data = data;
+ a.min = min;
+ a.max = max;
+ return a;
+ }
+
+ inline bool operator == ( const ConstArray< C >& a ) const { return data==a.data; }
+ inline bool operator != ( const ConstArray< C >& a ) const { return data!=a.data; }
+ inline bool operator == ( const C* c ) const { return data==c; }
+ inline bool operator != ( const C* c ) const { return data!=c; }
+ inline const C* operator -> ( void )
+ {
+ _assertBounds( 0 );
+ return data;
+ }
+ inline const C& operator[]( long long idx ) const
+ {
+ _assertBounds( idx );
+ return data[idx];
+ }
+ inline ConstArray operator + ( int idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray operator + ( long long idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray operator + ( unsigned int idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray operator + ( unsigned long long idx ) const
+ {
+ ConstArray a;
+ a.data = data+idx;
+ a.min = min-idx;
+ a.max = max-idx;
+ return a;
+ }
+ inline ConstArray& operator += ( int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator += ( long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator += ( unsigned int idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator += ( unsigned long long idx )
+ {
+ min -= idx;
+ max -= idx;
+ data += idx;
+ return (*this);
+ }
+ inline ConstArray& operator ++ ( void ) { return (*this) += 1; }
+ ConstArray operator - ( int idx ) const { return (*this) + (-idx); }
+ ConstArray operator - ( long long idx ) const { return (*this) + (-idx); }
+ ConstArray operator - ( unsigned int idx ) const { return (*this) + (-idx); }
+ ConstArray operator - ( unsigned long long idx ) const { return (*this) + (-idx); }
+ ConstArray& operator -= ( int idx ) { return (*this) += (-idx); }
+ ConstArray& operator -= ( long long idx ) { return (*this) += (-idx); }
+ ConstArray& operator -= ( unsigned int idx ) { return (*this) += (-idx); }
+ ConstArray& operator -= ( unsigned long long idx ) { return (*this) += (-idx); }
+ ConstArray& operator -- ( void ) { return (*this) -= 1; }
+ long long operator - ( const ConstArray& a ) const { return ( long long )( data - a.data ); }
+ long long operator - ( const Array< C >& a ) const { return ( long long )( data - a.pointer() ); }
+
+ const C* pointer( void ) const { return data; }
+ bool operator !( void ) { return data==NULL; }
+ operator bool( ) { return data!=NULL; }
+};
+
+#if FULL_ARRAY_DEBUG
+inline void PrintMemoryInfo( void ){ for( size_t i=0 ; i
+Array< C > memcpy( Array< C > destination , const void* source , size_t size )
+{
+ if( size>destination.maximum()*sizeof(C) )
+ {
+ fprintf( stderr , "Size of copy exceeds destination maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( &destination[0] , source , size );
+ return destination;
+}
+template< class C , class D >
+Array< C > memcpy( Array< C > destination , Array< D > source , size_t size )
+{
+ if( size>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of copy exceeds destination maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( &destination[0] , &source[0] , size );
+ return destination;
+}
+template< class C , class D >
+Array< C > memcpy( Array< C > destination , ConstArray< D > source , size_t size )
+{
+ if( size>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of copy exceeds destination maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( &destination[0] , &source[0] , size );
+ return destination;
+}
+template< class D >
+void* memcpy( void* destination , Array< D > source , size_t size )
+{
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( destination , &source[0] , size );
+ return destination;
+}
+template< class D >
+void* memcpy( void* destination , ConstArray< D > source , size_t size )
+{
+ if( size>source.maximum()*sizeof( D ) )
+ {
+ fprintf( stderr , "Size of copy exceeds source maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( source.maximum()*sizeof( D ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memcpy( destination , &source[0] , size );
+ return destination;
+}
+template< class C >
+Array< C > memset( Array< C > destination , int value , size_t size )
+{
+ if( size>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of set exceeds destination maximum: %lld > %lld\n" , ( long long )( size ) , ( long long )( destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( size ) memset( &destination[0] , value , size );
+ return destination;
+}
+
+template< class C >
+size_t fread( Array< C > destination , size_t eSize , size_t count , FILE* fp )
+{
+ if( count*eSize>destination.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of read exceeds source maximum: %lld > %lld\n" , ( long long )( count*eSize ) , ( long long )( destination.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ return fread( &destination[0] , eSize , count , fp );
+}
+template< class C >
+size_t fwrite( Array< C > source , size_t eSize , size_t count , FILE* fp )
+{
+ if( count*eSize>source.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of write exceeds source maximum: %lld > %lld\n" , ( long long )( count*eSize ) , ( long long )( source.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ return fwrite( &source[0] , eSize , count , fp );
+}
+template< class C >
+size_t fwrite( ConstArray< C > source , size_t eSize , size_t count , FILE* fp )
+{
+ if( count*eSize>source.maximum()*sizeof( C ) )
+ {
+ fprintf( stderr , "Size of write exceeds source maximum: %lld > %lld\n" , ( long long )( count*eSize ) , ( long long )( source.maximum()*sizeof( C ) ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ return fwrite( &source[0] , eSize , count , fp );
+}
+template< class C >
+void qsort( Array< C > base , size_t numElements , size_t elementSize , int (*compareFunction)( const void* , const void* ) )
+{
+ if( sizeof(C)!=elementSize )
+ {
+ fprintf( stderr , "Element sizes differ: %lld != %lld\n" , ( long long )( sizeof(C) ) , ( long long )( elementSize ) );
+ ASSERT( 0 );
+ exit( 0 );
+ }
+ if( base.minimum()>0 || base.maximum()
+struct BSplineElementCoefficients
+{
+ int coeffs[Degree+1];
+ BSplineElementCoefficients( void ){ memset( coeffs , 0 , sizeof( int ) * ( Degree+1 ) ); }
+ int& operator[]( int idx ){ return coeffs[idx]; }
+ const int& operator[]( int idx ) const { return coeffs[idx]; }
+};
+template< int Degree >
+struct BSplineElements : public std::vector< BSplineElementCoefficients< Degree > >
+{
+ static const int _off = (Degree+1)/2;
+ void _addLeft ( int offset , int boundary );
+ void _addRight( int offset , int boundary );
+public:
+ enum
+ {
+ NONE = 0,
+ DIRICHLET = -1,
+ NEUMANN = 1
+ };
+ // Coefficients are ordered as "/" "-" "\"
+ int denominator;
+
+ BSplineElements( void ) { denominator = 1; }
+ BSplineElements( int res , int offset , int boundary=NONE , int inset=0 );
+
+ void upSample( BSplineElements& high ) const;
+ void differentiate( BSplineElements< Degree-1 >& d ) const;
+
+ void print( FILE* fp=stdout ) const
+ {
+ for( int i=0 ; i >::size() ; i++ )
+ {
+ printf( "%d]" , i );
+ for( int j=0 ; j<=Degree ; j++ ) printf( " %d" , (*this)[i][j] );
+ printf( " (%d)\n" , denominator );
+ }
+ }
+};
+
+template< int Degree >
+class BSplineData
+{
+ int _boundaryType;
+ double _vvIntegrals[Degree+1][Degree+1];
+ double _vdIntegrals[Degree+1][Degree ];
+ double _dvIntegrals[Degree ][Degree+1];
+ double _ddIntegrals[Degree ][Degree ];
+
+public:
+ struct Integrator
+ {
+ struct IntegralTables
+ {
+ double vv_ccIntegrals[2*Degree+1][2*Degree+1] , vv_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ double dv_ccIntegrals[2*Degree+1][2*Degree+1] , dv_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ double vd_ccIntegrals[2*Degree+1][2*Degree+1] , vd_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ double dd_ccIntegrals[2*Degree+1][2*Degree+1] , dd_cpIntegrals[(2*Degree+1)*2][2*Degree+1];
+ };
+ std::vector< IntegralTables > iTables;
+ double dot( int depth , int off1 , int off2 , bool d1 , bool d2 , bool childParent=false ) const;
+ };
+ double dot( int depth1 , int off1 , int depth2 , int off2 , bool d1 , bool d2 , bool inset=false ) const;
+ void setIntegrator( Integrator& integrator , bool inset , bool useDotRatios=false ) const;
+ template< int Radius >
+ struct CenterEvaluator
+ {
+ struct ValueTables
+ {
+ double vValues[2*Degree+1][ 3*(2*Radius+1) ];
+ double dValues[2*Degree+1][ 3*(2*Radius+1) ];
+ };
+ std::vector< ValueTables > vTables;
+ double value( int depth , int off1 , int off2 , bool d , bool childParent=false ) const;
+ };
+ template< int Radius >
+ void setCenterEvaluator( CenterEvaluator< Radius >& evaluator , double smoothingRadius , double dSmoothingRadius, bool inset ) const;
+ double value( int depth , int off , double smoothingRadius , double s , bool d , bool inset=false ) const;
+ template< int Radius >
+ struct CornerEvaluator
+ {
+ struct ValueTables
+ {
+ double vValues[2*Degree+1][4*Radius+3];
+ double dValues[2*Degree+1][4*Radius+3];
+ };
+ std::vector< ValueTables > vTables;
+ double value( int depth , int off1 , int c1 , int off2 , bool d , bool childParent=false ) const;
+ };
+ template< int Radius >
+ void setCornerEvaluator( CornerEvaluator< Radius >& evaluator , double smoothingRadius , double dSmoothingRadius, bool inset ) const;
+
+ struct BSplineComponents
+ {
+ Polynomial< Degree > polys[Degree+1];
+ Polynomial< Degree >& operator[] ( int idx ) { return polys[idx]; }
+ const Polynomial< Degree >& operator[] ( int idx ) const { return polys[idx]; }
+ void printnl( void ) const { for( int d=0 ; d<=Degree ; d++ ) polys[d].printnl(); }
+ BSplineComponents scale( double s ) const { BSplineComponents b ; for( int d=0 ; d<=Degree ; d++ ) b[d] = polys[d].scale(s) ; return b; }
+ BSplineComponents shift( double s ) const { BSplineComponents b ; for( int d=0 ; d<=Degree ; d++ ) b[d] = polys[d].shift(s) ; return b; }
+ };
+
+ int depth;
+ size_t functionCount , sampleCount;
+ PPolynomial< Degree > baseFunction , leftBaseFunction , rightBaseFunction , leftRightBaseFunction;
+ PPolynomial< Degree-1 > dBaseFunction , dLeftBaseFunction , dRightBaseFunction , dLeftRightBaseFunction;
+ BSplineComponents baseBSpline , leftBSpline , rightBSpline , leftRightBSpline;
+ Pointer( PPolynomial< Degree > ) baseFunctions;
+ Pointer( BSplineComponents ) baseBSplines;
+
+ BSplineData( void );
+
+ const static int VV_DOT_FLAG = 1;
+ const static int DV_DOT_FLAG = 2;
+ const static int DD_DOT_FLAG = 4;
+ const static int VALUE_FLAG = 1;
+ const static int D_VALUE_FLAG = 2;
+ template< class Real >
+ struct DotTables
+ {
+ size_t functionCount;
+ Pointer( Real ) vvDotTable;
+ Pointer( Real ) dvDotTable;
+ Pointer( Real ) ddDotTable;
+
+ DotTables( void );
+ ~DotTables( void );
+
+ inline size_t Index( int i1 , int i2 ) const;
+ static inline size_t SymmetricIndex( int i1 , int i2 );
+ static inline int SymmetricIndex( int i1 , int i2 , size_t& index );
+ };
+ template< class Real >
+ struct ValueTables
+ {
+ size_t functionCount , sampleCount;
+ Pointer( Real ) valueTable;
+ Pointer( Real ) dValueTable;
+
+ ValueTables( void );
+ ~ValueTables( void );
+
+ inline size_t Index( int i1 , int i2 ) const;
+ void setSampleSpan( int idx , int& start , int& end , double smooth=0 ) const;
+ };
+ void set( int maxDepth , int boundaryType=BSplineElements< Degree >::NONE );
+ template< class Real >
+ typename BSplineData< Degree >::template DotTables< Real > getDotTables( int flags , bool useDotRatios=true , bool inset=false ) const;
+ template< class Real >
+ typename BSplineData< Degree >::template ValueTables< Real > getValueTables( int flags , double valueSmooth=0 , double normalSmooth=0 ) const;
+};
+
+template< int Degree1 , int Degree2 > void SetBSplineElementIntegrals( double integrals[Degree1+1][Degree2+1] );
+
+#include "BSplineData.inl"
+#endif // BSPLINE_DATA_INCLUDED
\ No newline at end of file
diff --git a/Src/BSplineData.inl b/Src/BSplineData.inl
new file mode 100644
index 0000000..ef9cb7b
--- /dev/null
+++ b/Src/BSplineData.inl
@@ -0,0 +1,733 @@
+/*
+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.
+*/
+
+/////////////////
+// BSplineData //
+/////////////////
+// Support[i]:
+// Odd: i +/- 0.5 * ( 1 + Degree )
+// i - 0.5 * ( 1 + Degree ) < 0
+// <=> i < 0.5 * ( 1 + Degree )
+// i + 0.5 * ( 1 + Degree ) > 0
+// <=> i > - 0.5 * ( 1 + Degree )
+// i + 0.5 * ( 1 + Degree ) > r
+// <=> i > r - 0.5 * ( 1 + Degree )
+// i - 0.5 * ( 1 + Degree ) < r
+// <=> i < r + 0.5 * ( 1 + Degree )
+// Even: i + 0.5 +/- 0.5 * ( 1 + Degree )
+// i - 0.5 * Degree < 0
+// <=> i < 0.5 * Degree
+// i + 1 + 0.5 * Degree > 0
+// <=> i > -1 - 0.5 * Degree
+// i + 1 + 0.5 * Degree > r
+// <=> i > r - 1 - 0.5 * Degree
+// i - 0.5 * Degree < r
+// <=> i < r + 0.5 * Degree
+template< int Degree > inline bool LeftOverlap( unsigned int depth , int offset )
+{
+ offset <<= 1;
+ if( Degree & 1 ) return (offset < 1+Degree) && (offset > -1-Degree );
+ else return (offset < Degree) && (offset > -2-Degree );
+}
+template< int Degree > inline bool RightOverlap( unsigned int depth , int offset )
+{
+ offset <<= 1;
+ int r = 1<<(depth+1);
+ if( Degree & 1 ) return (offset > 2-1-Degree) && (offset < 2+1+Degree );
+ else return (offset > 2-2-Degree) && (offset < 2+ Degree );
+}
+template< int Degree > inline int ReflectLeft( unsigned int depth , int offset )
+{
+ if( Degree&1 ) return -offset;
+ else return -1-offset;
+}
+template< int Degree > inline int ReflectRight( unsigned int depth , int offset )
+{
+ int r = 1<<(depth+1);
+ if( Degree&1 ) return r -offset;
+ else return r-1-offset;
+}
+
+template< int Degree >
+BSplineData< Degree >::BSplineData( void )
+{
+ functionCount = sampleCount = 0;
+ SetBSplineElementIntegrals< Degree , Degree >( _vvIntegrals );
+ SetBSplineElementIntegrals< Degree , Degree-1 >( _vdIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree >( _dvIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree-1 >( _ddIntegrals );
+}
+
+template< int Degree >
+double BSplineData< Degree >::Integrator::dot( int depth , int off1 , int off2 , bool d1 , bool d2 , bool childParent ) const
+{
+ if( depth<0 || depth>=int( iTables.size() ) ) return 0.;
+ const typename Integrator::IntegralTables& iTable = iTables[depth];
+ if( childParent )
+ {
+ int c = off1&1;
+ off1 >>= 1 , depth--;
+ int ii , d = off2-off1 , res = (1<=res || off2>=res || d<-Degree || d>Degree ) return 0;
+ if ( off1< Degree ) ii = off1;
+ else if( off1>=res-Degree ) ii = 2*Degree + off1 - (res-1);
+ else ii = Degree;
+ if ( d1 && d2 ) return iTable.dd_cpIntegrals[2*ii+c][d+Degree];
+ else if( d1 ) return iTable.dv_cpIntegrals[2*ii+c][d+Degree];
+ else if( d2 ) return iTable.vd_cpIntegrals[2*ii+c][d+Degree];
+ else return iTable.vv_cpIntegrals[2*ii+c][d+Degree];
+ }
+ else
+ {
+ int ii , d = off2-off1 , res = (1<=res || off2>=res || d<-Degree || d>Degree ) return 0;
+ if ( off1< Degree ) ii = off1;
+ else if( off1>=res-Degree ) ii = 2*Degree + off1 - (res-1);
+ else ii = Degree;
+ if ( d1 && d2 ) return iTable.dd_ccIntegrals[ii][d+Degree];
+ else if( d1 ) return iTable.dv_ccIntegrals[ii][d+Degree];
+ else if( d2 ) return iTable.vd_ccIntegrals[ii][d+Degree];
+ else return iTable.vv_ccIntegrals[ii][d+Degree];
+ }
+}
+template< int Degree >
+template< int Radius >
+double BSplineData< Degree >::CenterEvaluator< Radius >::value( int depth , int off1 , int off2 , bool d , bool childParent ) const
+{
+ if( depth<0 || depth>=int( vTables.size() ) ) return 0.;
+ if( childParent )
+ {
+ int c = off1&1;
+ off1 >>= 1 , depth--;
+ const typename CenterEvaluator::ValueTables& vTable = vTables[depth];
+ int ii , dd = off1-off2 , res = (1<=res || off2>=res || dd<-Radius || dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*3+2*c];
+ else return vTable.vValues[ii][(dd+Radius)*3+2*c];
+ }
+ else
+ {
+ const typename CenterEvaluator::ValueTables& vTable = vTables[depth];
+ int ii , dd = off1-off2 , res = (1<=res || off2>=res || dd<-Radius || dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*3+1];
+ else return vTable.vValues[ii][(dd+Radius)*3+1];
+ }
+}
+template< int Degree >
+template< int Radius >
+double BSplineData< Degree >::CornerEvaluator< Radius >::value( int depth , int off1 , int c1 , int off2 , bool d , bool childParent ) const
+{
+ if( c1<0 || c1>=2 )
+ {
+ fprintf( stderr , "[WARNING] Clamping corner to {0,1}\n" );
+ c1 = std::max< int >( 0 , std::min< int >( c1 , 1 ) );
+ }
+ if( depth<0 || depth>=int( vTables.size() ) ) return 0.;
+ if( childParent )
+ {
+ int c = off1&1;
+ off1 >>= 1 , depth--;
+ const typename CornerEvaluator::ValueTables& vTable = vTables[depth];
+ int ii , dd = off1-off2 , res = (1<=res || off2>=res || dd<-Radius || dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*2+c+c1];
+ else return vTable.vValues[ii][(dd+Radius)*2+c+c1];
+ }
+ else
+ {
+ const typename CornerEvaluator::ValueTables& vTable = vTables[depth];
+ int ii , dd = off1-off2 , res = (1<=res || off2>=res || dd<-Radius || dd>Radius ) return 0;
+ if ( off2< Degree ) ii = off2;
+ else if( off2>=res-Degree ) ii = 2*Degree + off2 - (res-1);
+ else ii = Degree;
+ if( d ) return vTable.dValues[ii][(dd+Radius)*2+2*c1];
+ else return vTable.vValues[ii][(dd+Radius)*2+2*c1];
+ }
+}
+template< int Degree >
+void BSplineData< Degree >::set( int maxDepth , int boundaryType )
+{
+ _boundaryType = boundaryType;
+
+ depth = maxDepth;
+ // [Warning] This assumes that the functions spacing is dual
+ functionCount = BinaryNode::CumulativeCenterCount( depth );
+ sampleCount = BinaryNode::CenterCount( depth ) + BinaryNode::CornerCount( depth );
+ baseFunctions = NewPointer< PPolynomial< Degree > >( functionCount );
+ baseBSplines = NewPointer< BSplineComponents >( functionCount );
+
+ baseFunction = PPolynomial< Degree >::BSpline();
+ for( int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree >::BSplineComponent( i ).shift( double(-(Degree+1)/2) + i - 0.5 );
+ dBaseFunction = baseFunction.derivative();
+ StartingPolynomial< Degree > sPolys[Degree+4];
+
+ for( int i=0 ; i=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1];
+ for( int j=0 ; j=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1].shift( 1 ) * _boundaryType;
+ for( int j=0 ; j=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1]; // The centered B-Spline
+ if( i>=2 && i<=Degree+2 ) sPolys[i].p += baseBSpline[i-2].shift( 1 ) * _boundaryType; // The right-shifted B-spline
+ for( int j=0 ; j
+double BSplineData< Degree >::dot( int depth1 , int off1 , int depth2 , int off2 , bool d1 , bool d2 , bool inset ) const
+{
+ const int _Degree1 = (d1 ? (Degree-1) : Degree) , _Degree2 = (d2 ? (Degree-1) : Degree);
+ int sums[ Degree+1 ][ Degree+1 ];
+
+ int depth = std::max< int >( depth1 , depth2 );
+
+ BSplineElements< Degree > b1( 1< b;
+ while( depth1 db1 , db2;
+ b1.differentiate( db1 ) , b2.differentiate( db2 );
+
+ int start1=-1 , end1=-1 , start2=-1 , end2=-1;
+ for( int i=0 ; i=end2 || start2>=end1 ) return 0.;
+ int start = std::max< int >( start1 , start2 ) , end = std::min< int >( end1 , end2 );
+ memset( sums , 0 , sizeof( sums ) );
+ for( int i=start ; i
+double BSplineData< Degree >::value( int depth , int off , double smoothingRadius , double s , bool d , bool inset ) const
+{
+ PPolynomial< Degree+1 > function;
+ PPolynomial< Degree > dFunction;
+
+ if( off<0 || off>=(1<0 ) function = baseFunctions[idx].MovingAverage( smoothingRadius );
+ else function = baseFunctions[idx];
+ dFunction = function.derivative();
+
+ if( d ) return dFunction(s);
+ else return function(s);
+}
+template< int Degree >
+void BSplineData< Degree >::setIntegrator( Integrator& integrator , bool inset , bool useDotRatios ) const
+{
+ integrator.iTables.resize( depth+1 );
+ for( int d=0 ; d<=depth ; d++ ) for( int i=0 ; i<=2*Degree ; i++ ) for( int j=-Degree ; j<=Degree ; j++ )
+ {
+ int res = 1<
+template< int Radius >
+void BSplineData< Degree >::setCenterEvaluator( CenterEvaluator< Radius >& evaluator , double smoothingRadius , double dSmoothingRadius , bool inset ) const
+{
+ evaluator.vTables.resize( depth+1 );
+ for( int d=0 ; d<=depth ; d++ ) for( int i=0 ; i<=2*Degree ; i++ ) for( int j=-Radius ; j<=Radius ; j++ ) for( int k=-1 ; k<=1 ; k++ )
+ {
+ int res = 1<
+template< int Radius >
+void BSplineData< Degree >::setCornerEvaluator( CornerEvaluator< Radius >& evaluator , double smoothingRadius , double dSmoothingRadius , bool inset ) const
+{
+ evaluator.vTables.resize( depth+1 );
+ for( int d=0 ; d<=depth ; d++ ) for( int i=0 ; i<=2*Degree ; i++ ) for( int j=-Radius ; j<=Radius ; j++ ) for( int k=0 ; k<=2 ; k++ )
+ {
+ int res = 1<
+template< class Real >
+BSplineData< Degree >::DotTables< Real >::DotTables( void )
+{
+ vvDotTable = NullPointer< Real >();
+ dvDotTable = NullPointer< Real >();
+ ddDotTable = NullPointer< Real >();
+}
+template< int Degree >
+template< class Real >
+BSplineData< Degree >::DotTables< Real >::~DotTables( void )
+{
+ DeletePointer( vvDotTable );
+ DeletePointer( dvDotTable );
+ DeletePointer( ddDotTable );
+}
+template< int Degree >
+template< class Real >
+inline size_t BSplineData< Degree >::DotTables< Real >::Index( int i1 , int i2 ) const { return size_t(i1)*functionCount + size_t(i2); }
+template< int Degree >
+template< class Real >
+inline size_t BSplineData< Degree >::DotTables< Real >::SymmetricIndex( int i1 , int i2 )
+{
+ size_t _i1 = i1 , _i2 = i2;
+ if( i1>i2 ) return ((_i1*_i1+i1)>>1)+_i2;
+ else return ((_i2*_i2+i2)>>1)+_i1;
+}
+template< int Degree >
+template< class Real >
+inline int BSplineData< Degree >::DotTables< Real >::SymmetricIndex( int i1 , int i2 , size_t& index )
+{
+ size_t _i1 = i1 , _i2 = i2;
+ if( i1>1)+_i1;
+ return 1;
+ }
+ else
+ {
+ index = ((_i1*_i1+_i1)>>1)+_i2;
+ return 0;
+ }
+}
+template< int Degree >
+template< class Real >
+typename BSplineData< Degree >::template DotTables< Real > BSplineData< Degree >::getDotTables( int flags , bool useDotRatios , bool inset ) const
+{
+ typename BSplineData< Degree >::template DotTables< Real > dTables;
+ dTables.functionCount = functionCount;
+
+ size_t size = ( functionCount*functionCount + functionCount )>>1;
+ size_t fullSize = functionCount*functionCount;
+ if( flags & VV_DOT_FLAG )
+ {
+ dTables.vvDotTable = NewPointer< Real >( size );
+ memset( dTables.vvDotTable , 0 , sizeof(Real)*size );
+ }
+ if( flags & DV_DOT_FLAG )
+ {
+ dTables.dvDotTable = NewPointer< Real >( fullSize );
+ memset( dTables.dvDotTable , 0 , sizeof(Real)*fullSize );
+ }
+ if( flags & DD_DOT_FLAG )
+ {
+ dTables.ddDotTable = NewPointer< Real >( size );
+ memset( dTables.ddDotTable , 0 , sizeof(Real)*size );
+ }
+ int vvSums[Degree+1][Degree+1];
+ int vdSums[Degree+1][Degree ];
+ int dvSums[Degree ][Degree+1];
+ int ddSums[Degree ][Degree ];
+ double vvIntegrals[Degree+1][Degree+1];
+ double vdIntegrals[Degree+1][Degree ];
+ double dvIntegrals[Degree ][Degree+1];
+ double ddIntegrals[Degree ][Degree ];
+ SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
+ SetBSplineElementIntegrals< Degree , Degree-1 >( vdIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree >( dvIntegrals );
+ SetBSplineElementIntegrals< Degree-1 , Degree-1 >( ddIntegrals );
+
+ for( int d1=0 ; d1<=depth ; d1++ ) for( int off1=0 ; off1<(1< b1( 1< db1;
+ b1.differentiate( db1 );
+ int start1 , end1;
+
+ start1 = -1 , end1 = -1;
+ for( int i=0 ; i=end1 || start1>=end2 ) continue;
+ start2 = std::max< int >( start1 , start2 );
+ end2 = std::min< int >( end1 , end2 );
+ if( d1==d2 && off2 b2( 1< db2;
+ b2.differentiate( db2 );
+
+ size_t idx = DotTables< Real >::SymmetricIndex( ii , jj );
+ size_t idx1 = DotTables< Real >::Index( ii , jj ) , idx2 = DotTables< Real >::Index( jj , ii );
+
+ memset( vvSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree+1 ) );
+ memset( vdSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree ) );
+ memset( dvSums , 0 , sizeof( int ) * ( Degree ) * ( Degree+1 ) );
+ memset( ddSums , 0 , sizeof( int ) * ( Degree ) * ( Degree ) );
+ for( int i=start2 ; i b;
+ b = b1;
+ b.upSample( b1 );
+ b1.differentiate( db1 );
+ start1 = -1;
+ for( int i=0 ; i
+template< class Real >
+BSplineData< Degree >::ValueTables< Real >::ValueTables( void )
+{
+ valueTable = NullPointer< Real >();
+ dValueTable = NullPointer< Real >();
+}
+template< int Degree >
+template< class Real >
+BSplineData< Degree >::ValueTables< Real >::~ValueTables( void )
+{
+ DeletePointer( valueTable );
+ DeletePointer( dValueTable );
+}
+template< int Degree >
+template< class Real >
+inline size_t BSplineData< Degree >::ValueTables< Real >::Index( int i1 , int i2 ) const { return size_t(i1)*functionCount + size_t(i2); }
+template< int Degree >
+template< class Real >
+typename BSplineData< Degree >::template ValueTables< Real > BSplineData< Degree >::getValueTables( int flags , double valueSmooth , double derivativeSmooth ) const
+{
+ typename BSplineData< Degree >::template ValueTables< Real > vTables;
+ vTables.functionCount = functionCount;
+ vTables.sampleCount = sampleCount;
+
+ if( flags & VALUE_FLAG ) vTables.valueTable = NewPointer< Real >( functionCount*sampleCount );
+ if( flags & D_VALUE_FLAG ) vTables.dValueTable = NewPointer< Real >( functionCount*sampleCount );
+ PPolynomial< Degree+1 > function;
+ PPolynomial< Degree > dFunction;
+ for( size_t i=0 ; i0 ) function=baseFunctions[i].MovingAverage( valueSmooth );
+ else function=baseFunctions[i];
+ if( derivativeSmooth>0 ) dFunction=baseFunctions[i].derivative().MovingAverage( derivativeSmooth );
+ else dFunction=baseFunctions[i].derivative();
+
+ for( size_t j=0 ; j
+template< class Real >
+void BSplineData< Degree >::ValueTables< Real >::setSampleSpan( int idx , int& start , int& end , double smooth ) const
+{
+ int d , off , res;
+ BinaryNode::DepthAndOffset( idx , d , off );
+ res = 1<_start && (start-1)/(sampleCount-1)<=_start
+ // => start > _start * (sampleCount-1 ) && start <= _start*(sampleCount-1) + 1
+ // => _start * (sampleCount-1) + 1 >= start > _start * (sampleCount-1)
+ start = int( floor( _start * (sampleCount-1) + 1 ) );
+ if( start<0 ) start = 0;
+ // (end)/(sampleCount-1)<_end && (end+1)/(sampleCount-1)>=_end
+ // => end < _end * (sampleCount-1 ) && end >= _end*(sampleCount-1) - 1
+ // => _end * (sampleCount-1) > end >= _end * (sampleCount-1) - 1
+ end = int( ceil( _end * (sampleCount-1) - 1 ) );
+ if( end>=int(sampleCount) ) end = int(sampleCount)-1;
+}
+
+
+/////////////////////
+// BSplineElements //
+/////////////////////
+template< int Degree >
+BSplineElements< Degree >::BSplineElements( int res , int offset , int boundary , int inset )
+{
+ denominator = 1;
+ std::vector< BSplineElementCoefficients< Degree > >::resize( res , BSplineElementCoefficients< Degree >() );
+
+ for( int i=0 ; i<=Degree ; i++ )
+ {
+ int idx = -_off + offset + i;
+ if( idx>=0 && idx
+void BSplineElements< Degree >::_addLeft( int offset , int boundary )
+{
+ int res = int( std::vector< BSplineElementCoefficients< Degree > >::size() );
+ bool set = false;
+ for( int i=0 ; i<=Degree ; i++ )
+ {
+ int idx = -_off + offset + i;
+ if( idx>=0 && idx
+void BSplineElements< Degree >::_addRight( int offset , int boundary )
+{
+ int res = int( std::vector< BSplineElementCoefficients< Degree > >::size() );
+ bool set = false;
+ for( int i=0 ; i<=Degree ; i++ )
+ {
+ int idx = -_off + offset + i;
+ if( idx>=0 && idx
+void BSplineElements< Degree >::upSample( BSplineElements< Degree >& high ) const
+{
+ fprintf( stderr , "[ERROR] B-spline up-sampling not supported for degree %d\n" , Degree );
+ exit( 0 );
+}
+template<>
+void BSplineElements< 1 >::upSample( BSplineElements< 1 >& high ) const
+{
+ high.resize( size()*2 );
+ high.assign( high.size() , BSplineElementCoefficients<1>() );
+ for( int i=0 ; i
+void BSplineElements< 2 >::upSample( BSplineElements< 2 >& high ) const
+{
+ // /----\
+ // / \
+ // / \ = 1 /--\ +3 /--\ +3 /--\ +1 /--\
+ // / \ / \ / \ / \ / \
+ // |----------| |----------| |----------| |----------| |----------|
+
+ high.resize( size()*2 );
+ high.assign( high.size() , BSplineElementCoefficients<2>() );
+ for( int i=0 ; i
+void BSplineElements< Degree >::differentiate( BSplineElements< Degree-1 >& d ) const
+{
+ d.resize( std::vector< BSplineElementCoefficients< Degree > >::size() );
+ d.assign( d.size() , BSplineElementCoefficients< Degree-1 >() );
+ for( int i=0 ; i >::size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
+ {
+ if( j-1>=0 ) d[i][j-1] -= (*this)[i][j];
+ if( j
+void SetBSplineElementIntegrals( double integrals[Degree1+1][Degree2+1] )
+{
+ for( int i=0 ; i<=Degree1 ; i++ )
+ {
+ Polynomial< Degree1 > p1 = Polynomial< Degree1 >::BSplineComponent( i );
+ for( int j=0 ; j<=Degree2 ; j++ )
+ {
+ Polynomial< Degree2 > p2 = Polynomial< Degree2 >::BSplineComponent( j );
+ integrals[i][j] = ( p1 * p2 ).integral( 0 , 1 );
+ }
+ }
+}
diff --git a/Src/BinaryNode.h b/Src/BinaryNode.h
new file mode 100644
index 0000000..96f3061
--- /dev/null
+++ b/Src/BinaryNode.h
@@ -0,0 +1,78 @@
+/*
+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 BINARY_NODE_INCLUDED
+#define BINARY_NODE_INCLUDED
+
+#define MSVC_2010_FIX 1
+
+
+class BinaryNode
+{
+public:
+ static inline int CenterCount( int depth ) { return 1< static inline Real CornerIndexPosition(int index,int maxDepth){ return Real(index)/(1< static inline Real Width(int depth){ return Real(1.0/(1< static inline void CenterAndWidth( int depth , int offset , Real& center , Real& width )
+ {
+ width=Real (1.0/(1< static inline void CenterAndWidth( int idx , Real& center , Real& width )
+ {
+ int depth , offset;
+ DepthAndOffset( idx , depth , offset );
+ CenterAndWidth( depth , offset , center , width );
+ }
+ static inline void DepthAndOffset( int idx , int& depth , int& offset )
+ {
+ int i=idx+1;
+#if MSVC_2010_FIX
+ depth = 0;
+#else // !MSVC_2010_FIX
+ depth = -1;
+#endif // MSVC_2010_FIX
+ while( i )
+ {
+ i >>= 1;
+ depth++;
+ }
+#if MSVC_2010_FIX
+ depth--;
+#endif // MSVC_2010_FIX
+ offset = ( idx+1 ) - (1<
+#include
+#include
+#include
+#include "CmdLineParser.h"
+
+
+#ifdef WIN32
+int strcasecmp(char* c1,char* c2){return _stricmp(c1,c2);}
+#endif
+
+cmdLineReadable::cmdLineReadable(const char* name)
+{
+ set=false;
+ this->name=new char[strlen(name)+1];
+ strcpy(this->name,name);
+}
+cmdLineReadable::~cmdLineReadable(void)
+{
+ if(name) delete[] name;
+ name=NULL;
+}
+int cmdLineReadable::read(char**,int){
+ set=true;
+ return 0;
+}
+void cmdLineReadable::writeValue(char* str)
+{
+ str[0] = 0;
+}
+
+////////////////
+// cmdLineInt //
+////////////////
+cmdLineInt::cmdLineInt(const char* name) : cmdLineReadable(name) {value=0;}
+cmdLineInt::cmdLineInt(const char* name,const int& v) : cmdLineReadable(name) {value=v;}
+int cmdLineInt::read(char** argv,int argc){
+ if(argc>0){
+ value=atoi(argv[0]);
+ set=true;
+ return 1;
+ }
+ else{return 0;}
+}
+void cmdLineInt::writeValue(char* str)
+{
+ sprintf(str,"%d",value);
+}
+
+//////////////////
+// cmdLineFloat //
+//////////////////
+cmdLineFloat::cmdLineFloat(const char* name) : cmdLineReadable(name) {value=0;}
+cmdLineFloat::cmdLineFloat(const char* name, const float& v) : cmdLineReadable(name) {value=v;}
+int cmdLineFloat::read(char** argv,int argc){
+ if(argc>0){
+ value=(float)atof(argv[0]);
+ set=true;
+ return 1;
+ }
+ else{return 0;}
+}
+void cmdLineFloat::writeValue(char* str)
+{
+ sprintf(str,"%f",value);
+}
+
+///////////////////
+// cmdLineString //
+///////////////////
+cmdLineString::cmdLineString(const char* name) : cmdLineReadable(name) {value=NULL;}
+cmdLineString::~cmdLineString(void)
+{
+ if(value) delete[] value;
+ value=NULL;
+}
+int cmdLineString::read(char** argv,int argc){
+ if(argc>0)
+ {
+ value=new char[strlen(argv[0])+1];
+ strcpy(value,argv[0]);
+ set=true;
+ return 1;
+ }
+ else{return 0;}
+}
+void cmdLineString::writeValue(char* str)
+{
+ sprintf(str,"%s",value);
+}
+
+////////////////////
+// cmdLineStrings //
+////////////////////
+cmdLineStrings::cmdLineStrings(const char* name,int Dim) : cmdLineReadable(name)
+{
+ this->Dim=Dim;
+ values=new char*[Dim];
+ for(int i=0;i=Dim)
+ {
+ for(int i=0;i 0)
+ {
+ if (argv[0][0] == '-' && argv[0][1]=='-')
+ {
+ for(i=0;iname))
+ {
+ argv++, argc--;
+ j=readable[i]->read(argv,argc);
+ argv+=j,argc-=j;
+ break;
+ }
+ }
+ if(i==num){
+ if(dumpError)
+ {
+ fprintf(stderr, "invalid option: %s\n",*argv);
+ fprintf(stderr, "possible options are:\n");
+ for(i=0;iname);
+ }
+ argv++, argc--;
+ }
+ }
+ else
+ {
+ if(dumpError)
+ {
+ fprintf(stderr, "invalid option: %s\n", *argv);
+ fprintf(stderr, " options must start with a \'--\'\n");
+ }
+ argv++, argc--;
+ }
+ }
+}
+char** ReadWords(const char* fileName,int& cnt)
+{
+ char** names;
+ char temp[500];
+ FILE* fp;
+
+ fp=fopen(fileName,"r");
+ if(!fp){return NULL;}
+ cnt=0;
+ while(fscanf(fp," %s ",temp)==1){cnt++;}
+ fclose(fp);
+
+ names=new char*[cnt];
+ if(!names){return NULL;}
+
+ fp=fopen(fileName,"r");
+ if(!fp){
+ delete[] names;
+ cnt=0;
+ return NULL;
+ }
+ cnt=0;
+ while(fscanf(fp," %s ",temp)==1){
+ names[cnt]=new char[strlen(temp)+1];
+ if(!names){
+ for(int j=0;j
+#include
+
+
+#ifdef WIN32
+int strcasecmp(char* c1,char* c2);
+#endif
+
+class cmdLineReadable{
+public:
+ bool set;
+ char* name;
+ cmdLineReadable(const char* name);
+ virtual ~cmdLineReadable(void);
+ virtual int read(char** argv,int argc);
+ virtual void writeValue(char* str);
+};
+
+class cmdLineInt : public cmdLineReadable {
+public:
+ int value;
+ cmdLineInt(const char* name);
+ cmdLineInt(const char* name,const int& v);
+ int read(char** argv,int argc);
+ void writeValue(char* str);
+};
+template
+class cmdLineIntArray : public cmdLineReadable {
+public:
+ int values[Dim];
+ cmdLineIntArray(const char* name);
+ cmdLineIntArray(const char* name,const int v[Dim]);
+ int read(char** argv,int argc);
+ void writeValue(char* str);
+};
+
+class cmdLineFloat : public cmdLineReadable {
+public:
+ float value;
+ cmdLineFloat(const char* name);
+ cmdLineFloat(const char* name,const float& f);
+ int read(char** argv,int argc);
+ void writeValue(char* str);
+};
+template
+class cmdLineFloatArray : public cmdLineReadable {
+public:
+ float values[Dim];
+ cmdLineFloatArray(const char* name);
+ cmdLineFloatArray(const char* name,const float f[Dim]);
+ int read(char** argv,int argc);
+ void writeValue(char* str);
+};
+class cmdLineString : public cmdLineReadable {
+public:
+ char* value;
+ cmdLineString(const char* name);
+ ~cmdLineString();
+ int read(char** argv,int argc);
+ void writeValue(char* str);
+};
+class cmdLineStrings : public cmdLineReadable {
+ int Dim;
+public:
+ char** values;
+ cmdLineStrings(const char* name,int Dim);
+ ~cmdLineStrings(void);
+ int read(char** argv,int argc);
+ void writeValue(char* str);
+};
+template
+class cmdLineStringArray : public cmdLineReadable {
+public:
+ char* values[Dim];
+ cmdLineStringArray(const char* name);
+ ~cmdLineStringArray();
+ int read(char** argv,int argc);
+ void writeValue(char* str);
+};
+
+// This reads the arguments in argc, matches them against "names" and sets
+// the values of "r" appropriately. Parameters start with "--"
+void cmdLineParse(int argc, char **argv,int num,cmdLineReadable** r,int dumpError=1);
+
+char* GetFileExtension(char* fileName);
+char* GetLocalFileName(char* fileName);
+char** ReadWords(const char* fileName,int& cnt);
+
+#include "CmdLineParser.inl"
+#endif // CMD_LINE_PARSER_INCLUDED
diff --git a/Src/CmdLineParser.inl b/Src/CmdLineParser.inl
new file mode 100644
index 0000000..eeded68
--- /dev/null
+++ b/Src/CmdLineParser.inl
@@ -0,0 +1,141 @@
+/* -*- C++ -*-
+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.
+*/
+
+/////////////////////
+// cmdLineIntArray //
+/////////////////////
+template
+cmdLineIntArray::cmdLineIntArray(const char* name) : cmdLineReadable(name)
+{
+ for(int i=0;i
+cmdLineIntArray::cmdLineIntArray(const char* name,const int v[Dim]) : cmdLineReadable(name)
+{
+ for(int i=0;i
+int cmdLineIntArray::read(char** argv,int argc)
+{
+ if(argc>=Dim)
+ {
+ for(int i=0;i
+void cmdLineIntArray::writeValue(char* str)
+{
+ char* temp=str;
+ for(int i=0;i
+cmdLineFloatArray::cmdLineFloatArray(const char* name) : cmdLineReadable(name)
+{
+ for(int i=0;i
+cmdLineFloatArray::cmdLineFloatArray(const char* name,const float f[Dim]) : cmdLineReadable(name)
+{
+ for(int i=0;i
+int cmdLineFloatArray::read(char** argv,int argc)
+{
+ if(argc>=Dim)
+ {
+ for(int i=0;i
+void cmdLineFloatArray::writeValue(char* str)
+{
+ char* temp=str;
+ for(int i=0;i
+cmdLineStringArray::cmdLineStringArray(const char* name) : cmdLineReadable(name)
+{
+ for(int i=0;i
+cmdLineStringArray::~cmdLineStringArray(void)
+{
+ for(int i=0;i
+int cmdLineStringArray::read(char** argv,int argc)
+{
+ if(argc>=Dim)
+ {
+ for(int i=0;i
+void cmdLineStringArray::writeValue(char* str)
+{
+ char* temp=str;
+ for(int i=0;i
+#include "Factor.h"
+int Factor(double a1,double a0,double roots[1][2],double EPS){
+ if(fabs(a1)<=EPS){return 0;}
+ roots[0][0]=-a0/a1;
+ roots[0][1]=0;
+ return 1;
+}
+int Factor(double a2,double a1,double a0,double roots[2][2],double EPS){
+ double d;
+ if(fabs(a2)<=EPS){return Factor(a1,a0,roots,EPS);}
+
+ d=a1*a1-4*a0*a2;
+ a1/=(2*a2);
+ if(d<0){
+ d=sqrt(-d)/(2*a2);
+ roots[0][0]=roots[1][0]=-a1;
+ roots[0][1]=-d;
+ roots[1][1]= d;
+ }
+ else{
+ d=sqrt(d)/(2*a2);
+ roots[0][1]=roots[1][1]=0;
+ roots[0][0]=-a1-d;
+ roots[1][0]=-a1+d;
+ }
+ return 2;
+}
+// Solution taken from: http://mathworld.wolfram.com/CubicFormula.html
+// and http://www.csit.fsu.edu/~burkardt/f_src/subpak/subpak.f90
+int Factor(double a3,double a2,double a1,double a0,double roots[3][2],double EPS){
+ double q,r,r2,q3;
+
+ if(fabs(a3)<=EPS){return Factor(a2,a1,a0,roots,EPS);}
+ a2/=a3;
+ a1/=a3;
+ a0/=a3;
+
+ q=-(3*a1-a2*a2)/9;
+ r=-(9*a2*a1-27*a0-2*a2*a2*a2)/54;
+ r2=r*r;
+ q3=q*q*q;
+
+ if(r20){return PI/2.0;}
+ else{return -PI/2.0;}
+ }
+ if(x>=0){return atan(y/x);}
+ else{
+ if(y>=0){return atan(y/x)+PI;}
+ else{return atan(y/x)-PI;}
+ }
+}
+double Angle(const double in[2]){
+ if((in[0]*in[0]+in[1]*in[1])==0.0){return 0;}
+ else{return ArcTan2(in[1],in[0]);}
+}
+void Sqrt(const double in[2],double out[2]){
+ double r=sqrt(sqrt(in[0]*in[0]+in[1]*in[1]));
+ double a=Angle(in)*0.5;
+ out[0]=r*cos(a);
+ out[1]=r*sin(a);
+}
+void Add(const double in1[2],const double in2[2],double out[2]){
+ out[0]=in1[0]+in2[0];
+ out[1]=in1[1]+in2[1];
+}
+void Subtract(const double in1[2],const double in2[2],double out[2]){
+ out[0]=in1[0]-in2[0];
+ out[1]=in1[1]-in2[1];
+}
+void Multiply(const double in1[2],const double in2[2],double out[2]){
+ out[0]=in1[0]*in2[0]-in1[1]*in2[1];
+ out[1]=in1[0]*in2[1]+in1[1]*in2[0];
+}
+void Divide(const double in1[2],const double in2[2],double out[2]){
+ double temp[2];
+ double l=in2[0]*in2[0]+in2[1]*in2[1];
+ temp[0]= in2[0]/l;
+ temp[1]=-in2[1]/l;
+ Multiply(in1,temp,out);
+}
+// Solution taken from: http://mathworld.wolfram.com/QuarticEquation.html
+// and http://www.csit.fsu.edu/~burkardt/f_src/subpak/subpak.f90
+int Factor(double a4,double a3,double a2,double a1,double a0,double roots[4][2],double EPS){
+ double R[2],D[2],E[2],R2[2];
+
+ if(fabs(a4)10e-8){
+ double temp1[2],temp2[2];
+ double p1[2],p2[2];
+
+ p1[0]=a3*a3*0.75-2.0*a2-R2[0];
+ p1[1]=0;
+
+ temp2[0]=((4.0*a3*a2-8.0*a1-a3*a3*a3)/4.0);
+ temp2[1]=0;
+ Divide(temp2,R,p2);
+
+ Add (p1,p2,temp1);
+ Subtract(p1,p2,temp2);
+
+ Sqrt(temp1,D);
+ Sqrt(temp2,E);
+ }
+ else{
+ R[0]=R[1]=0;
+ double temp1[2],temp2[2];
+ temp1[0]=roots[0][0]*roots[0][0]-4.0*a0;
+ temp1[1]=0;
+ Sqrt(temp1,temp2);
+ temp1[0]=a3*a3*0.75-2.0*a2+2.0*temp2[0];
+ temp1[1]= 2.0*temp2[1];
+ Sqrt(temp1,D);
+ temp1[0]=a3*a3*0.75-2.0*a2-2.0*temp2[0];
+ temp1[1]= -2.0*temp2[1];
+ Sqrt(temp1,E);
+ }
+
+ roots[0][0]=-a3/4.0+R[0]/2.0+D[0]/2.0;
+ roots[0][1]= R[1]/2.0+D[1]/2.0;
+
+ roots[1][0]=-a3/4.0+R[0]/2.0-D[0]/2.0;
+ roots[1][1]= R[1]/2.0-D[1]/2.0;
+
+ roots[2][0]=-a3/4.0-R[0]/2.0+E[0]/2.0;
+ roots[2][1]= -R[1]/2.0+E[1]/2.0;
+
+ roots[3][0]=-a3/4.0-R[0]/2.0-E[0]/2.0;
+ roots[3][1]= -R[1]/2.0-E[1]/2.0;
+ return 4;
+}
+
+int Solve(const double* eqns,const double* values,double* solutions,int dim){
+ int i,j,eIndex;
+ double v,m;
+ int *index=new int[dim];
+ int *set=new int[dim];
+ double* myEqns=new double[dim*dim];
+ double* myValues=new double[dim];
+
+ for(i=0;im){
+ m=fabs(myEqns[j*dim+i]);
+ eIndex=j;
+ }
+ }
+ if(eIndex==-1){
+ delete[] index;
+ delete[] myValues;
+ delete[] myEqns;
+ delete[] set;
+ return 0;
+ }
+ // The position in which the solution for the i-th variable can be found
+ index[i]=eIndex;
+ set[eIndex]=1;
+
+ // Normalize the equation
+ v=myEqns[eIndex*dim+i];
+ for(j=0;j
+class FunctionData{
+ bool useDotRatios;
+ int normalize;
+#if BOUNDARY_CONDITIONS
+ bool reflectBoundary;
+#endif // BOUNDARY_CONDITIONS
+public:
+ const static int DOT_FLAG = 1;
+ const static int D_DOT_FLAG = 2;
+ const static int D2_DOT_FLAG = 4;
+ const static int VALUE_FLAG = 1;
+ const static int D_VALUE_FLAG = 2;
+
+ int depth , res , res2;
+ Real *dotTable , *dDotTable , *d2DotTable;
+ Real *valueTables , *dValueTables;
+#if BOUNDARY_CONDITIONS
+ PPolynomial baseFunction , leftBaseFunction , rightBaseFunction;
+ PPolynomial dBaseFunction , dLeftBaseFunction , dRightBaseFunction;
+#else // !BOUNDARY_CONDITIONS
+ PPolynomial baseFunction;
+ PPolynomial dBaseFunction;
+#endif // BOUNDARY_CONDITIONS
+ PPolynomial* baseFunctions;
+
+ FunctionData(void);
+ ~FunctionData(void);
+
+ virtual void setDotTables(const int& flags);
+ virtual void clearDotTables(const int& flags);
+
+ virtual void setValueTables(const int& flags,const double& smooth=0);
+ virtual void setValueTables(const int& flags,const double& valueSmooth,const double& normalSmooth);
+ virtual void clearValueTables(void);
+
+ /********************************************************
+ * Sets the translates and scales of the basis function
+ * up to the prescribed depth
+ * the maximum depth
+ * the basis function
+ * how the functions should be scaled
+ * 0] Value at zero equals 1
+ * 1] Integral equals 1
+ * 2] L2-norm equals 1
+ * specifies if dot-products of derivatives
+ * should be pre-divided by function integrals
+ * spcifies if function space should be
+ * forced to be reflectively symmetric across the boundary
+ ********************************************************/
+#if BOUNDARY_CONDITIONS
+ void set( const int& maxDepth , const PPolynomial& F , const int& normalize , bool useDotRatios=true , bool reflectBoundary=false );
+#else // !BOUNDARY_CONDITIONS
+ void set(const int& maxDepth,const PPolynomial& F,const int& normalize , bool useDotRatios=true );
+#endif // BOUNDARY_CONDITIONS
+
+#if BOUNDARY_CONDITIONS
+ Real dotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const;
+ Real dDotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const;
+ Real d2DotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const;
+#else // !BOUNDARY_CONDITIONS
+ Real dotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 ) const;
+ Real dDotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 ) const;
+ Real d2DotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 ) const;
+#endif // BOUNDARY_CONDITIONS
+
+ static inline int SymmetricIndex( const int& i1 , const int& i2 );
+ static inline int SymmetricIndex( const int& i1 , const int& i2 , int& index );
+};
+
+
+#include "FunctionData.inl"
+#endif // FUNCTION_DATA_INCLUDED
\ No newline at end of file
diff --git a/Src/FunctionData.inl b/Src/FunctionData.inl
new file mode 100644
index 0000000..4e61b96
--- /dev/null
+++ b/Src/FunctionData.inl
@@ -0,0 +1,415 @@
+/*
+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.
+*/
+
+//////////////////
+// FunctionData //
+//////////////////
+template
+FunctionData::FunctionData(void)
+{
+ dotTable=dDotTable=d2DotTable=NULL;
+ valueTables=dValueTables=NULL;
+ res=0;
+}
+
+template
+FunctionData::~FunctionData(void)
+{
+ if(res)
+ {
+ if( dotTable) delete[] dotTable;
+ if( dDotTable) delete[] dDotTable;
+ if(d2DotTable) delete[] d2DotTable;
+ if( valueTables) delete[] valueTables;
+ if(dValueTables) delete[] dValueTables;
+ }
+ dotTable=dDotTable=d2DotTable=NULL;
+ valueTables=dValueTables=NULL;
+ res=0;
+}
+
+template
+#if BOUNDARY_CONDITIONS
+void FunctionData::set( const int& maxDepth , const PPolynomial& F , const int& normalize , bool useDotRatios , bool reflectBoundary )
+#else // !BOUNDARY_CONDITIONS
+void FunctionData::set(const int& maxDepth,const PPolynomial& F,const int& normalize , bool useDotRatios )
+#endif // BOUNDARY_CONDITIONS
+{
+ this->normalize = normalize;
+ this->useDotRatios = useDotRatios;
+#if BOUNDARY_CONDITIONS
+ this->reflectBoundary = reflectBoundary;
+#endif // BOUNDARY_CONDITIONS
+
+ depth = maxDepth;
+ res = BinaryNode::CumulativeCenterCount( depth );
+ res2 = (1<<(depth+1))+1;
+ baseFunctions = new PPolynomial[res];
+ // Scale the function so that it has:
+ // 0] Value 1 at 0
+ // 1] Integral equal to 1
+ // 2] Square integral equal to 1
+ switch( normalize )
+ {
+ case 2:
+ baseFunction=F/sqrt((F*F).integral(F.polys[0].start,F.polys[F.polyCount-1].start));
+ break;
+ case 1:
+ baseFunction=F/F.integral(F.polys[0].start,F.polys[F.polyCount-1].start);
+ break;
+ default:
+ baseFunction=F/F(0);
+ }
+ dBaseFunction = baseFunction.derivative();
+#if BOUNDARY_CONDITIONS
+ leftBaseFunction = baseFunction + baseFunction.shift( -1 );
+ rightBaseFunction = baseFunction + baseFunction.shift( 1 );
+ dLeftBaseFunction = leftBaseFunction.derivative();
+ dRightBaseFunction = rightBaseFunction.derivative();
+#endif // BOUNDARY_CONDITIONS
+ double c1,w1;
+ for( int i=0 ; i::CenterAndWidth( i , c1 , w1 );
+#if BOUNDARY_CONDITIONS
+ if( reflectBoundary )
+ {
+ int d , off;
+ BinaryNode< double >::DepthAndOffset( i , d , off );
+ if ( off==0 ) baseFunctions[i] = leftBaseFunction.scale( w1 ).shift( c1 );
+ else if( off==((1<
+void FunctionData::setDotTables( const int& flags )
+{
+ clearDotTables( flags );
+ int size;
+ size = ( res*res + res )>>1;
+ if( flags & DOT_FLAG )
+ {
+ dotTable = new Real[size];
+ memset( dotTable , 0 , sizeof(Real)*size );
+ }
+ if( flags & D_DOT_FLAG )
+ {
+ dDotTable = new Real[size];
+ memset( dDotTable , 0 , sizeof(Real)*size );
+ }
+ if( flags & D2_DOT_FLAG )
+ {
+ d2DotTable = new Real[size];
+ memset( d2DotTable , 0 , sizeof(Real)*size );
+ }
+ double t1 , t2;
+ t1 = baseFunction.polys[0].start;
+ t2 = baseFunction.polys[baseFunction.polyCount-1].start;
+ for( int i=0 ; i::CenterAndWidth( i , c1 , w1 );
+#if BOUNDARY_CONDITIONS
+ int d1 , d2 , off1 , off2;
+ BinaryNode< double >::DepthAndOffset( i , d1 , off1 );
+ int boundary1 = 0;
+ if ( reflectBoundary && off1==0 ) boundary1 = -1;
+ else if( reflectBoundary && off1==( (1<::CenterAndWidth( j , c2 , w2 );
+#if BOUNDARY_CONDITIONS
+ BinaryNode< double >::DepthAndOffset( j , d2 , off2 );
+ int boundary2 = 0;
+ if ( reflectBoundary && off2==0 ) boundary2 = -1;
+ else if( reflectBoundary && off2==( (1<1 ) start = 1;
+ if( end <0 ) end = 0;
+ if( end >1 ) end = 1;
+ }
+#endif // BOUNDARY_CONDITIONS
+
+ if( start< start1 ) start = start1;
+ if( end > end1 ) end = end1;
+ if( start>= end ) continue;
+
+#if BOUNDARY_CONDITIONS
+ Real dot = dotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
+#else // !BOUNDARY_CONDITIONS
+ Real dot = dotProduct( c1 , w1 , c2 , w2 );
+#endif // BOUNDARY_CONDITIONS
+ if( fabs(dot)<1e-15 ) continue;
+ if( flags & DOT_FLAG ) dotTable[idx]=dot;
+ if( useDotRatios )
+ {
+#if BOUNDARY_CONDITIONS
+ if( flags & D_DOT_FLAG ) dDotTable[idx] = -dDotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 ) / dot;
+ if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 ) / dot;
+#else // !BOUNDARY_CONDITIONS
+ if( flags & D_DOT_FLAG ) dDotTable[idx] = -dDotProduct(c1,w1,c2,w2)/dot;
+ if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,w1,c2,w2)/dot;
+#endif // BOUNDARY_CONDITIONS
+ }
+ else
+ {
+#if BOUNDARY_CONDITIONS
+ if( flags & D_DOT_FLAG ) dDotTable[idx] = dDotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
+ if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
+#else // !BOUNDARY_CONDTIONS
+ if( flags & D_DOT_FLAG ) dDotTable[idx] = dDotProduct(c1,w1,c2,w2);
+ if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,w1,c2,w2);
+#endif // BOUNDARY_CONDITIONS
+ }
+ }
+ }
+}
+template
+void FunctionData::clearDotTables( const int& flags )
+{
+ if((flags & DOT_FLAG) && dotTable)
+ {
+ delete[] dotTable;
+ dotTable=NULL;
+ }
+ if((flags & D_DOT_FLAG) && dDotTable)
+ {
+ delete[] dDotTable;
+ dDotTable=NULL;
+ }
+ if((flags & D2_DOT_FLAG) && d2DotTable)
+ {
+ delete[] d2DotTable;
+ d2DotTable=NULL;
+ }
+}
+template
+void FunctionData::setValueTables( const int& flags , const double& smooth )
+{
+ clearValueTables();
+ if( flags & VALUE_FLAG ) valueTables = new Real[res*res2];
+ if( flags & D_VALUE_FLAG ) dValueTables = new Real[res*res2];
+ PPolynomial function;
+ PPolynomial dFunction;
+ for( int i=0 ; i0)
+ {
+ function=baseFunctions[i].MovingAverage(smooth);
+ dFunction=baseFunctions[i].derivative().MovingAverage(smooth);
+ }
+ else
+ {
+ function=baseFunctions[i];
+ dFunction=baseFunctions[i].derivative();
+ }
+ for( int j=0 ; j
+void FunctionData::setValueTables(const int& flags,const double& valueSmooth,const double& normalSmooth){
+ clearValueTables();
+ if(flags & VALUE_FLAG){ valueTables=new Real[res*res2];}
+ if(flags & D_VALUE_FLAG){dValueTables=new Real[res*res2];}
+ PPolynomial function;
+ PPolynomial dFunction;
+ for(int i=0;i0) { function=baseFunctions[i].MovingAverage(valueSmooth);}
+ else { function=baseFunctions[i];}
+ if(normalSmooth>0) {dFunction=baseFunctions[i].derivative().MovingAverage(normalSmooth);}
+ else {dFunction=baseFunctions[i].derivative();}
+
+ for(int j=0;j
+void FunctionData::clearValueTables(void){
+ if( valueTables){delete[] valueTables;}
+ if(dValueTables){delete[] dValueTables;}
+ valueTables=dValueTables=NULL;
+}
+
+#if BOUNDARY_CONDITIONS
+template
+Real FunctionData::dotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
+{
+ const PPolynomial< Degree > *b1 , *b2;
+ if ( boundary1==-1 ) b1 = & leftBaseFunction;
+ else if( boundary1== 0 ) b1 = & baseFunction;
+ else if( boundary1== 1 ) b1 = &rightBaseFunction;
+ if ( boundary2==-1 ) b2 = & leftBaseFunction;
+ else if( boundary2== 0 ) b2 = & baseFunction;
+ else if( boundary2== 1 ) b2 = &rightBaseFunction;
+ double r=fabs( baseFunction.polys[0].start );
+ switch( normalize )
+ {
+ case 2:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
+ case 1:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
+ default:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
+ }
+}
+template
+Real FunctionData::dDotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
+{
+ const PPolynomial< Degree-1 > *b1;
+ const PPolynomial< Degree > *b2;
+ if ( boundary1==-1 ) b1 = & dLeftBaseFunction;
+ else if( boundary1== 0 ) b1 = & dBaseFunction;
+ else if( boundary1== 1 ) b1 = &dRightBaseFunction;
+ if ( boundary2==-1 ) b2 = & leftBaseFunction;
+ else if( boundary2== 0 ) b2 = & baseFunction;
+ else if( boundary2== 1 ) b2 = & rightBaseFunction;
+ double r=fabs(baseFunction.polys[0].start);
+ switch(normalize){
+ case 2:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
+ case 1:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
+ default:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
+ }
+}
+template
+Real FunctionData::d2DotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
+{
+ const PPolynomial< Degree-1 > *b1 , *b2;
+ if ( boundary1==-1 ) b1 = & dLeftBaseFunction;
+ else if( boundary1== 0 ) b1 = & dBaseFunction;
+ else if( boundary1== 1 ) b1 = &dRightBaseFunction;
+ if ( boundary2==-1 ) b2 = & dLeftBaseFunction;
+ else if( boundary2== 0 ) b2 = & dBaseFunction;
+ else if( boundary2== 1 ) b2 = &dRightBaseFunction;
+ double r=fabs(baseFunction.polys[0].start);
+ switch( normalize )
+ {
+ case 2:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
+ case 1:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
+ default:
+ return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
+ }
+}
+#else // !BOUNDARY_CONDITIONS
+template
+Real FunctionData::dotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
+ double r=fabs(baseFunction.polys[0].start);
+ switch( normalize )
+ {
+ case 2:
+ return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
+ case 1:
+ return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
+ default:
+ return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
+ }
+}
+template
+Real FunctionData::dDotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
+ double r=fabs(baseFunction.polys[0].start);
+ switch(normalize){
+ case 2:
+ return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
+ case 1:
+ return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
+ default:
+ return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
+ }
+}
+template
+Real FunctionData::d2DotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
+ double r=fabs(baseFunction.polys[0].start);
+ switch(normalize){
+ case 2:
+ return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
+ case 1:
+ return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
+ default:
+ return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
+ }
+}
+#endif // BOUNDARY_CONDITIONS
+template
+inline int FunctionData::SymmetricIndex( const int& i1 , const int& i2 )
+{
+ if( i1>i2 ) return ((i1*i1+i1)>>1)+i2;
+ else return ((i2*i2+i2)>>1)+i1;
+}
+template
+inline int FunctionData::SymmetricIndex( const int& i1 , const int& i2 , int& index )
+{
+ if( i1>1)+i1;
+ return 1;
+ }
+ else{
+ index = ((i1*i1+i1)>>1)+i2;
+ return 0;
+ }
+}
diff --git a/Src/Geometry.cpp b/Src/Geometry.cpp
new file mode 100644
index 0000000..d6fa728
--- /dev/null
+++ b/Src/Geometry.cpp
@@ -0,0 +1,122 @@
+/*
+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
+#ifdef _WIN32
+#include
+#endif // _WIN32
+
+///////////////////
+// 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 , _fp = fopen( _fileName , "w+b" );
+ else
+ {
+ strcpy( _fileName , "PR_XXXXXX" );
+#ifdef _WIN32
+ _mktemp( _fileName );
+ _fp = fopen( _fileName , "w+b" );
+#else // !_WIN32
+ _fp = fdopen( mkstemp( _fileName ) , "w+b" );
+#endif // _WIN32
+ tempFile = true;
+ }
+ 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/Geometry.h b/Src/Geometry.h
new file mode 100644
index 0000000..867e993
--- /dev/null
+++ b/Src/Geometry.h
@@ -0,0 +1,378 @@
+/*
+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 GEOMETRY_INCLUDED
+#define GEOMETRY_INCLUDED
+
+#include
+#include
+#include
+#include "Hash.h"
+
+template
+Real Random(void);
+
+template< class Real >
+struct Point3D
+{
+ Real coords[3];
+ Point3D( void ) { coords[0] = coords[1] = coords[2] = Real(0); }
+ template< class _Real > Point3D( _Real v0 , _Real v1 , _Real v2 ){ coords[0] = Real(v0) , coords[1] = Real(v1) , coords[2] = Real(v2); }
+ template< class _Real > Point3D( const Point3D< _Real >& p ){ coords[0] = Real( p[0] ) , coords[1] = Real( p[1] ) , coords[2] = Real( p[2] ); }
+ inline Real& operator[] ( int i ) { return coords[i]; }
+ inline const Real& operator[] ( int i ) const { return coords[i]; }
+ inline Point3D operator - ( void ) const { Point3D q ; q.coords[0] = -coords[0] , q.coords[1] = -coords[1] , q.coords[2] = -coords[2] ; return q; }
+
+ template< class _Real > inline Point3D& operator += ( Point3D< _Real > p ){ coords[0] += Real(p.coords[0]) , coords[1] += Real(p.coords[1]) , coords[2] += Real(p.coords[2]) ; return *this; }
+ template< class _Real > inline Point3D operator + ( Point3D< _Real > p ) const { Point3D q ; q.coords[0] = coords[0] + Real(p.coords[0]) , q.coords[1] = coords[1] + Real(p.coords[1]) , q.coords[2] = coords[2] + Real(p.coords[2]) ; return q; }
+ template< class _Real > inline Point3D& operator *= ( _Real r ) { coords[0] *= Real(r) , coords[1] *= Real(r) , coords[2] *= Real(r) ; return *this; }
+ template< class _Real > inline Point3D operator * ( _Real r ) const { Point3D q ; q.coords[0] = coords[0] * Real(r) , q.coords[1] = coords[1] * Real(r) , q.coords[2] = coords[2] * Real(r) ; return q; }
+
+ template< class _Real > inline Point3D& operator -= ( Point3D< _Real > p ){ return ( (*this)+=(-p) ); }
+ template< class _Real > inline Point3D operator - ( Point3D< _Real > p ) const { return (*this)+(-p); }
+ template< class _Real > inline Point3D& operator /= ( _Real r ){ return ( (*this)*=Real(1./r) ); }
+ template< class _Real > inline Point3D operator / ( _Real r ) const { return (*this) * ( Real(1.)/r ); }
+
+ static Real Dot( const Point3D< Real >& p1 , const Point3D< Real >& p2 ){ return p1.coords[0]*p2.coords[0] + p1.coords[1]*p2.coords[1] + p1.coords[2]*p2.coords[2]; }
+ template< class Real1 , class Real2 >
+ static Real Dot( const Point3D< Real1 >& p1 , const Point3D< Real2 >& p2 ){ return Real( p1.coords[0]*p2.coords[0] + p1.coords[1]*p2.coords[1] + p1.coords[2]*p2.coords[2] ); }
+};
+
+template< class Real >
+struct XForm3x3
+{
+ Real coords[3][3];
+ XForm3x3( void ) { for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) coords[i][j] = Real(0.); }
+ static XForm3x3 Identity( void )
+ {
+ XForm3x3 xForm;
+ xForm(0,0) = xForm(1,1) = xForm(2,2) = Real(1.);
+ return xForm;
+ }
+ Real& operator() ( int i , int j ){ return coords[i][j]; }
+ const Real& operator() ( int i , int j ) const { return coords[i][j]; }
+ template< class _Real > Point3D< _Real > operator * ( const Point3D< _Real >& p ) const
+ {
+ Point3D< _Real > q;
+ for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) q[i] += _Real( coords[j][i] * p[j] );
+ return q;
+ }
+ XForm3x3 operator * ( const XForm3x3& m ) const
+ {
+ XForm3x3 n;
+ for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) for( int k=0 ; k<3 ; k++ ) n.coords[i][j] += m.coords[i][k]*coords[k][j];
+ return n;
+ }
+ XForm3x3 transpose( void ) const
+ {
+ XForm3x3 xForm;
+ for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ ) xForm( i , j ) = coords[j][i];
+ return xForm;
+ }
+ Real subDeterminant( int i , int j ) const
+ {
+ int i1 = (i+1)%3 , i2 = (i+2)%3;
+ int j1 = (j+1)%3 , j2 = (j+2)%3;
+ return coords[i1][j1] * coords[i2][j2] - coords[i1][j2] * coords[i2][j1];
+ }
+ Real determinant( void ) const { return coords[0][0] * subDeterminant( 0 , 0 ) + coords[1][0] * subDeterminant( 1 , 0 ) + coords[2][0] * subDeterminant( 2 , 0 ); }
+ XForm3x3 inverse( void ) const
+ {
+ XForm3x3 xForm;
+ Real d = determinant();
+ for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ;j++ ) xForm.coords[j][i] = subDeterminant( i , j ) / d;
+ return xForm;
+ }
+};
+
+template< class Real >
+struct XForm4x4
+{
+ Real coords[4][4];
+ XForm4x4( void ) { for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ ) coords[i][j] = Real(0.); }
+ static XForm4x4 Identity( void )
+ {
+ XForm4x4 xForm;
+ xForm(0,0) = xForm(1,1) = xForm(2,2) = xForm(3,3) = Real(1.);
+ return xForm;
+ }
+ Real& operator() ( int i , int j ){ return coords[i][j]; }
+ const Real& operator() ( int i , int j ) const { return coords[i][j]; }
+ template< class _Real > Point3D< _Real > operator * ( const Point3D< _Real >& p ) const
+ {
+ Point3D< _Real > q;
+ for( int i=0 ; i<3 ; i++ )
+ {
+ for( int j=0 ; j<3 ; j++ ) q[i] += (_Real)( coords[j][i] * p[j] );
+ q[i] += (_Real)coords[3][i];
+ }
+ return q;
+ }
+ XForm4x4 operator * ( const XForm4x4& m ) const
+ {
+ XForm4x4 n;
+ for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ ) for( int k=0 ; k<4 ; k++ ) n.coords[i][j] += m.coords[i][k]*coords[k][j];
+ return n;
+ }
+ XForm4x4 transpose( void ) const
+ {
+ XForm4x4 xForm;
+ for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ ) xForm( i , j ) = coords[j][i];
+ return xForm;
+ }
+ Real subDeterminant( int i , int j ) const
+ {
+ XForm3x3< Real > xForm;
+ int ii[] = { (i+1)%4 , (i+2)%4 , (i+3)%4 } , jj[] = { (j+1)%4 , (j+2)%4 , (j+3)%4 };
+ for( int _i=0 ; _i<3 ; _i++ ) for( int _j=0 ; _j<3 ; _j++ ) xForm( _i , _j ) = coords[ ii[_i] ][ jj[_j] ];
+ return xForm.determinant();
+ }
+ Real determinant( void ) const { return coords[0][0] * subDeterminant( 0 , 0 ) - coords[1][0] * subDeterminant( 1 , 0 ) + coords[2][0] * subDeterminant( 2 , 0 ) - coords[3][0] * subDeterminant( 3 , 0 ); }
+ XForm4x4 inverse( void ) const
+ {
+ XForm4x4 xForm;
+ Real d = determinant();
+ for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ;j++ )
+ if( (i+j)%2==0 ) xForm.coords[j][i] = subDeterminant( i , j ) / d;
+ else xForm.coords[j][i] = -subDeterminant( i , j ) / d;
+ return xForm;
+ }
+};
+
+
+template
+Point3D RandomBallPoint(void);
+
+template
+Point3D RandomSpherePoint(void);
+
+template
+double Length(const Point3D& p);
+
+template
+double SquareLength(const Point3D& p);
+
+template
+double Distance(const Point3D& p1,const Point3D& p2);
+
+template
+double SquareDistance(const Point3D& p1,const Point3D& p2);
+
+template
+void CrossProduct(const Point3D& p1,const Point3D& p2,Point3D& p);
+
+
+class Edge{
+public:
+ double p[2][2];
+ double Length(void) const{
+ double d[2];
+ d[0]=p[0][0]-p[1][0];
+ d[1]=p[0][1]-p[1][1];
+
+ return sqrt(d[0]*d[0]+d[1]*d[1]);
+ }
+};
+class Triangle{
+public:
+ double p[3][3];
+ double Area(void) const{
+ double v1[3] , v2[3] , v[3];
+ for( int d=0 ; d<3 ; d++ )
+ {
+ v1[d] = p[1][d] - p[0][d];
+ v2[d] = p[2][d] - p[0][d];
+ }
+ v[0] = v1[1]*v2[2] - v1[2]*v2[1];
+ v[1] = -v1[0]*v2[2] + v1[2]*v2[0];
+ v[2] = v1[0]*v2[1] - v1[1]*v2[0];
+ return sqrt( v[0]*v[0] + v[1]*v[1] + v[2]*v[2] ) / 2;
+ }
+ double AspectRatio(void) const{
+ double d=0;
+ int i,j;
+ for(i=0;i<3;i++){
+ for(i=0;i<3;i++)
+ for(j=0;j<3;j++){d+=(p[(i+1)%3][j]-p[i][j])*(p[(i+1)%3][j]-p[i][j]);}
+ }
+ return Area()/d;
+ }
+
+};
+class CoredPointIndex
+{
+public:
+ int index;
+ char inCore;
+
+ int operator == (const CoredPointIndex& cpi) const {return (index==cpi.index) && (inCore==cpi.inCore);};
+ int operator != (const CoredPointIndex& cpi) const {return (index!=cpi.index) || (inCore!=cpi.inCore);};
+};
+class EdgeIndex{
+public:
+ int idx[2];
+};
+class CoredEdgeIndex
+{
+public:
+ CoredPointIndex idx[2];
+};
+class TriangleIndex{
+public:
+ int idx[3];
+};
+
+class TriangulationEdge
+{
+public:
+ TriangulationEdge(void);
+ int pIndex[2];
+ int tIndex[2];
+};
+
+class TriangulationTriangle
+{
+public:
+ TriangulationTriangle(void);
+ int eIndex[3];
+};
+
+template
+class Triangulation
+{
+public:
+
+ std::vector > points;
+ std::vector edges;
+ std::vector triangles;
+
+ int factor( int tIndex,int& p1,int& p2,int& p3);
+ double area(void);
+ double area( int tIndex );
+ double area( int p1 , int p2 , int p3 );
+ int flipMinimize( int eIndex);
+ int addTriangle( int p1 , int p2 , int p3 );
+
+protected:
+ hash_map edgeMap;
+ static long long EdgeIndex( int p1 , int p2 );
+ double area(const Triangle& t);
+};
+
+
+template
+void EdgeCollapse(const Real& edgeRatio,std::vector& triangles,std::vector< Point3D >& positions,std::vector >* normals);
+template
+void TriangleCollapse(const Real& edgeRatio,std::vector& triangles,std::vector >& positions,std::vector >* normals);
+
+struct CoredVertexIndex
+{
+ int idx;
+ bool inCore;
+};
+template< class Vertex >
+class CoredMeshData
+{
+public:
+ std::vector< Vertex > inCorePoints;
+ virtual void resetIterator( void ) = 0;
+
+ virtual int addOutOfCorePoint( const Vertex& p ) = 0;
+ virtual int addOutOfCorePoint_s( const Vertex& p ) = 0;
+ virtual int addPolygon_s( const std::vector< CoredVertexIndex >& vertices ) = 0;
+ virtual int addPolygon_s( const std::vector< int >& vertices ) = 0;
+
+ virtual int nextOutOfCorePoint( Vertex& p )=0;
+ virtual int nextPolygon( std::vector< CoredVertexIndex >& vertices ) = 0;
+
+ virtual int outOfCorePointCount(void)=0;
+ virtual int polygonCount( void ) = 0;
+};
+
+template< class Vertex >
+class CoredVectorMeshData : public CoredMeshData< Vertex >
+{
+ std::vector< Vertex > oocPoints;
+ std::vector< std::vector< int > > polygons;
+ int polygonIndex;
+ int oocPointIndex;
+public:
+ CoredVectorMeshData(void);
+
+ void resetIterator(void);
+
+ int addOutOfCorePoint( const Vertex& p );
+ int addOutOfCorePoint_s( const Vertex& p );
+ int addPolygon_s( const std::vector< CoredVertexIndex >& vertices );
+ int addPolygon_s( const std::vector< int >& vertices );
+
+ int nextOutOfCorePoint( Vertex& p );
+ int nextPolygon( std::vector< CoredVertexIndex >& vertices );
+
+ int outOfCorePointCount(void);
+ int polygonCount( void );
+};
+class BufferedReadWriteFile
+{
+ bool tempFile;
+ FILE* _fp;
+ char *_buffer , _fileName[1024];
+ size_t _bufferIndex , _bufferSize;
+public:
+ BufferedReadWriteFile( char* fileName=NULL , int bufferSize=(1<<20) );
+ ~BufferedReadWriteFile( void );
+ bool write( const void* data , size_t size );
+ bool read ( void* data , size_t size );
+ void reset( void );
+};
+template< class Vertex >
+class CoredFileMeshData : public CoredMeshData< Vertex >
+{
+ char pointFileName[1024] , polygonFileName[1024];
+ BufferedReadWriteFile *oocPointFile , *polygonFile;
+ int oocPoints , polygons;
+public:
+ CoredFileMeshData( void );
+ ~CoredFileMeshData( void );
+
+ void resetIterator( void );
+
+ int addOutOfCorePoint( const Vertex& p );
+ int addOutOfCorePoint_s( const Vertex& p );
+ int addPolygon_s( const std::vector< CoredVertexIndex >& vertices );
+ int addPolygon_s( const std::vector< int >& vertices );
+
+ int nextOutOfCorePoint( Vertex& p );
+ int nextPolygon( std::vector< CoredVertexIndex >& vertices );
+
+ int outOfCorePointCount( void );
+ int polygonCount( void );
+};
+#include "Geometry.inl"
+
+#endif // GEOMETRY_INCLUDED
diff --git a/Src/Geometry.inl b/Src/Geometry.inl
new file mode 100644
index 0000000..97be936
--- /dev/null
+++ b/Src/Geometry.inl
@@ -0,0 +1,590 @@
+/*
+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
+
+template
+Real Random(void){return Real(rand())/RAND_MAX;}
+
+template
+Point3D RandomBallPoint(void){
+ Point3D p;
+ while(1){
+ p.coords[0]=Real(1.0-2.0*Random());
+ p.coords[1]=Real(1.0-2.0*Random());
+ p.coords[2]=Real(1.0-2.0*Random());
+ double l=SquareLength(p);
+ if(l<=1){return p;}
+ }
+}
+template
+Point3D RandomSpherePoint(void){
+ Point3D p=RandomBallPoint();
+ Real l=Real(Length(p));
+ p.coords[0]/=l;
+ p.coords[1]/=l;
+ p.coords[2]/=l;
+ return p;
+}
+
+template
+double SquareLength(const Point3D& p){return p.coords[0]*p.coords[0]+p.coords[1]*p.coords[1]+p.coords[2]*p.coords[2];}
+
+template
+double Length(const Point3D& p){return sqrt(SquareLength(p));}
+
+template
+double SquareDistance(const Point3D& p1,const Point3D& p2){
+ return (p1.coords[0]-p2.coords[0])*(p1.coords[0]-p2.coords[0])+(p1.coords[1]-p2.coords[1])*(p1.coords[1]-p2.coords[1])+(p1.coords[2]-p2.coords[2])*(p1.coords[2]-p2.coords[2]);
+}
+
+template
+double Distance(const Point3D& p1,const Point3D& p2){return sqrt(SquareDistance(p1,p2));}
+
+template
+void CrossProduct(const Point3D& p1,const Point3D& p2,Point3D& p){
+ p.coords[0]= p1.coords[1]*p2.coords[2]-p1.coords[2]*p2.coords[1];
+ p.coords[1]=-p1.coords[0]*p2.coords[2]+p1.coords[2]*p2.coords[0];
+ p.coords[2]= p1.coords[0]*p2.coords[1]-p1.coords[1]*p2.coords[0];
+}
+template
+void EdgeCollapse(const Real& edgeRatio,std::vector& triangles,std::vector< Point3D >& positions,std::vector< Point3D >* normals){
+ int i,j,*remapTable,*pointCount,idx[3];
+ Point3D p[3],q[2],c;
+ double d[3],a;
+ double Ratio=12.0/sqrt(3.0); // (Sum of Squares Length / Area) for and equilateral triangle
+
+ remapTable=new int[positions.size()];
+ pointCount=new int[positions.size()];
+ for(i=0;i=0;i--){
+ for(j=0;j<3;j++){
+ idx[j]=triangles[i].idx[j];
+ while(remapTable[idx[j]] a*Ratio){
+ // Find the smallest edge
+ j=0;
+ if(d[1]=0;i--){
+ for(j=0;j<3;j++){
+ idx[j]=triangles[i].idx[j];
+ while(remapTable[idx[j]]
+void TriangleCollapse(const Real& edgeRatio,std::vector& triangles,std::vector< Point3D >& positions,std::vector< Point3D >* normals){
+ int i,j,*remapTable,*pointCount,idx[3];
+ Point3D p[3],q[2],c;
+ double d[3],a;
+ double Ratio=12.0/sqrt(3.0); // (Sum of Squares Length / Area) for and equilateral triangle
+
+ remapTable=new int[positions.size()];
+ pointCount=new int[positions.size()];
+ for(i=0;i=0;i--){
+ for(j=0;j<3;j++){
+ idx[j]=triangles[i].idx[j];
+ while(remapTable[idx[j]] a*Ratio){
+ // Find the smallest edge
+ j=0;
+ if(d[1]=0;i--){
+ for(j=0;j<3;j++){
+ idx[j]=triangles[i].idx[j];
+ while(remapTable[idx[j]]
+long long Triangulation::EdgeIndex( int p1 , int p2 )
+{
+ if(p1>p2) {return ((long long)(p1)<<32) | ((long long)(p2));}
+ else {return ((long long)(p2)<<32) | ((long long)(p1));}
+}
+
+template
+int Triangulation::factor(int tIndex,int& p1,int& p2,int & p3){
+ if(triangles[tIndex].eIndex[0]<0 || triangles[tIndex].eIndex[1]<0 || triangles[tIndex].eIndex[2]<0){return 0;}
+ if(edges[triangles[tIndex].eIndex[0]].tIndex[0]==tIndex){p1=edges[triangles[tIndex].eIndex[0]].pIndex[0];}
+ else {p1=edges[triangles[tIndex].eIndex[0]].pIndex[1];}
+ if(edges[triangles[tIndex].eIndex[1]].tIndex[0]==tIndex){p2=edges[triangles[tIndex].eIndex[1]].pIndex[0];}
+ else {p2=edges[triangles[tIndex].eIndex[1]].pIndex[1];}
+ if(edges[triangles[tIndex].eIndex[2]].tIndex[0]==tIndex){p3=edges[triangles[tIndex].eIndex[2]].pIndex[0];}
+ else {p3=edges[triangles[tIndex].eIndex[2]].pIndex[1];}
+ return 1;
+}
+template
+double Triangulation::area(int p1,int p2,int p3){
+ Point3D q1,q2,q;
+ for(int i=0;i<3;i++){
+ q1.coords[i]=points[p2].coords[i]-points[p1].coords[i];
+ q2.coords[i]=points[p3].coords[i]-points[p1].coords[i];
+ }
+ CrossProduct(q1,q2,q);
+ return Length(q);
+}
+template
+double Triangulation::area(int tIndex){
+ int p1,p2,p3;
+ factor(tIndex,p1,p2,p3);
+ return area(p1,p2,p3);
+}
+template
+double Triangulation::area(void){
+ double a=0;
+ for(int i=0;i
+int Triangulation::addTriangle(int p1,int p2,int p3){
+ hash_map::iterator iter;
+ int tIdx,eIdx,p[3];
+ p[0]=p1;
+ p[1]=p2;
+ p[2]=p3;
+ triangles.push_back(TriangulationTriangle());
+ tIdx=int(triangles.size())-1;
+
+ for(int i=0;i<3;i++)
+ {
+ long long e = EdgeIndex(p[i],p[(i+1)%3]);
+ iter=edgeMap.find(e);
+ if(iter==edgeMap.end())
+ {
+ TriangulationEdge edge;
+ edge.pIndex[0]=p[i];
+ edge.pIndex[1]=p[(i+1)%3];
+ edges.push_back(edge);
+ eIdx=int(edges.size())-1;
+ edgeMap[e]=eIdx;
+ edges[eIdx].tIndex[0]=tIdx;
+ }
+ else{
+ eIdx=edgeMap[e];
+ if(edges[eIdx].pIndex[0]==p[i]){
+ if(edges[eIdx].tIndex[0]<0){edges[eIdx].tIndex[0]=tIdx;}
+ else{printf("Edge Triangle in use 1\n");return 0;}
+ }
+ else{
+ if(edges[eIdx].tIndex[1]<0){edges[eIdx].tIndex[1]=tIdx;}
+ else{printf("Edge Triangle in use 2\n");return 0;}
+ }
+
+ }
+ triangles[tIdx].eIndex[i]=eIdx;
+ }
+ return tIdx;
+}
+template
+int Triangulation::flipMinimize(int eIndex){
+ double oldArea,newArea;
+ int oldP[3],oldQ[3],newP[3],newQ[3];
+ TriangulationEdge newEdge;
+
+ if(edges[eIndex].tIndex[0]<0 || edges[eIndex].tIndex[1]<0){return 0;}
+
+ if(!factor(edges[eIndex].tIndex[0],oldP[0],oldP[1],oldP[2])){return 0;}
+ if(!factor(edges[eIndex].tIndex[1],oldQ[0],oldQ[1],oldQ[2])){return 0;}
+
+ oldArea=area(oldP[0],oldP[1],oldP[2])+area(oldQ[0],oldQ[1],oldQ[2]);
+ int idxP,idxQ;
+ for(idxP=0;idxP<3;idxP++){
+ int i;
+ for(i=0;i<3;i++){if(oldP[idxP]==oldQ[i]){break;}}
+ if(i==3){break;}
+ }
+ for(idxQ=0;idxQ<3;idxQ++){
+ int i;
+ for(i=0;i<3;i++){if(oldP[i]==oldQ[idxQ]){break;}}
+ if(i==3){break;}
+ }
+ if(idxP==3 || idxQ==3){return 0;}
+ newP[0]=oldP[idxP];
+ newP[1]=oldP[(idxP+1)%3];
+ newP[2]=oldQ[idxQ];
+ newQ[0]=oldQ[idxQ];
+ newQ[1]=oldP[(idxP+2)%3];
+ newQ[2]=oldP[idxP];
+
+ newArea=area(newP[0],newP[1],newP[2])+area(newQ[0],newQ[1],newQ[2]);
+ if(oldArea<=newArea){return 0;}
+
+ // Remove the entry in the hash_table for the old edge
+ edgeMap.erase(EdgeIndex(edges[eIndex].pIndex[0],edges[eIndex].pIndex[1]));
+ // Set the new edge so that the zero-side is newQ
+ edges[eIndex].pIndex[0]=newP[0];
+ edges[eIndex].pIndex[1]=newQ[0];
+ // Insert the entry into the hash_table for the new edge
+ edgeMap[EdgeIndex(newP[0],newQ[0])]=eIndex;
+ // Update the triangle information
+ for(int i=0;i<3;i++){
+ int idx;
+ idx=edgeMap[EdgeIndex(newQ[i],newQ[(i+1)%3])];
+ triangles[edges[eIndex].tIndex[0]].eIndex[i]=idx;
+ if(idx!=eIndex){
+ if(edges[idx].tIndex[0]==edges[eIndex].tIndex[1]){edges[idx].tIndex[0]=edges[eIndex].tIndex[0];}
+ if(edges[idx].tIndex[1]==edges[eIndex].tIndex[1]){edges[idx].tIndex[1]=edges[eIndex].tIndex[0];}
+ }
+
+ idx=edgeMap[EdgeIndex(newP[i],newP[(i+1)%3])];
+ triangles[edges[eIndex].tIndex[1]].eIndex[i]=idx;
+ if(idx!=eIndex){
+ if(edges[idx].tIndex[0]==edges[eIndex].tIndex[0]){edges[idx].tIndex[0]=edges[eIndex].tIndex[1];}
+ if(edges[idx].tIndex[1]==edges[eIndex].tIndex[0]){edges[idx].tIndex[1]=edges[eIndex].tIndex[1];}
+ }
+ }
+ return 1;
+}
+/////////////////////////
+// CoredVectorMeshData //
+/////////////////////////
+template< class Vertex >
+CoredVectorMeshData< Vertex >::CoredVectorMeshData( void ) { oocPointIndex = polygonIndex = 0; }
+template< class Vertex >
+void CoredVectorMeshData< Vertex >::resetIterator ( void ) { oocPointIndex = polygonIndex = 0; }
+template< class Vertex >
+int CoredVectorMeshData< Vertex >::addOutOfCorePoint( const Vertex& p )
+{
+ oocPoints.push_back(p);
+ return int(oocPoints.size())-1;
+}
+template< class Vertex >
+int CoredVectorMeshData< Vertex >::addOutOfCorePoint_s( const Vertex& p )
+{
+ size_t sz;
+#pragma omp critical (CoredVectorMeshData_addOutOfCorePoint_s )
+ {
+ sz = oocPoints.size();
+ oocPoints.push_back(p);
+ }
+ return (int)sz;
+}
+template< class Vertex >
+int CoredVectorMeshData< Vertex >::addPolygon_s( const std::vector< int >& polygon )
+{
+ size_t sz;
+#pragma omp critical (CoredVectorMeshData_addPolygon_s)
+ {
+ sz = polygon.size();
+ polygons.push_back( polygon );
+ }
+ return (int)sz;
+}
+template< class Vertex >
+int CoredVectorMeshData< Vertex >::addPolygon_s( const std::vector< CoredVertexIndex >& vertices )
+{
+ std::vector< int > polygon( vertices.size() );
+ for( int i=0 ; i<(int)vertices.size() ; i++ )
+ if( vertices[i].inCore ) polygon[i] = vertices[i].idx;
+ else polygon[i] = -vertices[i].idx-1;
+ return addPolygon_s( polygon );
+}
+template< class Vertex >
+int CoredVectorMeshData< Vertex >::nextOutOfCorePoint( Vertex& p )
+{
+ if( oocPointIndex