diff --git a/.gitignore b/.gitignore
index 3fb3f818..73ae1738 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,3 +5,5 @@
/Bin/
*.opensdf
*.bat
+*.zip
+.vs/*
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 00000000..5f806177
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,50 @@
+cmake_minimum_required(VERSION 3.7)
+
+set(CMAKE_DISABLE_SOURCE_CHANGES ON)
+set(CMAKE_DISABLE_IN_SOURCE_BUILD ON)
+set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON)
+
+project(PoissonRecon LANGUAGES CXX)
+
+if(MSVC)
+ set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
+ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
+ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
+ set_property(GLOBAL PROPERTY USE_FOLDERS ON)
+endif()
+
+if((CMAKE_CXX_COMPILER_ID MATCHES "Clang") OR CMAKE_COMPILER_IS_GNUCC
+ OR CMAKE_COMPILER_IS_GNUCXX)
+ if(DEBUG)
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -Winline")
+ else()
+ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -DNDEBUG")
+ endif()
+endif()
+
+message(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
+
+if(NOT TARGET uninstall)
+ configure_file(
+ ${CMAKE_CURRENT_SOURCE_DIR}/cmake_uninstall.cmake.in
+ ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake IMMEDIATE @ONLY)
+
+ add_custom_target(
+ uninstall COMMAND ${CMAKE_COMMAND} -P
+ ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake)
+endif()
+
+export(TARGETS FILE ${CMAKE_BINARY_DIR}/${CMAKE_PROJECT_NAME}Targets.cmake)
+
+file(WRITE ${CMAKE_BINARY_DIR}/${CMAKE_PROJECT_NAME}Config.cmake
+ "include(CMakeFindDependencyMacro)\n"
+ "include(\${CMAKE_CURRENT_LIST_DIR}/${CMAKE_PROJECT_NAME}Targets.cmake)\n")
+
+add_subdirectory(modules)
+add_subdirectory(apps)
+
+install(FILES ${CMAKE_BINARY_DIR}/${CMAKE_PROJECT_NAME}Config.cmake
+ DESTINATION lib/cmake/${CMAKE_PROJECT_NAME})
+
+install(EXPORT ${CMAKE_PROJECT_NAME}Targets
+ DESTINATION lib/cmake/${CMAKE_PROJECT_NAME})
diff --git a/Makefile b/Makefile
deleted file mode 100644
index f71afaf5..00000000
--- a/Makefile
+++ /dev/null
@@ -1,71 +0,0 @@
-PR_TARGET=PoissonRecon
-SR_TARGET=SSDRecon
-ST_TARGET=SurfaceTrimmer
-PR_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cpp PoissonRecon.cpp
-SR_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cpp SSDRecon.cpp
-ST_SOURCE=CmdLineParser.cpp Factor.cpp Geometry.cpp MarchingCubes.cpp PlyFile.cpp SurfaceTrimmer.cpp
-
-CFLAGS += -fopenmp -Wno-deprecated -Wno-write-strings -std=c++11
-LFLAGS += -lgomp -lstdc++
-
-CFLAGS_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))))
-SR_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(SR_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)$(SR_TARGET)
-all: $(BIN)$(ST_TARGET)
-
-release: CFLAGS += $(CFLAGS_RELEASE)
-release: LFLAGS += $(LFLAGS_RELEASE)
-release: $(BIN)
-release: $(BIN)$(PR_TARGET)
-release: $(BIN)$(SR_TARGET)
-release: $(BIN)$(ST_TARGET)
-
-clean:
- rm -f $(BIN)$(PR_TARGET)
- rm -f $(BIN)$(SR_TARGET)
- rm -f $(BIN)$(ST_TARGET)
- rm -f $(PR_OBJECTS)
- rm -f $(SR_OBJECTS)
- rm -f $(ST_OBJECTS)
-
-$(BIN):
- $(MD) -p $(BIN)
-
-$(BIN)$(PR_TARGET): $(PR_OBJECTS)
- $(CXX) -o $@ $(PR_OBJECTS) $(LFLAGS)
-
-$(BIN)$(SR_TARGET): $(SR_OBJECTS)
- $(CXX) -o $@ $(SR_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
deleted file mode 100644
index ed7a3b4f..00000000
--- a/PoissonRecon.sln
+++ /dev/null
@@ -1,34 +0,0 @@
-
-Microsoft Visual Studio Solution File, Format Version 12.00
-# Visual Studio 14
-VisualStudioVersion = 14.0.25420.1
-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
-Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "SSDRecon", "SSDRecon.vcxproj", "{7838CA1E-8A39-4A2B-AC3D-3E25FEAEA2D6}"
-EndProject
-Global
- GlobalSection(SolutionConfigurationPlatforms) = preSolution
- Debug|x64 = Debug|x64
- Release|x64 = Release|x64
- EndGlobalSection
- GlobalSection(ProjectConfigurationPlatforms) = postSolution
- {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|x64.ActiveCfg = Release|x64
- {46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|x64.Build.0 = Release|x64
- {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|x64.ActiveCfg = Release|x64
- {99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|x64.Build.0 = Release|x64
- {7838CA1E-8A39-4A2B-AC3D-3E25FEAEA2D6}.Debug|x64.ActiveCfg = Debug|x64
- {7838CA1E-8A39-4A2B-AC3D-3E25FEAEA2D6}.Debug|x64.Build.0 = Debug|x64
- {7838CA1E-8A39-4A2B-AC3D-3E25FEAEA2D6}.Release|x64.ActiveCfg = Release|x64
- {7838CA1E-8A39-4A2B-AC3D-3E25FEAEA2D6}.Release|x64.Build.0 = Release|x64
- EndGlobalSection
- GlobalSection(SolutionProperties) = preSolution
- HideSolutionNode = FALSE
- EndGlobalSection
-EndGlobal
diff --git a/PoissonRecon.vcxproj b/PoissonRecon.vcxproj
deleted file mode 100644
index e5d1740e..00000000
--- a/PoissonRecon.vcxproj
+++ /dev/null
@@ -1,235 +0,0 @@
-
-
-
-
- Debug
- Win32
-
-
- Debug
- x64
-
-
- Release
- Win32
-
-
- Release
- x64
-
-
-
- PoissonRecon
- {46F87D0E-C53A-4F95-AB48-A5DBA8014340}
- PoissonRecon
- Win32Proj
- 8.1
-
-
-
- Application
- MultiByte
- true
- v140
-
-
- Application
- MultiByte
- v140
-
-
- Application
- MultiByte
- true
- v140
-
-
- Application
- MultiByte
- v140
- true
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- <_ProjectFileVersion>10.0.30319.1
- $(SolutionDir)\Bin\$(Platform)\$(Configuration)\
- $(SolutionDir)\Obj\$(Platform)\$(Configuration)\
- true
- $(SolutionDir)\Bin\$(Platform)\$(Configuration)\
- $(SolutionDir)\Obj\$(TargetName)\$(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
-
-
- 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
deleted file mode 100644
index 16466a46..00000000
--- a/PoissonRecon.vcxproj.filters
+++ /dev/null
@@ -1,143 +0,0 @@
-
-
-
-
- {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
-
-
-
-
- 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
-
-
- Include Files
-
-
- Include Files
-
-
-
\ No newline at end of file
diff --git a/README.md b/README.md
index 482b9742..f7db3246 100644
--- a/README.md
+++ b/README.md
@@ -1,55 +1,84 @@
-
Screened Poisson Surface Reconstruction (and Smoothed Signed Distance Reconstruction) Version 9.01
-
-links
-executables
-usage
-changes
-support
-
-
-LINKS
+Adaptive Multigrid Solvers (Version 12.00)
+
+links
+executables
+usage
+changes
+
+
+
+This code-base was born from the Poisson Surface Reconstruction code. It has evolved to support more general adaptive finite-elements systems:
-Papers:
-[Kazhdan, Bolitho, and Hoppe, 2006] ,
-[Calakli and Taubin, 2011] ,
-[Kazhdan and Hoppe, 2013]
-
-Executables:
-Win64
-Source Code:
-ZIP GitHub
-Older Versions:
-V9.0 ,
-V8.0 ,
-V7.0 ,
-V6.13a ,
-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
+ in spaces of arbitrary dimension,
+ discretized using finite-elements of arbitrary degree,
+ involving arbitrary function derivatives,
+ with both point-wise and integrated constraints.
-
-EXECUTABLES
-
-
-PoissonRecon
-
-
--in <input points >
- This string is the name of the file from which the point set will be read.
+
+LINKS
+
+Papers:
+[Kazhdan, Bolitho, and Hoppe, 2006] ,
+[Agarwala, 2007]
+[Calakli and Taubin, 2011] ,
+[Crane, Weischedel, and Wardetzky, 2013] ,
+[Kazhdan and Hoppe, 2013]
+
+Executables:
+Win64
+Source Code:
+ZIP GitHub
+Older Versions:
+V11.02 ,
+V11.01 ,
+V11.00 ,
+V10.07 ,
+V10.06 ,
+V10.05 ,
+V10.04 ,
+V10.03 ,
+V10.02 ,
+V10.01 ,
+V10.00 ,
+V9.011 ,
+V9.01 ,
+V9.0 ,
+V8.0 ,
+V7.0 ,
+V6.13a ,
+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
+
+
+EXECUTABLES
+
+
+
+
+PoissonRecon :
+Reconstructs a triangle mesh from a set of oriented 3D points by solving a Poisson system (solving a 3D Laplacian system with positional value constraints) [Kazhdan, Bolitho, and Hoppe, 2006] ,
+[Kazhdan and Hoppe, 2013]
+
+--in <input points >
+ This string is the name of the file from which the point set will be read.
If the file extension is .ply , the file should be in
-PLY format, giving the list of oriented
+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 .
@@ -60,119 +89,130 @@ Otherwise, the file should be an ascii file with groups of 6,
white space delimited, numbers: 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.
+ [--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
+ [--tree <output octree and coefficients >]
+ This string is the name of the file to which the the octree and solution coefficients are to be written.
+
+ [--grid <output grid >]
+ This string is the name of the file to which the sampled implicit function will be written.
+The file is written 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 ... bytes corresponding to the (single precision) floating point values
of the implicit function.
- [--degree <B-spline degree >]
- This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
-Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
+ [--degree <B-spline degree >]
+ This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
+Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
The default value for this parameter is 2.
- [--bType <boundary type >]
- This integer specifies the boundary type for the finite elements. Valid values are:
-
- 1 : Free boundary constraints
- 2 : Dirichlet boundary constraints
- 3 : Neumann boundary constraints
-
+ [--bType <boundary type >]
+ This integer specifies the boundary type for the finite elements. Valid values are:
+
+ 1 : Free boundary constraints
+ 2 : Dirichlet boundary constraints
+ 3 : Neumann boundary constraints
+
The default value for this parameter is 3 (Neumann).
- [--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
+ [--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 grid whose resolution is no larger than
+2^d x 2^d x ... Note that since the reconstructor adapts the octree to the
sampling density, the specified reconstruction depth is only an upper bound.
The default value for this parameter is 8.
- [--scale <scale factor >]
- This floating point value specifies the ratio between the diameter of the cube used for reconstruction
+ [--width <finest cell width >]
+ This floating point value specifies the target width of the finest level octree cells.
+This parameter is ignored if the --depth is also specified.
+
+ [--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.
The default value is 1.1.
- [--samplesPerNode <minimum number of samples >]
- This floating point value specifies the minimum number of sample points that should fall within an
+ [--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.
The default value is 1.0.
- [--pointWeight <interpolation weight >]
- This floating point value specifies the importance that interpolation of the point samples
+ [--pointWeight <interpolation weight >]
+ This floating point value specifies the importance that interpolation of the point samples
is given in the formulation of the screened Poisson equation.
The results of the original (unscreened) Poisson Reconstruction can be obtained by setting this value to 0.
-The default value for this parameter is 4.
+The default value for this parameter is twice the B-spline degree.
- [--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.
+ [--iters <Gauss-Seidel iterations per level >]
+ This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the octree hierarchy.
+The default value for this parameter is 8.
- [--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.
+ [--density ]
+ Enabling this flag tells the reconstructor to output the estimated depth values of the iso-surface vertices.
- [--iters <GS iters >]
- This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.
-The default value for this parameter is 8.
+ [--normals ]
+ Enabling this flag tells the reconstructor to output vertex normals, computed from the gradients of the implicit function.
- [--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.
-The default value for this parameter is 0.
+ [--colors ]
+ Enabling this flag tells the reconstructor to read in color values with the input points and extrapolate those to the vertices of the output.
- [--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.
-The default value for this parameter is 5.
+ [--data <pull factor >]
+ If --colors is specified, this floating point value specifies the relative importance
+of finer color estimates over lower ones.
+The default value for this parameter is 32.
- [--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 .
-The default value for this parameter is the value of the --depth parameter.
+ [--confidence <normal confidence exponent >]
+ This floating point value specifies the exponent to be applied to a point's confidence to adjust its weight. (A point's confidence is defined by the magnitude of its normal.)
+The default value for this parameter is 0.
- [--primalVoxel ]
- Enabling this flag when outputing to a voxel file has the reconstructor sample the implicit function at the corners of the grid, rather than the centers of the cells.
+ [--confidenceBias <normal confidence bias exponent >]
+ This floating point value specifies the exponent to be applied to a point's confidence to bias the resolution at which the sample contributes to the linear system. (Points with lower confidence are biased to contribute at coarser resolutions.)
+The default value for this parameter is 0.
- [--color <pull factor >]
- If specified, the reconstruction code assumes that the input is equipped with colors and will extrapolate
-the color values to the vertices of the reconstructed mesh. The floating point value specifies the relative importance
-of finer color estimates over lower ones. (In practice, we have found that a pull factor of 16 works well.)
+ [--primalGrid ]
+ Enabling this flag when outputing to a grid file has the reconstructor sample the implicit function at the corners of the grid, rather than the centers of the cells.
- [--density ]
- Enabling this flag tells the reconstructor to output the estimated depth values of the iso-surface vertices.
+ [--linearFit ]
+ Enabling this flag has the reconstructor use linear interpolation to estimate the positions of iso-vertices.
- [--linearFit ]
- Enabling this flag has the reconstructor use linear interpolation to estimate the positions of iso-vertices.
+ [--polygonMesh ]
+ Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes).
- [--polygonMesh ]
- Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes).
+ [--tempDir <temporary output directory >]
+ This string is the name of the directory to which temporary files will be written.
- [--threads <number of processing threads >]
- This integer specifies the number of threads across which the reconstruction
-algorithm should be parallelized.
+ [--threads <number of processing threads >]
+ This integer specifies the number of threads across which the algorithm should be parallelized.
The default value for this parameter is equal to the numer of (virtual) processors on the executing machine.
- [--verbose ]
- Enabling this flag provides a more verbose description of the running times and memory usages of
-individual components of the surface reconstructor.
-
-
-
+[--maxMemory <maximum memory usage (in GB) >]
+ If positive, this integer value specifies the peak memory utilization for running the reconstruction code (forcing the execution to terminate if the limit is exceeded).
+ [--performance ]
+ Enabling this flag provides running time and peak memory usage at the end of the execution.
-
-
-SSDRecon
-
-
--in <input points >
- This string is the name of the file from which the point set will be read.
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of individual components of the surface reconstructor.
+
+
+
+
+
+
+
+
+
+
+
+SSDRecon :
+Reconstructs a surface mesh from a set of oriented 3D points by solving for a Smooth Signed Distance function (solving a 3D bi-Laplacian system with positional value and gradient constraints) [Calakli and Taubin, 2011]
+
+--in <input points >
+ This string is the name of the file from which the point set will be read.
If the file extension is .ply , the file should be in
-PLY format, giving the list of oriented
+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 .
@@ -183,328 +223,787 @@ Otherwise, the file should be an ascii file with groups of 6,
white space delimited, numbers: 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.
+ [--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
+ [--tree <output octree and coefficients >]
+ This string is the name of the file to which the the octree and solution coefficients are to be written.
+
+ [--grid <output grid >]
+ This string is the name of the file to which the sampled implicit function will be written.
+The file 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 ... bytes corresponding to the (single precision) floating point values
of the implicit function.
- [--degree <B-spline degree >]
- This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
-Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
+ [--degree <B-spline degree >]
+ This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
+Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
The default value for this parameter is 2.
- [--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
+ [--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 grid whose resolution is no larger than
+2^d x 2^d x ... Note that since the reconstructor adapts the octree to the
sampling density, the specified reconstruction depth is only an upper bound.
The default value for this parameter is 8.
- [--scale <scale factor >]
- This floating point value specifies the ratio between the diameter of the cube used for reconstruction
+ [--width <finest cell width >]
+ This floating point value specifies the target width of the finest level octree cells.
+This parameter is ignored if the --depth is also specified.
+
+ [--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.
The default value is 1.1.
- [--samplesPerNode <minimum number of samples >]
- This floating point value specifies the minimum number of sample points that should fall within an
+ [--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.
The default value is 1.0.
- [--valueWeight <zero-crossing interpolation weight >]
- This floating point value specifies the importance that interpolation of the point samples
+ [--valueWeight <zero-crossing interpolation weight >]
+ This floating point value specifies the importance that interpolation of the point samples
is given in the formulation of the screened Smoothed Signed Distance Reconstruction.
-The default value for this parameter is 4.
+The default value for this parameter is 1.
- [--gradientWeight <normal interpolation weight >]
- This floating point value specifies the importance that interpolation of the points' normals
+ [--gradientWeight <normal interpolation weight >]
+ This floating point value specifies the importance that interpolation of the points' normals
is given in the formulation of the screened Smoothed Signed Distance Reconstruction.
-The default value for this parameter is 0.001.
+The default value for this parameter is 1.
- [--biLapWeight <bi-Laplacian weight weight >]
- This floating point value specifies the importance that the bi-Laplacian regularization
+ [--biLapWeight <bi-Laplacian weight weight >]
+ This floating point value specifies the importance that the bi-Laplacian regularization
is given in the formulation of the screened Smoothed Signed Distance Reconstruction.
-The default value for this parameter is 0.00001.
+The default value for this parameter is 1.
- [--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.
+ [--iters <GS iters >]
+ This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.
+The default value for this parameter is 8.
- [--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.
+ [--density ]
+ Enabling this flag tells the reconstructor to output the estimated depth values of the iso-surface vertices.
- [--iters <GS iters >]
- This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.
-The default value for this parameter is 8.
+ [--normals ]
+ Enabling this flag tells the reconstructor to output vertex normals, computed from the gradients of the implicit function.
- [--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.
-The default value for this parameter is 0.
+ [--colors ]
+ Enabling this flag tells the reconstructor to read in color values with the input points and extrapolate those to the vertices of the output.
- [--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.
-The default value for this parameter is 5.
+ [--data <pull factor >]
+ If --colors is specified, this floating point value specifies the relative importance
+of finer color estimates over lower ones.
+The default value for this parameter is 32.
- [--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 .
-The default value for this parameter is the value of the --depth parameter.
+ [--confidence <normal confidence exponent >]
+ This floating point value specifies the exponent to be applied to a point's confidence to adjust its weight. (A point's confidence is defined by the magnitude of its normal.)
+The default value for this parameter is 0.
- [--primalVoxel ]
- Enabling this flag when outputing to a voxel file has the reconstructor sample the implicit function at the corners of the grid, rather than the centers of the cells.
+ [--confidenceBias <normal confidence bias exponent >]
+ This floating point value specifies the exponent to be applied to a point's confidence to bias the resolution at which the sample contributes to the linear system. (Points with lower confidence are biased to contribute at coarser resolutions.)
+The default value for this parameter is 0.
- [--color <pull factor >]
- If specified, the reconstruction code assumes that the input is equipped with colors and will extrapolate
-the color values to the vertices of the reconstructed mesh. The floating point value specifies the relative importance
-of finer color estimates over lower ones. (In practice, we have found that a pull factor of 16 works well.)
+ [--primalGrid ]
+ Enabling this flag when outputing to a grid file has the reconstructor sample the implicit function at the corners of the grid, rather than the centers of the cells.
- [--density ]
- Enabling this flag tells the reconstructor to output the estimated depth values of the iso-surface vertices.
+ [--nonLinearFit ]
+ Enabling this flag has the reconstructor use quadratic interpolation to estimate the positions of iso-vertices.
- [--nonLinearFit ]
- Enabling this flag has the reconstructor use quadratic interpolation to estimate the positions of iso-vertices.
+ [--polygonMesh ]
+ Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes).
- [--polygonMesh ]
- Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes).
+ [--tempDir <temporary output directory >]
+ This string is the name of the directory to which temporary files will be written.
- [--threads <number of processing threads >]
- This integer specifies the number of threads across which the reconstruction
-algorithm should be parallelized.
+ [--threads <number of processing threads >]
+ This integer specifies the number of threads across which the algorithm should be parallelized.
The default value for this parameter is equal to the numer of (virtual) processors on the executing machine.
- [--verbose ]
- Enabling this flag provides a more verbose description of the running times and memory usages of
+ [--maxMemory <maximum memory usage (in GB) >]
+ If positive, this integer value specifies the peak memory utilization for running the reconstruction code (forcing the execution to terminate if the limit is exceeded).
+
+ [--performance ]
+ Enabling this flag provides running time and peak memory usage at the end of the execution.
+
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of
individual components of the surface reconstructor.
-
-
-
+
+
+
+
+
+
+
+
+
+PointInterpolant :
+Fits a function to a set of sample values (and gradients) by finding the coefficients of the function that minimize an energy composed of an interpolation term and Laplacian and bi-Laplacian smoothness terms
+
+
+--in <input sample positions and values >
+ This string is the name of the file from which the positions and values will be read.
+The file should be an ascii file with groups of Dim +1 or 2*Dim +1 (depending on whether gradients are provided or note)
+white space delimited, numbers: the coordinates of the point's position, followed by the value at that point (and gradient).
+No information about the number of samples should be specified.
+
+[--dim <dimension of the samples >]
+ This integerl value is the dimension of the samples.
+The default value is 2.
+
+[--tree <output octree and coefficients >]
+ This string is the name of the file to which the the octree and function coefficients are to be written.
+
+[--grid <output grid >]
+ This string is the name of the file to which the sampled implicit function will be written.
+The file 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 ... bytes corresponding to the (single precision) floating point values
+of the implicit function.
+
+[--degree <B-spline degree >]
+ This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
+Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
+The default value for this parameter is 2.
+
+ [--bType <boundary type >]
+ This integer specifies the boundary type for the finite elements. Valid values are:
+
+ 1 : Free boundary constraints
+ 2 : Dirichlet boundary constraints
+ 3 : Neumann boundary constraints
+
+The default value for this parameter is 1 (free).
+
+ [--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 grid whose resolution is no larger than
+2^d x 2^d x ... Note that since the reconstructor adapts the octree to the
+sampling density, the specified reconstruction depth is only an upper bound.
+The default value for this parameter is 8.
-
-
-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.
+ [--width <finest cell width >]
+
This floating point value specifies the target width of the finest level octree cells.
+This parameter is ignored if the --depth is also specified.
+
+
[--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.
+The default value is 1.1.
+
+
[--valueWeight <value interpolation weight >]
+
This floating point value specifies the importance that interpolation of the samples' values
+is given in the fitting of the function.
+The default value for this parameter is 1000.
+
+
[--gradientWeight <gradient interpolation weight >]
+
This floating point value specifies the importance that interpolation of the samples' gradients
+is given in the fitting of the function.
+The default value for this parameter is 1.
+This value is ignored unless gradient interpolation is specified.
+
+
[--lapWeight <Laplacian weight >]
+
This floating point value specifies the importance that Laplacian regularization
+is given in the fitting of the function.
+The default value for this parameter is 0.
+
+
[--biLapWeight <bi-Laplacian weight >]
+
This floating point value specifies the importance that bi-Laplacian regularization
+is given in the fitting of the function.
+The default value for this parameter is 1.
+
+
[--iters <GS iters >]
+
This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.
+The default value for this parameter is 8.
+
+
[--useGradients ]
+
Enabling this flag indicates that the input file contains gradients as well as sample values.
+
+
[--performance ]
+
Enabling this flag provides running time and peak memory usage at the end of the execution.
+
+
[--verbose ]
+
Enabling this flag provides a more verbose description of the running times and memory usages of
+individual components of the surface reconstructor.
+
+
+
+
+
+
+
+
+
+SurfaceTrimmer :
+Trims off parts of a triangle mesh with a per-vertex signal whose value falls below a threshold (used for removing parts of a reconstructed surface that are generated in low-sampling-density regions)
+
+--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.
The default value is 5.
- [--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.
+ [--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.
The default value 0.001.
- [--polygonMesh ]
- Enabling this flag tells the trimmer to output a polygon mesh (rather than triangulating the trimming results).
-
-
-
+[--polygonMesh ]
+ Enabling this flag tells the trimmer to output a polygon mesh (rather than triangulating the trimming results).
+
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of individual components of the surface reconstructor.
+
+
+
+
+
+
+
+
+
+
+
+ImageStitching :
+Stitches together a composite of image tiles into a seamless panorama by solving for the correction term (solving a 2D Laplacian system) [Agarwala, 2007]
+
+--in <input composite image > <input label image >
+ This pair of strings give the name of the composite image file and the associated label file.
+All pixels in the composite that come from the same source should be assigned the same color in the label file.
+PNG and JPG files are supported (though only PNG should be used for the label file as it is lossless).
+
+ [--out <output image >]
+ This string is the name of the file to which the stitched image will be written.
+PNG and JPG files are supported.
+
+ [--degree <B-spline degree >]
+ This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
+Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
+The default value for this parameter is 1.
+
+ [--wScl <successive under-relaxation scale >]
+ This floating point value specifies the scale for the adapted successive under-relaxation used to remove ringing.
+The default value 0.125.
+
+ [--wExp <successive under-relaxation exponent >]
+ This floating point value specifies the exponent for the adapted successive under-relaxation used to remove ringing.
+The default value 6.
+
+ [--iters <GS iters >]
+ This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.
+The default value for this parameter is 8.
+
+ [--threads <number of processing threads >]
+ This integer specifies the number of threads across which the algorithm should be parallelized.
+The default value for this parameter is equal to the numer of (virtual) processors on the executing machine.
-
-USAGE
+ [--maxMemory <maximum memory usage (in GB) >]
+ If positive, this integer value specifies the peak memory utilization for running the code (forcing the execution to terminate if the limit is exceeded).
+
+ [--performance ]
+ Enabling this flag provides running time and peak memory usage at the end of the execution.
+
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of
+individual components of the image stitcher.
+
+
+
+
+
+
+
+
+
+
+
+EDTInHeat :
+Computes the unsigned Euclidean Distance Transform of a triangle mesh (solving two 3D Laplacian systems) [Crane, Weischedel, and Wardetzky, 2013]
+
+--in <input mesh >
+ This string is the name of the file from which the triangle mesh will be read.
+The file is assumed to be in PLY format.
+
+ [--out <output octree and coefficients >]
+ This string is the name of the file to which the the octree and solution coefficients are to be written.
+
+ [--degree <B-spline degree >]
+ This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
+Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
+The default value for this parameter is 1.
+
+ [--depth <edt depth >]
+ This integer is the maximum depth of the tree that will be used for computing the Euclidean Distance Transform.
+Running at depth d corresponds to solving on a grid whose resolution is no larger than
+2^d x 2^d x ...
+The default value for this parameter is 8.
+
+ [--scale <scale factor >]
+ This floating point value specifies the ratio between the diameter of the cube used for computing the EDT
+and the diameter of the mesh's bounding cube.
+The default value is 2.
+
+ [--diffusion <diffusion time >]
+ This floating point value specifies the time-scale for the initial heat diffusion.
+The default value is 0.0005.
+
+ [--valueWeight <zero-crossing interpolation weight >]
+ This floating point value specifies the importance that the EDT evaluate to zero at points on the input mesh is given.
+The default value for this parameter is 0.01.
+
+ [--wScl <successive under-relaxation scale >]
+ This floating point value specifies the scale for the adapted successive under-relaxation used to remove ringing.
+The default value 0.125.
+
+ [--wExp <successive under-relaxation exponent >]
+ This floating point value specifies the exponent for the adapted successive under-relaxation used to remove ringing.
+The default value 6.
+
+ [--iters <GS iters >]
+ This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.
+The default value for this parameter is 8.
+
+ [--threads <number of processing threads >]
+ This integer specifies the number of threads across which the algorithm should be parallelized.
+The default value for this parameter is equal to the numer of (virtual) processors on the executing machine.
+
+ [--maxMemory <maximum memory usage (in GB) >]
+ If positive, this integer value specifies the peak memory utilization for running the code (forcing the execution to terminate if the limit is exceeded).
+
+ [--performance ]
+ Enabling this flag provides running time and peak memory usage at the end of the execution.
+
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of individual components of the EDT computation.
+
+
+
+
+
+
+
+
+
+
+
+AdaptiveTreeVisualization :
+Extracts iso-surfaces and a sampling on a regular grid from an implicit function represented over an adapted tree
+
+--in <input tree and coefficients >
+ This string is the name of the file from which the tree and implicit functions coefficients are to be read.
+
+[--samples <input sample positions >]
+ This string is the name of the file from which sampling positions are to be read.
+The file should be an ascii file with groups of Dim white space delimited, numbers giving the coordinates of the sampling points' position.
+No information about the number of samples should be specified.
+
+
+[--grid <output value grid >]
+ This string is the name of the file to which the sampling of the implicit along a regular grid will be written.
+The file is written out in binary, with the first 4 bytes corresponding to the (integer) sampling resolution, R ,
+and the next 4 x R ^D bytes corresponding to the (single precision) floating point values of the implicit function. (Here, D is the dimension.)
+
+ [--primalGrid ]
+ Enabling this flag when outputing a grid file samples the implicit function at the corners of the grid, rather than the centers of the cells.
+
+
+ [--mesh <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.
+This is only supported for dimension 3.
+
+ [--iso <iso-value for mesh extraction >]
+ This floating point value specifies the iso-value at which the implicit surface is to be extracted.
+The default value is 0.
+
+ [--nonLinearFit ]
+ Enabling this flag has the reconstructor use quadratic interpolation to estimate the positions of iso-vertices.
+
+ [--polygonMesh ]
+ Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes).
+
+ [--flip ]
+ Enabling this flag flips the orientations of the output triangles.
+
+ [--threads <number of processing threads >]
+ This integer specifies the number of threads across which the algorithm should be parallelized.
+The default value for this parameter is equal to the numer of (virtual) processors on the executing machine.
+
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of
+individual components of the visualizer.
+
+
+
+
+
+
+
+
+
+
+ChunkPly :
+Decomposes a single mesh/point-set file into a set of chunks with prescribed bounding box widths.
+
+--in <input ply file >
+ This string is the name of the file containing the geometry that is to be chunked.
+
+ [--out <output ply file name/header >]
+ This string is the name of the file/header to which the chunks should be written. If the width of the chunk is W , the file containing the geometry inside the cube [W ·i ,W ·(i+1 )]×[W ·j ,W ·(j+1 )]×[W ·k ,W ·(k+1 )] will be named <output header>.i.j.k.ply .
+
+ [--width <chunk width > ]
+ This floating point value specifies the width of the cubes used for chunking.
+The default value for this parameter is -1 , indicating that the input should be written to a single ouput. (In this case the value of the --out parameter is the name of the single file to which the output is written.
+
+ [--radius <padding radius > ]
+ This floating point value specifies the size of the padding region used, as a fraction of the total width of the cube.
+The default value for this parameter is 0 , indicating that no padding should be used.
+
+ [--verbose ]
+ Enabling this flag provides a more verbose description of the running times and memory usages of
+individual components of the visualizer.
+
+
+
+
+
+
+
+USAGE EXAMPLES (WITH SAMPLE DATA)
+
+
+
+
+
+PoissonRecon / SSDRecon / SurfaceTrimmer / ChunkPly
+
For testing purposes, three point sets are provided:
-
+
- Eagle :
-A set of 796,825 oriented point samples with color (represented in PLY format) was obtained in the EPFL Scanning 3D Statues from Photos course.
-
-The original Poisson Reconstruction algorithm can be invoked by calling:
-% PoissonRecon --in eagle.points.ply --out eagle.unscreened.ply --depth 10 --pointWeight 0
-using the --pointWeight 0 argument to disable the screening.
-
- By default, screening is enabled so the call:
-% PoissonRecon --in eagle.points.ply --out eagle.screened.ply --depth 10
-produces a reconstruction that more faithfully fits the input point positions.
-
- A reconstruction of the eagle that extrapolates the color values from the input samples can be obtained by calling:
-% PoissonRecon --in eagle.points.ply --out eagle.screened.color.ply --depth 10 --color 16
-using the --color 16 to indicate both that color should be used, and the extent to which finer color estimates should be preferenced over coarser estimates.
-
- Finally, a reconstruction the eagle that does not close up the holes can be obtained by first calling:
-% PoissonRecon --in eagle.points.ply --out eagle.screened.color.ply --depth 10 --color 16 --density
-using the --density flag to indicate that density estimates should be output with the vertices of the mesh, and then calling:
-% SurfaceTrimmer --in eagle.screened.color.ply --out eagle.screened.color.trimmed.ply --trim 7
-to remove all subsets of the surface where the sampling density corresponds to a depth smaller than 7.
-
-
-
-
-
-
-
-
-
-Unscreened
- Screened
- Screened + Color
- Screened + Color + Trimmed
-
-
-
- Bunny :
+ Horse :
+A set of 100,000 oriented point samples (represented in ASCII format) was obtained by sampling a virtual horse model with a sampling density proportional to curvature, giving a set of non-uniformly distributed points.
+The surface of the model can be reconstructed by calling the either Poisson surface reconstructor:
+% PoissonRecon --in horse.npts --out horse.ply --depth 10
+or the SSD surface reconstructor
+% SSDRecon --in horse.npts --out horse.ply --depth 10
+
+
+ Bunny :
A set of 362,271 oriented point samples (represented in PLY format) was obtained by merging the data from the original Stanford Bunny
-range scans . The orientation of the sample points was estimated
+range scans . The orientation of the sample points was estimated
using the connectivity information within individual range scans.
-The surface of the model can be reconstructed by calling the surface reconstructor as follows:
-% PoissonRecon --in bunny.points.ply --out bunny.ply --depth 10
-
- Horse :
-A set of 100,000 oriented point samples (represented in ASCII format) was obtained by sampling a virtual horse model with a sampling density proportional to curvature, giving a set of non-uniformly distributed points.
-The surface of the model can be reconstructed by calling the surface reconstructor as follows:
-% 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.
-As an examples, the reconstructed bunny can be converted into the ASCII mesh format as follows:
-% ply2mesh.pl bunny.ply > bunny.m
-
-
-CHANGES
-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 ]
-
-
-Version 7.0 :
-
- Added functionality to support color extrapolation if present in the input.
- Modified a bug with the way in which sample contributions were scaled.
-
-
-Version 8.0 :
-
- Added support for different degree B-splines.
+The original (unscreened) Poisson reconstruction can be obtained by setting the point interpolation weight to zero:
+% PoissonRecon --in bunny.points.ply --out bunny.ply --depth 10 --pointWeight 0
+By default, the Poisson surface reconstructor uses degree-2 B-splines. A more efficient reconstruction can be obtained using degree-1 B-splines:
+% PoissonRecon --in bunny.points.ply --out bunny.ply --depth 10 --pointWeight 0 --degree 1
+(The SSD reconstructor requires B-splines of degree at least 2 since second derivatives are required to formulate the bi-Laplacian energy.)
+
+
+ Eagle :
+A set of 796,825 oriented point samples with color (represented in PLY format) was obtained in the EPFL Scanning 3D Statues from Photos course.
+A reconstruction of the eagle that extrapolates the color values from the input samples can be obtained by calling:
+% PoissonRecon --in eagle.points.ply --out eagle.pr.color.ply --depth 10 --colors
+using the --colors flag to indicate that color extrapolation should be used.
+A reconstruction of the eagle that does not close up the holes can be obtained by first calling:
+% SSDRecon --in eagle.points.ply --out eagle.ssd.color.ply --depth 10 --colors --density
+using the --density flag to indicate that density estimates should be output with the vertices of the mesh, and then calling:
+% SurfaceTrimmer --in eagle.ssd.color.ply --out eagle.ssd.color.trimmed.ply --trim 7
+to remove all subsets of the surface where the sampling density corresponds to a depth smaller than 7.
+This reconstruction can be chunked into cubes of size 4×4×4 by calling:
+% ChunkPly --in eagle.ssd.color.trimmed.ply --out eagle.ssd.color.trimmed.chnks --width 4
+which partitions the reconstruction into 11 pieces.
+
+
+
+
+
+
+
+
+
+
+
+
+
+PointInterpolant / AdaptiveTreeVisualization
+
+For testing purposes, a pair of point-sets is provided:
+
+
+ fitting samples :
+A set of 1000 random 2D samples from within the square [-1,1,]x[-1,1] along with the evaluation of the quadratic f(x,y)=x*x+y*y at each sample point (represented in ASCII format).
+ evaluation samples :
+A set of 4 2D positions at which the fit function is to be evaluated (represented in ASCII format).
+
+
+The function fitting the input samples can be by calling the point interpolant:
+% PointInterpolant --in quadratic.2D.fitting.samples --out quadratic.2D.tree --dim 2
+Then, the reconstructed function can be evaluated at the evaluation samples by calling the adaptive tree visualization:
+% AdaptiveTreeVisualization --in quadratic.2D.tree --samples quadratic.2D.evaluation.samples
+This will output the evaluation positions and values:
+0 0 1.33836e-05
+0.5 0 0.25001
+0.5 0.5 0.500006
+2 2 nan
+Note that because the last evaluation position, (2,2), is outside the bounding box of the fitting samples, the function cannot be evaluated at this point and a value of "nan" is output.
+
+
+
+
+
+
+
+
+ImageStitching
+
+For testing purposes, two panoramas are provided: Jaffa (23794 x 9492 pixels) and OldRag (87722 x 12501 pixels).
+
+A seamless panorama can be obtained by running:
+% ImageSitching --in pixels.png labels.png --out out.png
+
+
+
+
+
+
+
+
+
+
+EDTInHeat / AdaptiveTreeVisualization
+
+The Euclidean Distance Tranform of the reconstructed horse can be obtained by running:
+% EDTInHeat --in horse.ply --out horse.edt --depth 9
+Then, the visualization code can be used to extract iso-surfaces from the implicit function.
+To obtain a visualization near the input surface, use an iso-value close to zero:
+% AdaptiveTreeVisualization.exe --in horse.edt --mesh horse_0.01_.ply --iso 0.01 --flip
+(By default, the surface is aligned so that the outward facing normal aligns with the negative gradient. Hence, specifying the --flip flag is used to re-orient the surface.)
+To obtain a visualization closer to the boundary of the bounding box, use an iso-value close to zero:
+% AdaptiveTreeVisualization.exe --in horse.edt --mesh horse_0.25_.ply --iso 0.25 --flip
+(Since the default --scale is 2, a value of 0.25 should still give a surface that is contained within the bounding box.)
+To obtain a sampling of the implicit function over a regular grid:
+% AdaptiveTreeVisualization.exe --in horse.edt --grid horse.grid
+
+
+
+
+
+
+
+
+
+HISTORY OF CHANGES
+
+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 ]
+
+
+Version 7.0 :
+
+ Added functionality to support color extrapolation if present in the input.
+ Modified a bug with the way in which sample contributions were scaled.
+
+
+Version 8.0 :
+
+ Added support for different degree B-splines.
(Note that as the B-spline degree is a template parameter, only degree 1 through 4 are supported.
-If higher order degrees are desired, additional template parameters can be easily added in the body of the Execute function inside of PoissonRecon.cpp .
+If higher order degrees are desired, additional template parameters can be easily added in the body of the Execute function inside of PoissonRecon.cpp .
Similarly, to reduce compilation times, support for specific degrees can be removed.)
- Added the --primalVoxel flag to support to extraction of a voxel grid using primal sampling.
- Changed the implementation of the voxel sampling so that computation is now linear, rather than log-linear, in the number of samples.
-
-
-Version 9.0 :
-
- Added support for free boundary conditions.
- Extended the solver to support more general linear systems. This makes it possible to use the same framework to implement the Smoothed Signed Distance Reconstruction of Calakli and Taubin (2011).
- Modified the implementation of density estimation and input representation. This tends to define a slightly larger system. On its own, this results in slightly increased running-time/footprint for full-res reconstructions, but provides a substantially faster implementation when the output complexity is smaller than the input.
-
-
-Version 9.01 :
-
- Reverted the density estimation to behave as in Version 8.0.
-
-
-
-SUPPORT
+ Added the --primalGrid flag to support to extraction of a grid using primal sampling.
+ Changed the implementation of the grid sampling so that computation is now linear, rather than log-linear, in the number of samples.
+
+
+Version 9.0 :
+
+ Added support for free boundary conditions.
+ Extended the solver to support more general linear systems. This makes it possible to use the same framework to implement the Smoothed Signed Distance Reconstruction of Calakli and Taubin (2011).
+ Modified the implementation of density estimation and input representation. This tends to define a slightly larger system. On its own, this results in slightly increased running-time/footprint for full-res reconstructions, but provides a substantially faster implementation when the output complexity is smaller than the input.
+
+
+Version 9.01 :
+
+ Reverted the density estimation to behave as in Version 8.0.
+
+
+Version 9.011 :
+
+ Added a parameter for specifying the temporary directory.
+
+
+Version 10.00 :
+
+ The code has been reworked to support arbitrary dimensions, finite elements of arbitrary degree, generally SPD systems in the evaluated/integrated values and derivatives of the functions, etc.
+ For the reconstruction code, added the --width flag which allows the system to compute the depth of the octree given a target depth for the finest resolution nodes.
+ For the reconstruction code, fixed a bug in the handling of the confidence encoded in the lengths of the normals. In addition, added the flags --confidence and --confidenceBias which allow the user more control of how confidence is used to affect the contribution of a sample.
+
+
+Version 10.01 :
+
+ Modified the reconstruction code to facilitate interpolation of other input-sample quantities, in addition to color.
+
+
+Version 10.02 :
+
+ Set the default value for --degree in PoissonRecon to 1 and change the definitiion of DATA_DEGREE to 0 for sharper color interpolation.
+
+
+Version 10.03 :
+
+ Cleaned up memory leaks and fixed a bug causing ImageStitching and EDTInHeat to SEGFAULT on Linux.
+
+
+Version 10.04 :
+
+ Replaced the ply I/O code with an object-oriented implementation.
+ Updated the code to support compilation under gcc version 4.8.
+
+
+Version 10.05 :
+
+ Added cleaner support for warning and error handling.
+ Minor bug fixes.
+ Added a --inCore flag that enables keeping the pointset in memory instead of streaming it in from disk.
+
+
+Version 10.06 :
+
+ Improved performance.
+ Modified PoissonRecon and SSDRecon to support processing of 2D point sets.
+ Modified the 2D implementations of PoissonRecon, SSDRecon, and AdaptiveTreeVisualization to support ouput to .jpg and .png image files.
+
+
+Version 10.07 :
+
+ Removed a bug that would cause memory access errors when some slices were empty.
+
+
+Version 11.00 :
+
+ Added support for processing point-sets so large that 32-bit indices for octrees are not sufficient. (Enabled by defining the preprocessor variable BIG_DATA in the file PreProcessor.h .
+ Added C++11 parallelism for compilers that do not support OpenMP.
+ Added the code for ChunkPly which breaks up large meshes and/or point-sets into chunks.
+
+
+Version 11.01 :
+
+ Fixed bug with _mktemp that caused the code to crash on Windows machine with more than 26 cores.
+
+
+Version 11.02 :
+
+ Added error handling for numerical imprecision issues arrising when too many samples fall into a leaf node.
+
+
+Version 12.00 :
+
+ Added functionality enabling AdaptiveTreeVisualization to output the values of a function at prescribed sample positions.
+ Added the implementation of PointInterpolant that fits a function to a discrete set of sample values.
+
+
+
+
+
+
+SUPPORT
This work genersouly supported by NSF grants #0746039 and #1422325.
-
-
-HOME
diff --git a/SSDRecon.vcxproj b/SSDRecon.vcxproj
deleted file mode 100644
index 5ddf8894..00000000
--- a/SSDRecon.vcxproj
+++ /dev/null
@@ -1,235 +0,0 @@
-
-
-
-
- Debug
- Win32
-
-
- Debug
- x64
-
-
- Release
- Win32
-
-
- Release
- x64
-
-
-
- SSDRecon
- {7838CA1E-8A39-4A2B-AC3D-3E25FEAEA2D6}
- PoissonRecon
- Win32Proj
- 8.1
-
-
-
- Application
- MultiByte
- true
- v140
-
-
- Application
- MultiByte
- v140
-
-
- Application
- MultiByte
- true
- v140
-
-
- Application
- MultiByte
- v140
- true
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- <_ProjectFileVersion>10.0.30319.1
- $(SolutionDir)\Bin\$(Platform)\$(Configuration)\
- $(SolutionDir)\Obj\$(Platform)\$(Configuration)\
- true
- $(SolutionDir)\Bin\$(Platform)\$(Configuration)\
- $(SolutionDir)\Obj\$(TargetName)\$(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
-
-
- 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/SSDRecon.vcxproj.filters b/SSDRecon.vcxproj.filters
deleted file mode 100644
index ce331d6a..00000000
--- a/SSDRecon.vcxproj.filters
+++ /dev/null
@@ -1,143 +0,0 @@
-
-
-
-
- {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
-
-
-
-
- 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
-
-
- Include Files
-
-
- Include Files
-
-
-
\ No newline at end of file
diff --git a/Src/Array.inl b/Src/Array.inl
deleted file mode 100644
index 6247f528..00000000
--- a/Src/Array.inl
+++ /dev/null
@@ -1,662 +0,0 @@
-/*
-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.data;
- 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; }
- inline Array operator++( int ){ Array< C > temp = (*this) ; (*this) +=1 ; return temp; }
- 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; }
- inline Array operator--( int ){ Array< C > temp = (*this) ; (*this) -=1 ; return temp; }
- 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; }
- inline ConstArray operator++( int ){ ConstArray< C > temp = (*this) ; (*this) +=1 ; return temp; }
- 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; }
- inline ConstArray operator--( int ){ ConstArray< C > temp = (*this) ; (*this) -=1 ; return temp; }
- 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() 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 CornerAndWidth( int depth , int offset , Real& corner , Real& width ){ width = Real(1.0/(1< static inline void CenterAndWidth( int idx , Real& center , Real& width )
- {
- int depth , offset;
- CenterDepthAndOffset( idx , depth , offset );
- CenterAndWidth( depth , offset , center , width );
- }
- template< class Real > static inline void CornerAndWidth( int idx , Real& corner , Real& width )
- {
- int depth , offset;
- CornerDepthAndOffset( idx , depth , offset );
- CornerAndWidth( depth , offset , corner , width );
- }
- static inline void CenterDepthAndOffset( int idx , int& depth , int& offset )
- {
- offset = idx , depth = 0;
- while( offset>=(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
deleted file mode 100644
index eeded680..00000000
--- a/Src/CmdLineParser.inl
+++ /dev/null
@@ -1,141 +0,0 @@
-/* -*- 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/Geometry.cpp b/Src/Geometry.cpp
deleted file mode 100644
index e6d03d88..00000000
--- a/Src/Geometry.cpp
+++ /dev/null
@@ -1,123 +0,0 @@
-/*
-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( const char* fileName , const char* fileHeader , int bufferSize )
-{
- _bufferIndex = 0;
- _bufferSize = bufferSize;
- if( fileName ) strcpy( _fileName , fileName ) , tempFile = false , _fp = fopen( _fileName , "w+b" );
- else
- {
- if( fileHeader && strlen(fileHeader) ) sprintf( _fileName , "%sXXXXXX" , fileHeader );
- else strcpy( _fileName , "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
deleted file mode 100644
index ea8c951c..00000000
--- a/Src/Geometry.h
+++ /dev/null
@@ -1,414 +0,0 @@
-/*
-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
-#include
-
-template
-Real Random(void);
-
-template< class Real >
-struct Point3D
-{
- Real coords[3];
- Point3D( void ) { coords[0] = coords[1] = coords[2] = Real(0); }
- Point3D( Real v ) { coords[0] = coords[1] = coords[2] = v; }
- 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< class Real >
-struct OrientedPoint3D
-{
- Point3D< Real > p , n;
- OrientedPoint3D( Point3D< Real > pp=Point3D< Real >() , Point3D< Real > nn=Point3D< Real >() ) : p(pp) , n(nn) { ; }
- template< class _Real > OrientedPoint3D( const OrientedPoint3D< _Real >& p ) : OrientedPoint3D( Point3D< Real >( p.p ) , Point3D< Real >( p.n ) ){ ; }
-
- template< class _Real > inline OrientedPoint3D& operator += ( OrientedPoint3D< _Real > _p ){ p += _p.p , n += _p.n ; return *this; }
- template< class _Real > inline OrientedPoint3D operator + ( OrientedPoint3D< _Real > _p ) const { return OrientedPoint3D< Real >( p+_p.p , n+_p.n ); }
- template< class _Real > inline OrientedPoint3D& operator *= ( _Real r ) { p *= r , n *= r ; return *this; }
- template< class _Real > inline OrientedPoint3D operator * ( _Real r ) const { return OrientedPoint3D< Real >( p*r , n*r ); }
-
- template< class _Real > inline OrientedPoint3D& operator -= ( OrientedPoint3D< _Real > p ){ return ( (*this)+=(-p) ); }
- template< class _Real > inline OrientedPoint3D operator - ( OrientedPoint3D< _Real > p ) const { return (*this)+(-p); }
- template< class _Real > inline OrientedPoint3D& operator /= ( _Real r ){ return ( (*this)*=Real(1./r) ); }
- template< class _Real > inline OrientedPoint3D operator / ( _Real r ) const { return (*this) * ( Real(1.)/r ); }
-};
-
-template< class Data , class Real >
-struct ProjectiveData
-{
- Data data;
- Real weight;
- ProjectiveData( Data d=Data(0) , Real w=Real(0) ) : data(d) , weight(w) { ; }
- operator Data (){ return weight!=0 ? data/weight : data*weight; }
- ProjectiveData& operator += ( const ProjectiveData& p ){ data += p.data , weight += p.weight ; return *this; }
- ProjectiveData& operator -= ( const ProjectiveData& p ){ data -= p.data , weight -= p.weight ; return *this; }
- ProjectiveData& operator *= ( Real s ){ data *= s , weight *= s ; return *this; }
- ProjectiveData& operator /= ( Real s ){ data /= s , weight /= s ; return *this; }
- ProjectiveData operator + ( const ProjectiveData& p ) const { return ProjectiveData( data+p.data , weight+p.weight ); }
- ProjectiveData operator - ( const ProjectiveData& p ) const { return ProjectiveData( data-p.data , weight-p.weight ); }
- ProjectiveData operator * ( Real s ) const { return ProjectiveData( data*s , weight*s ); }
- ProjectiveData operator / ( Real s ) const { return ProjectiveData( data/s , weight/s ); }
-};
-
-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:
- std::unordered_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( const char* fileName=NULL , const char* fileHeader="" , 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( const char* fileHeader="" );
- ~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
deleted file mode 100644
index ad2be93b..00000000
--- a/Src/Geometry.inl
+++ /dev/null
@@ -1,591 +0,0 @@
-/*
-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)
-{
- std::unordered_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
-int CoredVectorMeshData< Vertex >::nextPolygon( std::vector< CoredVertexIndex >& vertices )
-{
- if( polygonIndex& polygon = polygons[ polygonIndex++ ];
- vertices.resize( polygon.size() );
- for( int i=0 ; i
-int CoredVectorMeshData< Vertex >::outOfCorePointCount(void){return int(oocPoints.size());}
-template< class Vertex >
-int CoredVectorMeshData< Vertex >::polygonCount( void ) { return int( polygons.size() ); }
-
-///////////////////////
-// CoredFileMeshData //
-///////////////////////
-template< class Vertex >
-CoredFileMeshData< Vertex >::CoredFileMeshData( const char* fileHeader )
-{
- oocPoints = polygons = 0;
-
- oocPointFile = new BufferedReadWriteFile( NULL , fileHeader );
- polygonFile = new BufferedReadWriteFile( NULL , fileHeader );
-}
-template< class Vertex >
-CoredFileMeshData< Vertex >::~CoredFileMeshData( void )
-{
- delete oocPointFile;
- delete polygonFile;
-}
-template< class Vertex >
-void CoredFileMeshData< Vertex >::resetIterator ( void )
-{
- oocPointFile->reset();
- polygonFile->reset();
-}
-template< class Vertex >
-int CoredFileMeshData< Vertex >::addOutOfCorePoint( const Vertex& p )
-{
- oocPointFile->write( &p , sizeof( Vertex ) );
- oocPoints++;
- return oocPoints-1;
-}
-template< class Vertex >
-int CoredFileMeshData< Vertex >::addOutOfCorePoint_s( const Vertex& p )
-{
- int sz;
-#pragma omp critical (CoredFileMeshData_addOutOfCorePoint_s)
- {
- sz = oocPoints;
- oocPointFile->write( &p , sizeof( Vertex ) );
- oocPoints++;
- }
- return sz;
-}
-template< class Vertex >
-int CoredFileMeshData< Vertex >::addPolygon_s( const std::vector< int >& vertices )
-{
- int sz , vSize = (int)vertices.size();
-#pragma omp critical (CoredFileMeshData_addPolygon_s )
- {
- sz = polygons;
- polygonFile->write( &vSize , sizeof(int) );
- polygonFile->write( &vertices[0] , sizeof(int) * vSize );
- polygons++;
- }
- return sz;
-}
-template< class Vertex >
-int CoredFileMeshData< 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 CoredFileMeshData< Vertex >::nextOutOfCorePoint( Vertex& p )
-{
- if( oocPointFile->read( &p , sizeof( Vertex ) ) ) return 1;
- else return 0;
-}
-template< class Vertex >
-int CoredFileMeshData< Vertex >::nextPolygon( std::vector< CoredVertexIndex >& vertices )
-{
- int pSize;
- if( polygonFile->read( &pSize , sizeof(int) ) )
- {
- std::vector< int > polygon( pSize );
- if( polygonFile->read( &polygon[0] , sizeof(int)*pSize ) )
- {
- vertices.resize( pSize );
- for( int i=0 ; i
-int CoredFileMeshData< Vertex >::outOfCorePointCount( void ){ return oocPoints; }
-template< class Vertex >
-int CoredFileMeshData< Vertex >::polygonCount( void ) { return polygons; }
diff --git a/Src/MAT.inl b/Src/MAT.inl
deleted file mode 100644
index 5106659f..00000000
--- a/Src/MAT.inl
+++ /dev/null
@@ -1,217 +0,0 @@
-/*
-Copyright (c) 2007, Michael Kazhdan
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without modification,
-are permitted provided that the following conditions are met:
-
-Redistributions of source code must retain the above copyright notice, this list of
-conditions and the following disclaimer. Redistributions in binary form must reproduce
-the above copyright notice, this list of conditions and the following disclaimer
-in the documentation and/or other materials provided with the distribution.
-
-Neither the name of the Johns Hopkins University nor the names of its contributors
-may be used to endorse or promote products derived from this software without specific
-prior written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
-EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
-OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
-SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
-INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
-TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
-BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
-ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
-DAMAGE.
-*/
-//////////////////////////////
-// MinimalAreaTriangulation //
-//////////////////////////////
-template
-MinimalAreaTriangulation::MinimalAreaTriangulation(void)
-{
- bestTriangulation=NULL;
- midPoint=NULL;
-}
-template
-MinimalAreaTriangulation::~MinimalAreaTriangulation(void)
-{
- if(bestTriangulation)
- delete[] bestTriangulation;
- bestTriangulation=NULL;
- if(midPoint)
- delete[] midPoint;
- midPoint=NULL;
-}
-template
-void MinimalAreaTriangulation::GetTriangulation(const std::vector >& vertices,std::vector& triangles)
-{
- if(vertices.size()==3)
- {
- triangles.resize(1);
- triangles[0].idx[0]=0;
- triangles[0].idx[1]=1;
- triangles[0].idx[2]=2;
- return;
- }
- else if(vertices.size()==4)
- {
- TriangleIndex tIndex[2][2];
- Real area[2];
-
- area[0]=area[1]=0;
- triangles.resize(2);
-
- tIndex[0][0].idx[0]=0;
- tIndex[0][0].idx[1]=1;
- tIndex[0][0].idx[2]=2;
- tIndex[0][1].idx[0]=2;
- tIndex[0][1].idx[1]=3;
- tIndex[0][1].idx[2]=0;
-
- tIndex[1][0].idx[0]=0;
- tIndex[1][0].idx[1]=1;
- tIndex[1][0].idx[2]=3;
- tIndex[1][1].idx[0]=3;
- tIndex[1][1].idx[1]=1;
- tIndex[1][1].idx[2]=2;
-
- Point3D n,p1,p2;
- for(int i=0;i<2;i++)
- for(int j=0;j<2;j++)
- {
- p1=vertices[tIndex[i][j].idx[1]]-vertices[tIndex[i][j].idx[0]];
- p2=vertices[tIndex[i][j].idx[2]]-vertices[tIndex[i][j].idx[0]];
- CrossProduct(p1,p2,n);
- area[i] += Real( Length(n) );
- }
- if(area[0]>area[1])
- {
- triangles[0]=tIndex[1][0];
- triangles[1]=tIndex[1][1];
- }
- else
- {
- triangles[0]=tIndex[0][0];
- triangles[1]=tIndex[0][1];
- }
- return;
- }
- if(bestTriangulation)
- delete[] bestTriangulation;
- if(midPoint)
- delete[] midPoint;
- bestTriangulation=NULL;
- midPoint=NULL;
- size_t eCount=vertices.size();
- bestTriangulation=new Real[eCount*eCount];
- midPoint=new int[eCount*eCount];
- for(size_t i=0;i
-Real MinimalAreaTriangulation::GetArea(const std::vector >& vertices)
-{
- if(bestTriangulation)
- delete[] bestTriangulation;
- if(midPoint)
- delete[] midPoint;
- bestTriangulation=NULL;
- midPoint=NULL;
- int eCount=vertices.size();
- bestTriangulation=new double[eCount*eCount];
- midPoint=new int[eCount*eCount];
- for(int i=0;i
-void MinimalAreaTriangulation::GetTriangulation(const size_t& i,const size_t& j,const std::vector >& vertices,std::vector& triangles)
-{
- TriangleIndex tIndex;
- size_t eCount=vertices.size();
-#ifdef BRUNO_LEVY_FIX
- int ii=(int)i;
- if( i=ii )
- return;
- ii=midPoint[i*eCount+j];
- if( ii>=0 )
- {
- tIndex.idx[0] = int( i );
- tIndex.idx[1] = int( j );
- tIndex.idx[2] = int( ii );
- triangles.push_back(tIndex);
- GetTriangulation(i,ii,vertices,triangles);
- GetTriangulation(ii,j,vertices,triangles);
- }
-}
-
-template
-Real MinimalAreaTriangulation::GetArea(const size_t& i,const size_t& j,const std::vector >& vertices)
-{
- Real a=FLT_MAX,temp;
- size_t eCount=vertices.size();
- size_t idx=i*eCount+j;
- size_t ii=i;
- if(i=ii)
- {
- bestTriangulation[idx]=0;
- return 0;
- }
- if(midPoint[idx]!=-1)
- return bestTriangulation[idx];
- int mid=-1;
- for(size_t r=j+1;r p,p1,p2;
- p1=vertices[i]-vertices[rr];
- p2=vertices[j]-vertices[rr];
- CrossProduct(p1,p2,p);
- temp = Real( Length(p) );
- if(bestTriangulation[idx1]>=0)
- {
- temp+=bestTriangulation[idx1];
- if(temp>a)
- continue;
- if(bestTriangulation[idx2]>0)
- temp+=bestTriangulation[idx2];
- else
- temp+=GetArea(rr,j,vertices);
- }
- else
- {
- if(bestTriangulation[idx2]>=0)
- temp+=bestTriangulation[idx2];
- else
- temp+=GetArea(rr,j,vertices);
- if(temp>a)
- continue;
- temp+=GetArea(i,rr,vertices);
- }
-
- if(temp
-#include "MarchingCubes.h"
-
-////////////
-// Square //
-////////////
-int Square::AntipodalCornerIndex(int idx){
- int x,y;
- FactorCornerIndex(idx,x,y);
- return CornerIndex( (x+1)%2 , (y+1)%2 );
-}
-int Square::CornerIndex( int x , int y ){ return (y<<1)|x; }
-void Square::FactorCornerIndex( int idx , int& x , int& y ){ x=(idx>>0)&1 , y=(idx>>1)&1; }
-int Square::EdgeIndex( int orientation , int i )
-{
- switch( orientation )
- {
- case 0: // x
- if( !i ) return 0; // (0,0) -> (1,0)
- else return 2; // (0,1) -> (1,1)
- case 1: // y
- if( !i ) return 3; // (0,0) -> (0,1)
- else return 1; // (1,0) -> (1,1)
- };
- return -1;
-}
-void Square::FactorEdgeIndex(int idx,int& orientation,int& i){
- switch(idx){
- case 0: case 2:
- orientation=0;
- i=idx/2;
- return;
- case 1: case 3:
- orientation=1;
- i=((idx/2)+1)%2;
- return;
- };
-}
-void Square::EdgeCorners(int idx,int& c1,int& c2){
- int orientation,i;
- FactorEdgeIndex(idx,orientation,i);
- switch(orientation){
- case 0:
- c1 = CornerIndex(0,i);
- c2 = CornerIndex(1,i);
- break;
- case 1:
- c1 = CornerIndex(i,0);
- c2 = CornerIndex(i,1);
- break;
- };
-}
-int Square::ReflectEdgeIndex(int idx,int edgeIndex){
- int orientation=edgeIndex%2;
- int o,i;
- FactorEdgeIndex(idx,o,i);
- if(o!=orientation){return idx;}
- else{return EdgeIndex(o,(i+1)%2);}
-}
-int Square::ReflectCornerIndex(int idx,int edgeIndex){
- int orientation=edgeIndex%2;
- int x,y;
- FactorCornerIndex(idx,x,y);
- switch(orientation){
- case 0: return CornerIndex((x+1)%2,y);
- case 1: return CornerIndex(x,(y+1)%2);
- };
- return -1;
-}
-
-
-
-//////////
-// Cube //
-//////////
-int Cube::CornerIndex( int x , int y , int z ){ return (z<<2)|(y<<1)|x; }
-void Cube::FactorCornerIndex( int idx , int& x , int& y , int& z ){ x = (idx>>0)&1 , y = (idx>>1)&1 , z = (idx>>2)&1; }
-int Cube::EdgeIndex(int orientation,int i,int j){return (i | (j<<1))|(orientation<<2);}
-void Cube::FactorEdgeIndex( int idx , int& orientation , int& i , int &j )
-{
- orientation=idx>>2;
- i = (idx&1);
- j = (idx&2)>>1;
-}
-int Cube::FaceIndex( int x , int y , int z )
-{
- if ( x<0 ) return 0;
- else if( x>0 ) return 1;
- else if( y<0 ) return 2;
- else if( y>0 ) return 3;
- else if( z<0 ) return 4;
- else if( z>0 ) return 5;
- else return -1;
-}
-int Cube::FaceIndex( int dir , int offSet ){ return (dir<<1)|offSet; }
-
-void Cube::FactorFaceIndex( int idx , int& x , int& y , int& z )
-{
- x=y=z=0;
- switch( idx )
- {
- case 0: x=-1; break;
- case 1: x= 1; break;
- case 2: y=-1; break;
- case 3: y= 1; break;
- case 4: z=-1; break;
- case 5: z= 1; break;
- };
-}
-void Cube::FactorFaceIndex( int idx , int& dir , int& offSet )
-{
- dir = idx>>1;
- offSet=idx &1;
-}
-bool Cube::IsEdgeCorner( int cIndex , int e )
-{
- int o , i , j;
- FactorEdgeIndex( e , o , i , j );
- switch( o )
- {
- case 0: return (cIndex && 2)==(i<<1) && (cIndex && 4)==(j<<2);
- case 1: return (cIndex && 1)==(i<<0) && (cIndex && 4)==(j<<2);
- case 2: return (cIndex && 4)==(i<<2) && (cIndex && 2)==(j<<1);
- default: return false;
- }
-}
-bool Cube::IsFaceCorner( int cIndex , int f )
-{
- int dir , off;
- FactorFaceIndex( f , dir , off );
- return ( cIndex & (1< (1,0)
-1} // (1,0) -> (1,1)
-2} // (0,1) -> (1,1)
-3} // (0,0) -> (0,1)
-*/
-const int MarchingSquares::edgeMask[1< -> ->
- 9, // 1 -> 0 -> (0,0) -> 0,3 -> 9
- 3, // 2 -> 1 -> (1,0) -> 0,1 -> 3
- 10, // 3 -> 0,1 -> (0,0) (1,0) -> 1,3 -> 10
- 12, // 4 -> 2 -> (0,1) -> 2,3 -> 12
- 5, // 5 -> 0,2 -> (0,0) (0,1) -> 0,2 -> 5
- 15, // 6 -> 1,2 -> (1,0) (0,1) -> 0,1,2,3 -> 15
- 6, // 7 -> 0,1,2 -> (0,0) (1,0) (0,1) -> 1,2 -> 6
- 6, // 8 -> 3 -> (1,1) -> 1,2 -> 6
- 15, // 9 -> 0,3 -> (0,0) (1,1) -> 0,1,2,3 -> 15
- 5, // 10 -> 1,3 -> (1,0) (1,1) -> 0,2 -> 5
- 12, // 11 -> 0,1,3 -> (0,0) (1,0) (1,1) -> 2,3 -> 12
- 10, // 12 -> 2,3 -> (0,1) (1,1) -> 1,3 -> 10
- 3, // 13 -> 0,2,3 -> (0,0) (0,1) (1,1) -> 0,1 -> 3
- 9, // 14 -> 1,2,3 -> (1,0) (0,1) (1,1) -> 0,3 -> 9
- 0, // 15 -> 0,1,2,3 -> (0,0) (1,0) (0,1) (1,1) ->
-};
-#if NEW_ORDERING
-/*
-0} // (0,0) -> (1,0)
-1} // (1,0) -> (1,1)
-2} // (0,1) -> (1,1)
-3} // (0,0) -> (0,1)
-*/
-const int MarchingSquares::edges[1<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(1,i,j)];}}}
- else if (y<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,0,j)];}}}
- else if (y>0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,1,j)];}}}
- else if (z<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,j,0)];}}}
- else if (z>0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,j,1)];}}}
- if (v[0][0] < iso) idx |= 1;
- if (v[1][0] < iso) idx |= 2;
- if (v[1][1] < iso) idx |= 4;
- if (v[0][1] < iso) idx |= 8;
- return idx;
-}
-bool MarchingCubes::IsAmbiguous( const double v[Cube::CORNERS] , double isoValue , int faceIndex ){ return MarchingSquares::IsAmbiguous( GetFaceIndex( v , isoValue , faceIndex ) ); }
-bool MarchingCubes::HasRoots( const double v[Cube::CORNERS] , double isoValue , int faceIndex ){ return MarchingSquares::HasRoots( GetFaceIndex( v , isoValue , faceIndex ) ); }
-bool MarchingCubes::HasRoots( const double v[Cube::CORNERS] , double isoValue ){ return HasRoots( GetIndex( v , isoValue ) ); }
-bool MarchingCubes::HasRoots( unsigned char mcIndex ){ return !(mcIndex==0 || mcIndex==255); }
-int MarchingCubes::AddTriangles( const double v[Cube::CORNERS] , double iso , Triangle* isoTriangles )
-{
- unsigned char idx;
- int ntriang=0;
- Triangle tri;
-
- idx=GetIndex(v,iso);
-
- /* Cube is entirely in/out of the surface */
- if (!edgeMask[idx]) return 0;
-
- /* Find the vertices where the surface intersects the cube */
- int i,j,ii=1;
- for(i=0;i<12;i++){
- if(edgeMask[idx] & ii){SetVertex(i,v,iso);}
- ii<<=1;
- }
- /* Create the triangle */
- for( i=0 ; triangles[idx][i]!=-1 ; i+=3 )
- {
- for(j=0;j<3;j++){
- tri.p[0][j]=vertexList[triangles[idx][i+0]][j];
- tri.p[1][j]=vertexList[triangles[idx][i+1]][j];
- tri.p[2][j]=vertexList[triangles[idx][i+2]][j];
- }
- isoTriangles[ntriang++]=tri;
- }
- return ntriang;
-}
-
-int MarchingCubes::AddTriangleIndices(const double v[Cube::CORNERS],double iso,int* isoIndices){
- unsigned char idx;
- int ntriang=0;
-
- idx=GetIndex(v,iso);
-
- /* Cube is entirely in/out of the surface */
- if (!edgeMask[idx]) return 0;
-
- /* Create the triangle */
- for(int i=0;triangles[idx][i]!=-1;i+=3){
- for(int j=0;j<3;j++){isoIndices[i+j]=triangles[idx][i+j];}
- ntriang++;
- }
- return ntriang;
-}
-
-void MarchingCubes::SetVertex( int e , const double values[Cube::CORNERS] , double iso )
-{
- double t;
- int o , i1 , i2;
- Cube::FactorEdgeIndex( e , o , i1 , i2 );
- switch( o )
- {
- case 0:
- t = Interpolate( values[ Cube::CornerIndex( 0 , i1 , i2 ) ] - iso , values[ Cube::CornerIndex( 1 , i1 , i2 ) ] - iso );
- vertexList[e][0] = t , vertexList[e][1] = i1 , vertexList[e][2] = i2;
- break;
- case 1:
- t = Interpolate( values[ Cube::CornerIndex( i1 , 0 , i2 ) ] - iso , values[ Cube::CornerIndex( i1 , 1 , i2 ) ] - iso );
- vertexList[e][0] = i1 , vertexList[e][1] = t , vertexList[e][2] = i2;
- break;
- case 2:
- t = Interpolate( values[ Cube::CornerIndex( i1 , i2 , 0 ) ] - iso , values[ Cube::CornerIndex( i1 , i2 , 1 ) ] - iso );
- vertexList[e][0] = i1 , vertexList[e][1] = i2 , vertexList[e][2] = t;
- break;
- }
-}
-double MarchingCubes::Interpolate( double v1 , double v2 ) { return v1/(v1-v2); }
-
-
-///////////////////////////////////
-unsigned char MarchingCubes::GetIndex(const float v[Cube::CORNERS],float iso){
- unsigned char idx=0;
- if (v[Cube::CornerIndex(0,0,0)] < iso) idx |= 1;
- if (v[Cube::CornerIndex(1,0,0)] < iso) idx |= 2;
- if (v[Cube::CornerIndex(1,1,0)] < iso) idx |= 4;
- if (v[Cube::CornerIndex(0,1,0)] < iso) idx |= 8;
- if (v[Cube::CornerIndex(0,0,1)] < iso) idx |= 16;
- if (v[Cube::CornerIndex(1,0,1)] < iso) idx |= 32;
- if (v[Cube::CornerIndex(1,1,1)] < iso) idx |= 64;
- if (v[Cube::CornerIndex(0,1,1)] < iso) idx |= 128;
- return idx;
-}
-unsigned char MarchingCubes::GetFaceIndex( const float values[Cube::CORNERS] , float iso , int faceIndex )
-{
- int i,j,x,y,z;
- unsigned char idx=0;
- double v[2][2];
- Cube::FactorFaceIndex(faceIndex,x,y,z);
- if (x<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(0,i,j)];}}}
- else if (x>0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(1,i,j)];}}}
- else if (y<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,0,j)];}}}
- else if (y>0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,1,j)];}}}
- else if (z<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,j,0)];}}}
- else if (z>0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=values[Cube::CornerIndex(i,j,1)];}}}
- if (v[0][0] < iso) idx |= 1;
- if (v[1][0] < iso) idx |= 2;
- if (v[1][1] < iso) idx |= 4;
- if (v[0][1] < iso) idx |= 8;
- return idx;
-}
-unsigned char MarchingCubes::GetFaceIndex( unsigned char mcIndex , int faceIndex )
-{
- int i,j,x,y,z;
- unsigned char idx=0;
- int v[2][2];
- Cube::FactorFaceIndex(faceIndex,x,y,z);
- if (x<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=mcIndex&(1<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=mcIndex&(1<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=mcIndex&(1<0){for(i=0;i<2;i++){for(j=0;j<2;j++){v[i][j]=mcIndex&(1<
-#include "Geometry.h"
-
-#define NEW_ORDERING 1
-
-class Square
-{
-public:
- const static unsigned int CORNERS=4 , EDGES=4 , FACES=1;
- static int CornerIndex (int x,int y);
- static int AntipodalCornerIndex(int idx);
- static void FactorCornerIndex (int idx,int& x,int& y);
- static int EdgeIndex (int orientation,int i);
- static void FactorEdgeIndex (int idx,int& orientation,int& i);
-
- static int ReflectCornerIndex (int idx,int edgeIndex);
- static int ReflectEdgeIndex (int idx,int edgeIndex);
-
- static void EdgeCorners(int idx,int& c1,int &c2);
-};
-
-class Cube{
-public:
- const static unsigned int CORNERS=8 , EDGES=12 , FACES=6;
-
- static int CornerIndex ( int x , int y , int z );
- static void FactorCornerIndex ( int idx , int& x , int& y , int& z );
- static int EdgeIndex ( int orientation , int i , int j );
- static void FactorEdgeIndex ( int idx , int& orientation , int& i , int &j);
- static int FaceIndex ( int dir , int offSet );
- static int FaceIndex ( int x , int y , int z );
- static void FactorFaceIndex ( int idx , int& x , int &y , int& z );
- static void FactorFaceIndex ( int idx , int& dir , int& offSet );
-
- static int AntipodalCornerIndex ( int idx );
- static int FaceReflectCornerIndex ( int idx , int faceIndex );
- static int FaceReflectEdgeIndex ( int idx , int faceIndex );
- static int FaceReflectFaceIndex ( int idx , int faceIndex );
- static int EdgeReflectCornerIndex ( int idx , int edgeIndex );
- static int EdgeReflectEdgeIndex ( int edgeIndex );
-
- static int FaceAdjacentToEdges ( int eIndex1 , int eIndex2 );
- static void FacesAdjacentToEdge ( int eIndex , int& f1Index , int& f2Index );
-
- static void EdgeCorners( int idx , int& c1 , int &c2 );
- static void FaceCorners( int idx , int& c1 , int &c2 , int& c3 , int& c4 );
-
- static bool IsEdgeCorner( int cIndex , int e );
- static bool IsFaceCorner( int cIndex , int f );
-};
-
-class MarchingSquares
-{
- static double Interpolate(double v1,double v2);
- static void SetVertex(int e,const double values[Square::CORNERS],double iso);
-public:
- const static unsigned int MAX_EDGES=2;
- static const int edgeMask[1<
-struct MemoryInfo
-{
- static size_t Usage( void )
- {
- HANDLE h = GetCurrentProcess();
- PROCESS_MEMORY_COUNTERS pmc;
- return GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ? pmc.WorkingSetSize : 0;
- }
-};
-
-#else // !_WIN32 && !_WIN64
-
-#ifndef __APPLE__ // Linux variants
-
-#include
-#include
-
-class MemoryInfo
-{
- public:
- static size_t Usage(void)
- {
- FILE* f = fopen("/proc/self/stat","rb");
-
- int d;
- long ld;
- unsigned long lu;
- unsigned long long llu;
- char s[1024];
- char c;
-
- int pid;
- unsigned long vm;
-
- int n = fscanf(f, "%d %s %c %d %d %d %d %d %lu %lu %lu %lu %lu %lu %lu %ld %ld %ld %ld %d %ld %llu %lu %ld %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %lu %d %d %lu %lu"
- ,&pid ,s ,&c ,&d ,&d ,&d ,&d ,&d ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&ld ,&ld ,&ld ,&ld ,&d ,&ld ,&llu ,&vm ,&ld ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&lu ,&d ,&d ,&lu ,&lu );
-
- fclose(f);
-/*
-pid %d
-comm %s
-state %c
-ppid %d
-pgrp %d
-session %d
-tty_nr %d
-tpgid %d
-flags %lu
-minflt %lu
-cminflt %lu
-majflt %lu
-cmajflt %lu
-utime %lu
-stime %lu
-cutime %ld
-cstime %ld
-priority %ld
-nice %ld
-0 %ld
-itrealvalue %ld
-starttime %lu
-vsize %lu
-rss %ld
-rlim %lu
-startcode %lu
-endcode %lu
-startstack %lu
-kstkesp %lu
-kstkeip %lu
-signal %lu
-blocked %lu
-sigignore %lu
-sigcatch %lu
-wchan %lu
-nswap %lu
-cnswap %lu
-exit_signal %d
-processor %d
-rt_priority %lu (since kernel 2.5.19)
-policy %lu (since kernel 2.5.19)
-*/
- return vm;
- }
-
-};
-#else // __APPLE__: has no "/proc" pseudo-file system
-
-// Thanks to David O'Gwynn for providing this fix.
-// This comes from a post by Michael Knight:
-//
-// http://miknight.blogspot.com/2005/11/resident-set-size-in-mac-os-x.html
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-
-void getres(task_t task, unsigned long *rss, unsigned long *vs)
-{
- struct task_basic_info t_info;
- mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT;
-
- task_info(task, TASK_BASIC_INFO, (task_info_t)&t_info, &t_info_count);
- *rss = t_info.resident_size;
- *vs = t_info.virtual_size;
-}
-
-class MemoryInfo
-{
- public:
- static size_t Usage(void)
- {
- unsigned long rss, vs, psize;
- task_t task = MACH_PORT_NULL;
-
- if (task_for_pid(current_task(), getpid(), &task) != KERN_SUCCESS)
- abort();
- getres(task, &rss, &vs);
- return rss;
- }
-
-};
-
-#endif // !__APPLE__
-
-#endif // _WIN32 || _WIN64
-
-#endif // MEMORY_USAGE_INCLUDE
diff --git a/Src/MultiGridOctreeData.Evaluation.inl b/Src/MultiGridOctreeData.Evaluation.inl
deleted file mode 100644
index c965effb..00000000
--- a/Src/MultiGridOctreeData.Evaluation.inl
+++ /dev/null
@@ -1,1151 +0,0 @@
-/*
-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.
-*/
-
-template< class Real >
-template< int FEMDegree , BoundaryType BType>
-void Octree< Real >::_Evaluator< FEMDegree , BType >::set( LocalDepth depth )
-{
- static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd;
- static const int RightPointSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart;
-
- BSplineEvaluationData< FEMDegree , BType >::SetEvaluator( evaluator , depth );
- if( depth>0 ) BSplineEvaluationData< FEMDegree , BType >::SetChildEvaluator( childEvaluator , depth-1 );
- int center = ( 1<>1;
-
- // First set the stencils for the current depth
- for( int x=-LeftPointSupportRadius ; x<=RightPointSupportRadius ; x++ ) for( int y=-LeftPointSupportRadius ; y<=RightPointSupportRadius ; y++ ) for( int z=-LeftPointSupportRadius ; z<=RightPointSupportRadius ; z++ )
- {
- int fIdx[] = { center+x , center+y , center+z };
-
- // The cell stencil
- {
- double vv[3] , dv[3];
- for( int dd=0 ; dd( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
-
- //// The face stencil
- for( int f=0 ; f( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
-
- //// The edge stencil
- for( int e=0 ; e( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
-
- //// The corner stencil
- for( int c=0 ; c( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
- }
-
- // Now set the stencils for the parents
- for( int child=0 ; child( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
-
- //// The face stencil
- for( int f=0 ; f( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
-
- //// The edge stencil
- for( int e=0 ; e( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
-
- //// The corner stencil
- for( int c=0 ; c( dv[0] * vv[1] * vv[2] , vv[0] * dv[1] * vv[2] , vv[0] * vv[1] * dv[2] );
- }
- }
- }
- if( _bsData ) delete _bsData;
- _bsData = new BSplineData< FEMDegree , BType >( depth );
-}
-template< class Real >
-template< class V , int FEMDegree , BoundaryType BType >
-V Octree< Real >::_getValue( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , Point3D< Real > p , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator ) const
-{
- static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize;
- static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart;
- static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd;
- static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd;
- static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart;
-
- if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] getValue assumes leaf node\n" );
- V value(0);
-
- while( GetGhostFlag( node ) )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
-
- for( int i=0 ; i _s ; Real _w;
- _startAndWidth( _n , _s , _w );
- int _fIdx[3];
- functionIndex< FEMDegree , BType >( _n , _fIdx );
- for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) );
- value +=
- solution[ _n->nodeData.nodeIndex ] *
- (Real)
- (
- evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) *
- evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) *
- evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- );
- }
- }
- node = node->parent;
- }
-
- LocalDepth d = _localDepth( node );
-
- for( int dd=0 ; dd<3 ; dd++ )
- if ( p[dd]==0 ) p[dd] = (Real)(0.+1e-6);
- else if( p[dd]==1 ) p[dd] = (Real)(1.-1e-6);
-
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
-
- for( int i=0 ; i _s ; Real _w;
- _startAndWidth( _n , _s , _w );
- int _fIdx[3];
- functionIndex< FEMDegree , BType >( _n , _fIdx );
- for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) );
- value +=
- solution[ _n->nodeData.nodeIndex ] *
- (Real)
- (
- evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) *
- evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) *
- evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- );
- }
- }
- if( d>0 )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int i=0 ; i _s ; Real _w;
- _startAndWidth( _n , _s , _w );
- int _fIdx[3];
- functionIndex< FEMDegree , BType >( _n , _fIdx );
- for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) );
- value +=
- coarseSolution[ _n->nodeData.nodeIndex ] *
- (Real)
- (
- evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) *
- evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) *
- evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- );
- }
- }
- }
- }
- return value;
-}
-template< class Real >
-template< int FEMDegree , BoundaryType BType >
-std::pair< Real , Point3D< Real > > Octree< Real >::_getValueAndGradient( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , Point3D< Real > p , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator ) const
-{
- static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize;
- static const int LeftSupportRadius = -BSplineSupportSizes< FEMDegree >::SupportStart;
- static const int RightSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd;
- static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd;
- static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart;
-
- if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] _getValueAndGradient assumes leaf node\n" );
- Real value(0);
- Point3D< Real > gradient;
-
- while( GetGhostFlag( node ) )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
-
- for( int i=0 ; i _s; Real _w;
- _startAndWidth( _n , _s , _w );
- int _fIdx[3];
- functionIndex< FEMDegree , BType >( _n , _fIdx );
- for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) );
- value +=
- solution[ _n->nodeData.nodeIndex ] *
- (Real)
- (
- evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- );
- gradient +=
- Point3D< Real >
- (
- evaluator._bsData->dBaseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) ,
- evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->dBaseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) ,
- evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->dBaseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- ) * solution[ _n->nodeData.nodeIndex ];
- }
- }
- node = node->parent;
- }
-
-
- LocalDepth d = _localDepth( node );
-
- for( int dd=0 ; dd<3 ; dd++ )
- if ( p[dd]==0 ) p[dd] = (Real)(0.+1e-6);
- else if( p[dd]==1 ) p[dd] = (Real)(1.-1e-6);
-
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
-
- for( int i=0 ; i _s ; Real _w;
- _startAndWidth( _n , _s , _w );
- int _fIdx[3];
- functionIndex< FEMDegree , BType >( _n , _fIdx );
- for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) );
- value +=
- solution[ _n->nodeData.nodeIndex ] *
- (Real)
- (
- evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- );
- gradient +=
- Point3D< Real >
- (
- evaluator._bsData->dBaseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) ,
- evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->dBaseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) ,
- evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->dBaseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- ) * solution[ _n->nodeData.nodeIndex ];
- }
- }
- if( d>0 )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int i=0 ; i _s ; Real _w;
- _startAndWidth( _n , _s , _w );
- int _fIdx[3];
- functionIndex< FEMDegree , BType >( _n , _fIdx );
- for( int dd=0 ; dd<3 ; dd++ ) _pIdx[dd] = std::max< int >( 0 , std::min< int >( SupportSize-1 , LeftSupportRadius + (int)floor( ( p[dd]-_s[dd] ) / _w ) ) );
- value +=
- coarseSolution[ _n->nodeData.nodeIndex ] *
- (Real)
- (
- evaluator._bsData->baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- );
- gradient +=
- Point3D< Real >
- (
- evaluator._bsData->dBaseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) ,
- evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData->dBaseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData-> baseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] ) ,
- evaluator._bsData-> baseBSplines[ _fIdx[0] ][ _pIdx[0] ]( p[0] ) * evaluator._bsData-> baseBSplines[ _fIdx[1] ][ _pIdx[1] ]( p[1] ) * evaluator._bsData->dBaseBSplines[ _fIdx[2] ][ _pIdx[2] ]( p[2] )
- ) * coarseSolution[ _n->nodeData.nodeIndex ];
- }
- }
- }
- }
- return std::pair< Real , Point3D< Real > >( value , gradient );
-}
-template< class Real >
-template< class V , int FEMDegree , BoundaryType BType >
-V Octree< Real >::_getCenterValue( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const
-{
- static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize;
- static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd;
- static const int RightPointSupportRadius = - BSplineEvaluationData< FEMDegree , BType >::SupportStart;
-
- if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] getCenterValue assumes leaf node\n" );
- V value(0);
- LocalDepth d = _localDepth( node );
-
- if( isInterior )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
- for( int i=0 ; inodeData.nodeIndex ] * Real( evaluator.cellStencil( i , j , k ) );
- }
- if( d>0 )
- {
- int _corner = int( node - node->parent->children );
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int i=0 ; inodeData.nodeIndex] * Real( evaluator.cellStencils[_corner]( i , j , k ) );
- }
- }
- }
- else
- {
- LocalOffset cIdx;
- _localDepthAndOffset( node , d , cIdx );
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
-
- for( int i=0 ; inodeData.nodeIndex ] *
- Real(
- evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false )
- );
- }
- }
- if( d>0 )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int i=0 ; inodeData.nodeIndex ] *
- Real(
- evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false )
- );
- }
- }
- }
- }
- return value;
-}
-template< class Real >
-template< int FEMDegree , BoundaryType BType >
-std::pair< Real , Point3D< Real > > Octree< Real >::_getCenterValueAndGradient( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const
-{
- static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize;
- static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd;
- static const int RightPointSupportRadius = - BSplineEvaluationData< FEMDegree , BType >::SupportStart;
-
- if( IsActiveNode( node->children ) ) fprintf( stderr , "[WARNING] getCenterValueAndGradient assumes leaf node\n" );
- Real value(0);
- Point3D< Real > gradient;
- LocalDepth d = _localDepth( node );
-
- if( isInterior )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
- for( int i=0 ; inodeData.nodeIndex ];
- gradient += Point3D< Real >( evaluator.dCellStencil( i , j , k ) ) * solution[ n->nodeData.nodeIndex ];
- }
- }
- if( d>0 )
- {
- int _corner = int( node - node->parent->children );
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int i=0 ; inodeData.nodeIndex];
- gradient += Point3D< Real >( evaluator.dCellStencils[_corner]( i , j , k ) ) * coarseSolution[n->nodeData.nodeIndex];
- }
- }
- }
- }
- else
- {
- LocalOffset cIdx;
- _localDepthAndOffset( node , d , cIdx );
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
-
- for( int i=0 ; inodeData.nodeIndex ];
- gradient +=
- Point3D< Real >
- (
- evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , true ) * evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false ) ,
- evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , true ) * evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false ) ,
- evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , true )
- ) * solution[ n->nodeData.nodeIndex ];
- }
- }
- if( d>0 )
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int i=0 ; inodeData.nodeIndex ];
- gradient +=
- Point3D< Real >
- (
- evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , true ) * evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false ) ,
- evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , true ) * evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false ) ,
- evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) * evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) * evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , true )
- ) * coarseSolution[ n->nodeData.nodeIndex ];
- }
- }
- }
- }
- return std::pair< Real , Point3D< Real > >( value , gradient );
-}
-template< class Real >
-template< class V , int FEMDegree , BoundaryType BType >
-V Octree< Real >::_getEdgeValue( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int edge , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const
-{
- static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize;
- static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd;
- static const int RightPointSupportRadius = -BSplineEvaluationData< FEMDegree , BType >::SupportStart;
- V value(0);
- LocalDepth d ; LocalOffset cIdx;
- _localDepthAndOffset( node , d , cIdx );
- int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize;
- int orientation , i1 , i2;
- Cube::FactorEdgeIndex( edge , orientation , i1 , i2 );
- switch( orientation )
- {
- case 0:
- cIdx[1] += i1 , cIdx[2] += i2;
- if( i1 ) startY++ ; else endY--;
- if( i2 ) startZ++ ; else endZ--;
- break;
- case 1:
- cIdx[0] += i1 , cIdx[2] += i2;
- if( i1 ) startX++ ; else endX--;
- if( i2 ) startZ++ ; else endZ--;
- break;
- case 2:
- cIdx[0] += i1 , cIdx[1] += i2;
- if( i1 ) startX++ ; else endX--;
- if( i2 ) startY++ ; else endY--;
- break;
- }
-
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , d );
- for( int x=startX ; xnodeData.nodeIndex ] * evaluator.edgeStencil[edge]( x , y , z );
- else
- {
- LocalDepth _d ; LocalOffset fIdx;
- _localDepthAndOffset( _node , _d , fIdx );
- switch( orientation )
- {
- case 0:
- value +=
- solution[ _node->nodeData.nodeIndex ] *
- Real(
- evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false )
- );
- break;
- case 1:
- value +=
- solution[ _node->nodeData.nodeIndex ] *
- Real(
- evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false )
- );
- break;
- case 2:
- value +=
- solution[ _node->nodeData.nodeIndex ] *
- Real(
- evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false )
- );
- break;
- }
- }
- }
- }
- }
- if( d>0 )
- {
- int _corner = int( node - node->parent->children );
- int _cx , _cy , _cz;
- Cube::FactorCornerIndex( _corner , _cx , _cy , _cz );
- // If the corner/child indices don't match, then the sample position is in the interior of the
- // coarser cell and so the full support resolution should be used.
- switch( orientation )
- {
- case 0:
- if( _cy!=i1 ) startY = 0 , endY = SupportSize;
- if( _cz!=i2 ) startZ = 0 , endZ = SupportSize;
- break;
- case 1:
- if( _cx!=i1 ) startX = 0 , endX = SupportSize;
- if( _cz!=i2 ) startZ = 0 , endZ = SupportSize;
- break;
- case 2:
- if( _cx!=i1 ) startX = 0 , endX = SupportSize;
- if( _cy!=i2 ) startY = 0 , endY = SupportSize;
- break;
- }
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int x=startX ; xnodeData.nodeIndex ] * evaluator.edgeStencils[_corner][edge]( x , y , z );
- else
- {
- LocalDepth _d ; LocalOffset fIdx;
- _localDepthAndOffset( _node , _d , fIdx );
- switch( orientation )
- {
- case 0:
- value +=
- coarseSolution[ _node->nodeData.nodeIndex ] *
- Real(
- evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false )
- );
- break;
- case 1:
- value +=
- coarseSolution[ _node->nodeData.nodeIndex ] *
- Real(
- evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false )
- );
- break;
- case 2:
- value +=
- coarseSolution[ _node->nodeData.nodeIndex ] *
- Real(
- evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false )
- );
- break;
- }
- }
- }
- }
- }
- return Real( value );
-}
-template< class Real >
-template< int FEMDegree , BoundaryType BType >
-std::pair< Real , Point3D< Real > > Octree< Real >::_getEdgeValueAndGradient( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int edge , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const
-{
- static const int SupportSize = BSplineEvaluationData< FEMDegree , BType >::SupportSize;
- static const int LeftPointSupportRadius = BSplineEvaluationData< FEMDegree , BType >::SupportEnd;
- static const int RightPointSupportRadius = -BSplineEvaluationData< FEMDegree , BType >::SupportStart;
- double value = 0;
- Point3D< double > gradient;
- LocalDepth d ; LocalOffset cIdx;
- _localDepthAndOffset( node , d , cIdx );
-
- int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize;
- int orientation , i1 , i2;
- Cube::FactorEdgeIndex( edge , orientation , i1 , i2 );
- switch( orientation )
- {
- case 0:
- cIdx[1] += i1 , cIdx[2] += i2;
- if( i1 ) startY++ ; else endY--;
- if( i2 ) startZ++ ; else endZ--;
- break;
- case 1:
- cIdx[0] += i1 , cIdx[2] += i2;
- if( i1 ) startX++ ; else endX--;
- if( i2 ) startZ++ ; else endZ--;
- break;
- case 2:
- cIdx[0] += i1 , cIdx[1] += i2;
- if( i1 ) startX++ ; else endX--;
- if( i2 ) startY++ ; else endY--;
- break;
- }
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
- for( int x=startX ; xnodeData.nodeIndex ];
- gradient += evaluator.dEdgeStencil[edge]( x , y , z ) * solution[ _node->nodeData.nodeIndex ];
- }
- else
- {
- LocalDepth _d ; LocalOffset fIdx;
- _localDepthAndOffset( _node , _d , fIdx );
-
- double vv[3] , dv[3];
- switch( orientation )
- {
- case 0:
- vv[0] = evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , false );
- vv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false );
- vv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false );
- dv[0] = evaluator.evaluator.centerValue( fIdx[0] , cIdx[0] , true );
- dv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , true );
- dv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , true );
- break;
- case 1:
- vv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false );
- vv[1] = evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , false );
- vv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false );
- dv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , true );
- dv[1] = evaluator.evaluator.centerValue( fIdx[1] , cIdx[1] , true );
- dv[2] = evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , true );
- break;
- case 2:
- vv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false );
- vv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false );
- vv[2] = evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , false );
- dv[0] = evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , true );
- dv[1] = evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , true );
- dv[2] = evaluator.evaluator.centerValue( fIdx[2] , cIdx[2] , true );
- break;
- }
- value += solution[ _node->nodeData.nodeIndex ] * vv[0] * vv[1] * vv[2];
- gradient += Point3D< double >( dv[0]*vv[1]*vv[2] , vv[0]*dv[1]*vv[2] , vv[0]*vv[1]*dv[2] ) * solution[ _node->nodeData.nodeIndex ];
- }
- }
- }
- }
- if( d>0 )
- {
- int _corner = int( node - node->parent->children );
- int _cx , _cy , _cz;
- Cube::FactorCornerIndex( _corner , _cx , _cy , _cz );
- // If the corner/child indices don't match, then the sample position is in the interior of the
- // coarser cell and so the full support resolution should be used.
- switch( orientation )
- {
- case 0:
- if( _cy!=i1 ) startY = 0 , endY = SupportSize;
- if( _cz!=i2 ) startZ = 0 , endZ = SupportSize;
- break;
- case 1:
- if( _cx!=i1 ) startX = 0 , endX = SupportSize;
- if( _cz!=i2 ) startZ = 0 , endZ = SupportSize;
- break;
- case 2:
- if( _cx!=i1 ) startX = 0 , endX = SupportSize;
- if( _cy!=i2 ) startY = 0 , endY = SupportSize;
- break;
- }
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- for( int x=startX ; xnodeData.nodeIndex ];
- gradient += evaluator.dEdgeStencils[_corner][edge]( x , y , z ) * coarseSolution[ _node->nodeData.nodeIndex ];
- }
- else
- {
- LocalDepth _d ; LocalOffset fIdx;
- _localDepthAndOffset( _node , _d , fIdx );
- double vv[3] , dv[3];
- switch( orientation )
- {
- case 0:
- vv[0] = evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , false );
- vv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false );
- vv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false );
- dv[0] = evaluator.childEvaluator.centerValue( fIdx[0] , cIdx[0] , true );
- dv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , true );
- dv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , true );
- break;
- case 1:
- vv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false );
- vv[1] = evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , false );
- vv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false );
- dv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , true );
- dv[1] = evaluator.childEvaluator.centerValue( fIdx[1] , cIdx[1] , true );
- dv[2] = evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , true );
- break;
- case 2:
- vv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false );
- vv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false );
- vv[2] = evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , false );
- dv[0] = evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , true );
- dv[1] = evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , true );
- dv[2] = evaluator.childEvaluator.centerValue( fIdx[2] , cIdx[2] , true );
- break;
- }
- value += coarseSolution[ _node->nodeData.nodeIndex ] * vv[0] * vv[1] * vv[2];
- gradient += Point3D< double >( dv[0]*vv[1]*vv[2] , vv[0]*dv[1]*vv[2] , vv[0]*vv[1]*dv[2] ) * coarseSolution[ _node->nodeData.nodeIndex ];
- }
- }
- }
- }
- return std::pair< Real , Point3D< Real > >( Real( value ) , Point3D< Real >( gradient ) );
-}
-
-template< class Real >
-template< class V , int FEMDegree , BoundaryType BType >
-V Octree< Real >::_getCornerValue( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int corner , const DenseNodeData< V , FEMDegree >& solution , const DenseNodeData< V , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const
-{
- static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize;
- static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd;
- static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart;
-
- V value(0);
- LocalDepth d ; LocalOffset cIdx;
- _localDepthAndOffset( node , d , cIdx );
-
- int cx , cy , cz;
- int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize;
- Cube::FactorCornerIndex( corner , cx , cy , cz );
- cIdx[0] += cx , cIdx[1] += cy , cIdx[2] += cz;
- {
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
- if( cx==0 ) endX--;
- else startX++;
- if( cy==0 ) endY--;
- else startY++;
- if( cz==0 ) endZ--;
- else startZ++;
- if( isInterior )
- for( int x=startX ; xnodeData.nodeIndex ] * Real( evaluator.cornerStencil[corner]( x , y , z ) );
- }
- else
- for( int x=startX ; xnodeData.nodeIndex ] *
- Real(
- evaluator.evaluator.cornerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.evaluator.cornerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.evaluator.cornerValue( fIdx[2] , cIdx[2] , false )
- );
- }
- }
- }
- if( d>0 )
- {
- int _corner = int( node - node->parent->children );
- int _cx , _cy , _cz;
- Cube::FactorCornerIndex( _corner , _cx , _cy , _cz );
- // If the corner/child indices don't match, then the sample position is in the interior of the
- // coarser cell and so the full support resolution should be used.
- if( cx!=_cx ) startX = 0 , endX = SupportSize;
- if( cy!=_cy ) startY = 0 , endY = SupportSize;
- if( cz!=_cz ) startZ = 0 , endZ = SupportSize;
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- if( isInterior )
- for( int x=startX ; xnodeData.nodeIndex ] * Real( evaluator.cornerStencils[_corner][corner]( x , y , z ) );
- }
- else
- for( int x=startX ; xnodeData.nodeIndex ] *
- Real(
- evaluator.childEvaluator.cornerValue( fIdx[0] , cIdx[0] , false ) *
- evaluator.childEvaluator.cornerValue( fIdx[1] , cIdx[1] , false ) *
- evaluator.childEvaluator.cornerValue( fIdx[2] , cIdx[2] , false )
- );
- }
- }
- }
- return Real( value );
-}
-template< class Real >
-template< int FEMDegree , BoundaryType BType >
-std::pair< Real , Point3D< Real > > Octree< Real >::_getCornerValueAndGradient( const ConstPointSupportKey< FEMDegree >& neighborKey , const TreeOctNode* node , int corner , const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , const _Evaluator< FEMDegree , BType >& evaluator , bool isInterior ) const
-{
- static const int SupportSize = BSplineSupportSizes< FEMDegree >::SupportSize;
- static const int LeftPointSupportRadius = BSplineSupportSizes< FEMDegree >::SupportEnd;
- static const int RightPointSupportRadius = - BSplineSupportSizes< FEMDegree >::SupportStart;
-
- double value = 0;
- Point3D< double > gradient;
- LocalDepth d ; LocalOffset cIdx;
- _localDepthAndOffset( node , d , cIdx );
-
- int cx , cy , cz;
- int startX = 0 , endX = SupportSize , startY = 0 , endY = SupportSize , startZ = 0 , endZ = SupportSize;
- Cube::FactorCornerIndex( corner , cx , cy , cz );
- cIdx[0] += cx , cIdx[1] += cy , cIdx[2] += cz;
- {
- if( cx==0 ) endX--;
- else startX++;
- if( cy==0 ) endY--;
- else startY++;
- if( cz==0 ) endZ--;
- else startZ++;
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node );
- if( isInterior )
- for( int x=startX ; xnodeData.nodeIndex ] * evaluator.cornerStencil[corner]( x , y , z ) , gradient += evaluator.dCornerStencil[corner]( x , y , z ) * solution[ _node->nodeData.nodeIndex ];
- }
- else
- for( int x=startX ; xnodeData.nodeIndex ] * v[0] * v[1] * v[2];
- gradient += Point3D< double >( dv[0]*v[1]*v[2] , v[0]*dv[1]*v[2] , v[0]*v[1]*dv[2] ) * solution[ _node->nodeData.nodeIndex ];
- }
- }
- }
- if( d>0 )
- {
- int _corner = int( node - node->parent->children );
- int _cx , _cy , _cz;
- Cube::FactorCornerIndex( _corner , _cx , _cy , _cz );
- if( cx!=_cx ) startX = 0 , endX = SupportSize;
- if( cy!=_cy ) startY = 0 , endY = SupportSize;
- if( cz!=_cz ) startZ = 0 , endZ = SupportSize;
- const typename TreeOctNode::ConstNeighbors< SupportSize >& neighbors = _neighbors< LeftPointSupportRadius , RightPointSupportRadius >( neighborKey , node->parent );
- if( isInterior )
- for( int x=startX ; xnodeData.nodeIndex ] * evaluator.cornerStencils[_corner][corner]( x , y , z ) , gradient += evaluator.dCornerStencils[_corner][corner]( x , y , z ) * coarseSolution[ _node->nodeData.nodeIndex ];
- }
- else
- for( int x=startX ; xnodeData.nodeIndex ] * v[0] * v[1] * v[2];
- gradient += Point3D< double >( dv[0]*v[1]*v[2] , v[0]*dv[1]*v[2] , v[0]*v[1]*dv[2] ) * coarseSolution[ _node->nodeData.nodeIndex ];
- }
- }
- }
- return std::pair< Real , Point3D< Real > >( Real( value ) , Point3D< Real >( gradient ) );
-}
-template< class Real >
-template< int Degree , BoundaryType BType >
-Octree< Real >::MultiThreadedEvaluator< Degree , BType >::MultiThreadedEvaluator( const Octree< Real >* tree , const DenseNodeData< Real , Degree >& coefficients , int threads ) : _coefficients( coefficients ) , _tree( tree )
-{
- _threads = std::max< int >( 1 , threads );
- _neighborKeys.resize( _threads );
- _coarseCoefficients = _tree->template coarseCoefficients< Real , Degree , BType >( _coefficients );
- _evaluator.set( _tree->_maxDepth );
- for( int t=0 ; t<_threads ; t++ ) _neighborKeys[t].set( tree->_localToGlobal( _tree->_maxDepth ) );
-}
-template< class Real >
-template< int Degree , BoundaryType BType >
-Real Octree< Real >::MultiThreadedEvaluator< Degree , BType >::value( Point3D< Real > p , int thread , const TreeOctNode* node )
-{
- if( !node ) node = _tree->leaf( p );
- ConstPointSupportKey< Degree >& nKey = _neighborKeys[thread];
- nKey.getNeighbors( node );
- return _tree->template _getValue< Real , Degree >( nKey , node , p , _coefficients , _coarseCoefficients , _evaluator );
-}
-template< class Real >
-template< int Degree , BoundaryType BType >
-std::pair< Real , Point3D< Real > > Octree< Real >::MultiThreadedEvaluator< Degree , BType >::valueAndGradient( Point3D< Real > p , int thread , const TreeOctNode* node )
-{
- if( !node ) node = _tree->leaf( p );
- ConstPointSupportKey< Degree >& nKey = _neighborKeys[thread];
- nKey.getNeighbors( node );
- return _tree->template _getValueAndGradient< Degree >( nKey , node , p , _coefficients , _coarseCoefficients , _evaluator );
-}
diff --git a/Src/MultiGridOctreeData.IsoSurface.inl b/Src/MultiGridOctreeData.IsoSurface.inl
deleted file mode 100644
index 2be26dc5..00000000
--- a/Src/MultiGridOctreeData.IsoSurface.inl
+++ /dev/null
@@ -1,1106 +0,0 @@
-/*
-Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without modification,
-are permitted provided that the following conditions are met:
-
-Redistributions of source code must retain the above copyright notice, this list of
-conditions and the following disclaimer. Redistributions in binary form must reproduce
-the above copyright notice, this list of conditions and the following disclaimer
-in the documentation and/or other materials provided with the distribution.
-
-Neither the name of the Johns Hopkins University nor the names of its contributors
-may be used to endorse or promote products derived from this software without specific
-prior written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
-EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
-OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
-SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
-INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
-TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
-BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
-ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
-DAMAGE.
-*/
-
-#include "Octree.h"
-#include "MyTime.h"
-#include "MemoryUsage.h"
-#include "MAT.h"
-
-template< class Real >
-template< class Vertex >
-Octree< Real >::_SliceValues< Vertex >::_SliceValues( void )
-{
- _oldCCount = _oldECount = _oldFCount = _oldNCount = 0;
- cornerValues = NullPointer( Real ) ; cornerGradients = NullPointer( Point3D< Real > ) ; cornerSet = NullPointer( char );
- edgeKeys = NullPointer( long long ) ; edgeSet = NullPointer( char );
- faceEdges = NullPointer( _FaceEdges ) ; faceSet = NullPointer( char );
- mcIndices = NullPointer( char );
-}
-template< class Real >
-template< class Vertex >
-Octree< Real >::_SliceValues< Vertex >::~_SliceValues( void )
-{
- _oldCCount = _oldECount = _oldFCount = _oldNCount = 0;
- FreePointer( cornerValues ) ; FreePointer( cornerGradients ) ; FreePointer( cornerSet );
- FreePointer( edgeKeys ) ; FreePointer( edgeSet );
- FreePointer( faceEdges ) ; FreePointer( faceSet );
- FreePointer( mcIndices );
-}
-template< class Real >
-template< class Vertex >
-void Octree< Real >::_SliceValues< Vertex >::reset( bool nonLinearFit )
-{
- faceEdgeMap.clear() , edgeVertexMap.clear() , vertexPairMap.clear();
-
- if( _oldNCount0 ) mcIndices = AllocPointer< char >( _oldNCount );
- }
- if( _oldCCount0 )
- {
- cornerValues = AllocPointer< Real >( _oldCCount );
- if( nonLinearFit ) cornerGradients = AllocPointer< Point3D< Real > >( _oldCCount );
- cornerSet = AllocPointer< char >( _oldCCount );
- }
- }
- if( _oldECount( _oldECount );
- edgeSet = AllocPointer< char >( _oldECount );
- }
- if( _oldFCount( _oldFCount );
- faceSet = AllocPointer< char >( _oldFCount );
- }
-
- if( sliceData.cCount>0 ) memset( cornerSet , 0 , sizeof( char ) * sliceData.cCount );
- if( sliceData.eCount>0 ) memset( edgeSet , 0 , sizeof( char ) * sliceData.eCount );
- if( sliceData.fCount>0 ) memset( faceSet , 0 , sizeof( char ) * sliceData.fCount );
-}
-template< class Real >
-template< class Vertex >
-Octree< Real >::_XSliceValues< Vertex >::_XSliceValues( void )
-{
- _oldECount = _oldFCount = 0;
- edgeKeys = NullPointer( long long ) ; edgeSet = NullPointer( char );
- faceEdges = NullPointer( _FaceEdges ) ; faceSet = NullPointer( char );
-}
-template< class Real >
-template< class Vertex >
-Octree< Real >::_XSliceValues< Vertex >::~_XSliceValues( void )
-{
- _oldECount = _oldFCount = 0;
- FreePointer( edgeKeys ) ; FreePointer( edgeSet );
- FreePointer( faceEdges ) ; FreePointer( faceSet );
-}
-template< class Real >
-template< class Vertex >
-void Octree< Real >::_XSliceValues< Vertex >::reset( void )
-{
- faceEdgeMap.clear() , edgeVertexMap.clear() , vertexPairMap.clear();
-
- if( _oldECount( _oldECount );
- edgeSet = AllocPointer< char >( _oldECount );
- }
- if( _oldFCount( _oldFCount );
- faceSet = AllocPointer< char >( _oldFCount );
- }
- if( xSliceData.eCount>0 ) memset( edgeSet , 0 , sizeof( char ) * xSliceData.eCount );
- if( xSliceData.fCount>0 ) memset( faceSet , 0 , sizeof( char ) * xSliceData.fCount );
-}
-
-template< class Real >
-template< int FEMDegree , BoundaryType BType , int WeightDegree , int ColorDegree , class Vertex >
-void Octree< Real >::getMCIsoSurface( const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , const DenseNodeData< Real , FEMDegree >& solution , Real isoValue , CoredMeshData< Vertex >& mesh , bool nonLinearFit , bool addBarycenter , bool polygonMesh )
-{
- if( FEMDegree==1 && nonLinearFit ) fprintf( stderr , "[WARNING] First order B-Splines do not support non-linear interpolation\n" ) , nonLinearFit = false;
-
- BSplineData< ColorDegree , BOUNDARY_NEUMANN >* colorBSData = NULL;
- if( colorData ) colorBSData = new BSplineData< ColorDegree , BOUNDARY_NEUMANN >( _maxDepth );
- DenseNodeData< Real , FEMDegree > coarseSolution( _sNodesEnd(_maxDepth-1) );
- memset( &coarseSolution[0] , 0 , sizeof(Real)*_sNodesEnd( _maxDepth-1) );
-#pragma omp parallel for num_threads( threads )
- for( int i=_sNodesBegin(0) ; i<_sNodesEnd(_maxDepth-1) ; i++ ) coarseSolution[i] = solution[i];
- for( LocalDepth d=1 ; d<_maxDepth ; d++ ) _upSample< Real , FEMDegree , BType >( d , coarseSolution );
- memoryUsage();
-
- std::vector< _Evaluator< FEMDegree , BType > > evaluators( _maxDepth+1 );
- for( LocalDepth d=0 ; d<=_maxDepth ; d++ ) evaluators[d].set( d );
-
- int vertexOffset = 0;
-
- std::vector< _SlabValues< Vertex > > slabValues( _maxDepth+1 );
-
- // Initialize the back slice
- for( LocalDepth d=_maxDepth ; d>=0 ; d-- )
- {
- _sNodes.setSliceTableData ( slabValues[d]. sliceValues(0). sliceData , _localToGlobal( d ) , 0 + _localInset( d ) , threads );
- _sNodes.setSliceTableData ( slabValues[d]. sliceValues(1). sliceData , _localToGlobal( d ) , 1 + _localInset( d ) , threads );
- _sNodes.setXSliceTableData( slabValues[d].xSliceValues(0).xSliceData , _localToGlobal( d ) , 0 + _localInset( d ) , threads );
- slabValues[d].sliceValues (0).reset( nonLinearFit );
- slabValues[d].sliceValues (1).reset( nonLinearFit );
- slabValues[d].xSliceValues(0).reset( );
- }
- for( LocalDepth d=_maxDepth ; d>=0 ; d-- )
- {
- // Copy edges from finer
- if( d<_maxDepth ) _copyFinerSliceIsoEdgeKeys( d , 0 , slabValues , threads );
- _setSliceIsoCorners( solution , coarseSolution , isoValue , d , 0 , slabValues , evaluators[d] , threads );
- _setSliceIsoVertices< WeightDegree , ColorDegree >( colorBSData , densityWeights , colorData , isoValue , d , 0 , vertexOffset , mesh , slabValues , threads );
- _setSliceIsoEdges( d , 0 , slabValues , threads );
- }
-
- // Iterate over the slices at the finest level
- for( int slice=0 ; slice<( 1<<_maxDepth ) ; slice++ )
- {
- // Process at all depths that contain this slice
- LocalDepth d ; int o;
- for( d=_maxDepth , o=slice+1 ; d>=0 ; d-- , o>>=1 )
- {
- // Copy edges from finer (required to ensure we correctly track edge cancellations)
- if( d<_maxDepth )
- {
- _copyFinerSliceIsoEdgeKeys( d , o , slabValues , threads );
- _copyFinerXSliceIsoEdgeKeys( d , o-1 , slabValues , threads );
- }
-
- // Set the slice values/vertices
- _setSliceIsoCorners( solution , coarseSolution , isoValue , d , o , slabValues , evaluators[d] , threads );
- _setSliceIsoVertices< WeightDegree , ColorDegree >( colorBSData , densityWeights , colorData , isoValue , d , o , vertexOffset , mesh , slabValues , threads );
- _setSliceIsoEdges( d , o , slabValues , threads );
-
- // Set the cross-slice edges
- _setXSliceIsoVertices< WeightDegree , ColorDegree >( colorBSData , densityWeights , colorData , isoValue , d , o-1 , vertexOffset , mesh , slabValues , threads );
- _setXSliceIsoEdges( d , o-1 , slabValues , threads );
-
- // Add the triangles
- _setIsoSurface( d , o-1 , slabValues[d].sliceValues(o-1) , slabValues[d].sliceValues(o) , slabValues[d].xSliceValues(o-1) , mesh , polygonMesh , addBarycenter , vertexOffset , threads );
-
- if( o&1 ) break;
- }
-
- for( d=_maxDepth , o=slice+1 ; d>=0 ; d-- , o>>=1 )
- {
- // Initialize for the next pass
- if( o<(1<<(d+1)) )
- {
- _sNodes.setSliceTableData( slabValues[d].sliceValues(o+1).sliceData , _localToGlobal( d ) , o+1 + _localInset( d ) , threads );
- _sNodes.setXSliceTableData( slabValues[d].xSliceValues(o).xSliceData , _localToGlobal( d ) , o + _localInset( d ) , threads );
- slabValues[d].sliceValues(o+1).reset( nonLinearFit );
- slabValues[d].xSliceValues(o).reset();
- }
- if( o&1 ) break;
- }
- }
- memoryUsage();
- if( colorBSData ) delete colorBSData;
-}
-
-
-template< class Real >
-template< class Vertex , int FEMDegree , BoundaryType BType >
-void Octree< Real >::_setSliceIsoCorners( const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , Real isoValue , LocalDepth depth , int slice , std::vector< _SlabValues< Vertex > >& slabValues , const _Evaluator< FEMDegree , BType >& evaluator , int threads )
-{
- if( slice>0 ) _setSliceIsoCorners( solution , coarseSolution , isoValue , depth , slice , 1 , slabValues , evaluator , threads );
- if( slice<(1<
-template< class Vertex , int FEMDegree , BoundaryType BType >
-void Octree< Real >::_setSliceIsoCorners( const DenseNodeData< Real , FEMDegree >& solution , const DenseNodeData< Real , FEMDegree >& coarseSolution , Real isoValue , LocalDepth depth , int slice , int z , std::vector< _SlabValues< Vertex > >& slabValues , const struct _Evaluator< FEMDegree , BType >& evaluator , int threads )
-{
- typename Octree::template _SliceValues< Vertex >& sValues = slabValues[depth].sliceValues( slice );
- std::vector< ConstPointSupportKey< FEMDegree > > neighborKeys( std::max< int >( 1 , threads ) );
- for( size_t i=0 ; i& neighborKey = neighborKeys[ omp_get_thread_num() ];
- TreeOctNode* leaf = _sNodes.treeNodes[i];
- if( !IsActiveNode( leaf->children ) )
- {
- const typename SortedTreeNodes::SquareCornerIndices& cIndices = sValues.sliceData.cornerIndices( leaf );
-
- bool isInterior = _isInteriorlySupported< FEMDegree >( leaf->parent );
- neighborKey.getNeighbors( leaf );
-
- for( int x=0 ; x<2 ; x++ ) for( int y=0 ; y<2 ; y++ )
- {
- int cc = Cube::CornerIndex( x , y , z );
- int fc = Square::CornerIndex( x , y );
- int vIndex = cIndices[fc];
- if( !sValues.cornerSet[vIndex] )
- {
- if( sValues.cornerGradients )
- {
- std::pair< Real , Point3D< Real > > p = _getCornerValueAndGradient( neighborKey , leaf , cc , solution , coarseSolution , evaluator , isInterior );
- sValues.cornerValues[vIndex] = p.first , sValues.cornerGradients[vIndex] = p.second;
- }
- else sValues.cornerValues[vIndex] = _getCornerValue( neighborKey , leaf , cc , solution , coarseSolution , evaluator , isInterior );
- sValues.cornerSet[vIndex] = 1;
- }
- squareValues[fc] = sValues.cornerValues[ vIndex ];
- TreeOctNode* node = leaf;
- LocalDepth _depth = depth;
- int _slice = slice;
- while( _isValidSpaceNode( node->parent ) && (node-node->parent->children)==cc )
- {
- node = node->parent , _depth-- , _slice >>= 1;
- typename Octree::template _SliceValues< Vertex >& _sValues = slabValues[_depth].sliceValues( _slice );
- const typename SortedTreeNodes::SquareCornerIndices& _cIndices = _sValues.sliceData.cornerIndices( node );
- int _vIndex = _cIndices[fc];
- _sValues.cornerValues[_vIndex] = sValues.cornerValues[vIndex];
- if( _sValues.cornerGradients ) _sValues.cornerGradients[_vIndex] = sValues.cornerGradients[vIndex];
- _sValues.cornerSet[_vIndex] = 1;
- }
- }
- sValues.mcIndices[ i - sValues.sliceData.nodeOffset ] = MarchingSquares::GetIndex( squareValues , isoValue );
- }
- }
-}
-
-template< class Real >
-template< int WeightDegree , int ColorDegree , BoundaryType BType , class Vertex >
-void Octree< Real >::_setSliceIsoVertices( const BSplineData< ColorDegree , BType >* colorBSData , const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Point3D< Real > , Real > , ColorDegree >* colorData , Real isoValue , LocalDepth depth , int slice , int& vOffset , CoredMeshData< Vertex >& mesh , std::vector< _SlabValues< Vertex > >& slabValues , int threads )
-{
- if( slice>0 ) _setSliceIsoVertices< WeightDegree , ColorDegree >( colorBSData , densityWeights , colorData , isoValue , depth , slice , 1 , vOffset , mesh , slabValues , threads );
- if( slice<(1<