From 9230936ab47e9f3d0dbf0e732d3c57fa44bca059 Mon Sep 17 00:00:00 2001
From: mkazhdan <misha@cs.jhu.edu>
Date: Fri, 19 Jul 2019 00:25:18 -0400
Subject: [PATCH] Version 12.00: Added functionality to fit a function to set
 of sample values (PointInterpolant) Added functionality to sample a function
 at prescribed locations (AdaptiveTreeVisualization)

---
 AdaptiveSolvers.sln               |  73 +++
 Makefile                          |  16 +
 PointInterpolant.vcxproj          | 167 ++++++
 README.md                         | 145 +++++-
 Src/AdaptiveTreeVisualization.cpp |  65 ++-
 Src/EDTInHeat.cpp                 |   2 +-
 Src/FEMTree.WeightedSamples.inl   |   2 +-
 Src/FEMTree.h                     |  14 +-
 Src/FEMTree.inl                   |  93 ++--
 Src/Image.h                       |   2 +-
 Src/Ply.h                         |   6 +-
 Src/PlyFile.h                     |   6 +-
 Src/PointInterpolant.cpp          | 820 ++++++++++++++++++++++++++++++
 Src/PointStream.inl               |   4 +-
 Src/PointStreamData.h             |   6 +-
 Src/PoissonRecon.cpp              |  28 +-
 Src/PreProcessor.h                |   2 +-
 Src/SSDRecon.cpp                  |  27 +-
 Src/SparseMatrix.h                |   6 +-
 19 files changed, 1391 insertions(+), 93 deletions(-)
 create mode 100644 PointInterpolant.vcxproj
 create mode 100644 Src/PointInterpolant.cpp

diff --git a/AdaptiveSolvers.sln b/AdaptiveSolvers.sln
index 26b9c71..b696e93 100644
--- a/AdaptiveSolvers.sln
+++ b/AdaptiveSolvers.sln
@@ -102,32 +102,105 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Inline Files", "Inline File
 EndProject
 Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "ChunkPLY", "ChunkPLY.vcxproj", "{95B998F8-A73B-40AB-AF35-C38F2CB89480}"
 EndProject
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "PointInterpolant", "PointInterpolant.vcxproj", "{D5516BBE-9613-4AE3-BE4D-E754995F944C}"
+EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|x64 = Debug|x64
+		Debug|x86 = Debug|x86
 		Release|x64 = Release|x64
+		Release|x86 = Release|x86
 	EndGlobalSection
 	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Debug|x64.ActiveCfg = Debug|x64
+		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Debug|x64.Build.0 = Debug|x64
+		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Debug|x86.ActiveCfg = Debug|Win32
+		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Debug|x86.Build.0 = Debug|Win32
 		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Release|x64.ActiveCfg = Release|x64
 		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Release|x64.Build.0 = Release|x64
+		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Release|x86.ActiveCfg = Release|Win32
+		{D3D173AB-D306-4179-BEC4-95CE1B14E647}.Release|x86.Build.0 = Release|Win32
+		{477765F8-C16A-406B-807F-1302DAE74EBA}.Debug|x64.ActiveCfg = Debug|x64
+		{477765F8-C16A-406B-807F-1302DAE74EBA}.Debug|x64.Build.0 = Debug|x64
+		{477765F8-C16A-406B-807F-1302DAE74EBA}.Debug|x86.ActiveCfg = Debug|Win32
+		{477765F8-C16A-406B-807F-1302DAE74EBA}.Debug|x86.Build.0 = Debug|Win32
 		{477765F8-C16A-406B-807F-1302DAE74EBA}.Release|x64.ActiveCfg = Release|x64
 		{477765F8-C16A-406B-807F-1302DAE74EBA}.Release|x64.Build.0 = Release|x64
+		{477765F8-C16A-406B-807F-1302DAE74EBA}.Release|x86.ActiveCfg = Release|Win32
+		{477765F8-C16A-406B-807F-1302DAE74EBA}.Release|x86.Build.0 = Release|Win32
+		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Debug|x64.ActiveCfg = Debug|x64
+		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Debug|x64.Build.0 = Debug|x64
+		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Debug|x86.ActiveCfg = Debug|Win32
+		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Debug|x86.Build.0 = Debug|Win32
 		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Release|x64.ActiveCfg = Release|x64
 		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Release|x64.Build.0 = Release|x64
+		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Release|x86.ActiveCfg = Release|Win32
+		{B5899B32-FAC2-477E-99AA-86736B97F2FC}.Release|x86.Build.0 = Release|Win32
+		{0BE55595-4080-4265-82AF-51108EC956B2}.Debug|x64.ActiveCfg = Debug|x64
+		{0BE55595-4080-4265-82AF-51108EC956B2}.Debug|x64.Build.0 = Debug|x64
+		{0BE55595-4080-4265-82AF-51108EC956B2}.Debug|x86.ActiveCfg = Debug|Win32
+		{0BE55595-4080-4265-82AF-51108EC956B2}.Debug|x86.Build.0 = Debug|Win32
 		{0BE55595-4080-4265-82AF-51108EC956B2}.Release|x64.ActiveCfg = Release|x64
 		{0BE55595-4080-4265-82AF-51108EC956B2}.Release|x64.Build.0 = Release|x64
+		{0BE55595-4080-4265-82AF-51108EC956B2}.Release|x86.ActiveCfg = Release|Win32
+		{0BE55595-4080-4265-82AF-51108EC956B2}.Release|x86.Build.0 = Release|Win32
+		{E09D1615-1036-4F43-909A-2E751FC947F6}.Debug|x64.ActiveCfg = Debug|x64
+		{E09D1615-1036-4F43-909A-2E751FC947F6}.Debug|x64.Build.0 = Debug|x64
+		{E09D1615-1036-4F43-909A-2E751FC947F6}.Debug|x86.ActiveCfg = Debug|Win32
+		{E09D1615-1036-4F43-909A-2E751FC947F6}.Debug|x86.Build.0 = Debug|Win32
 		{E09D1615-1036-4F43-909A-2E751FC947F6}.Release|x64.ActiveCfg = Release|x64
 		{E09D1615-1036-4F43-909A-2E751FC947F6}.Release|x64.Build.0 = Release|x64
+		{E09D1615-1036-4F43-909A-2E751FC947F6}.Release|x86.ActiveCfg = Release|Win32
+		{E09D1615-1036-4F43-909A-2E751FC947F6}.Release|x86.Build.0 = Release|Win32
+		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Debug|x64.ActiveCfg = Debug|x64
+		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Debug|x64.Build.0 = Debug|x64
+		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Debug|x86.ActiveCfg = Debug|Win32
+		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Debug|x86.Build.0 = Debug|Win32
 		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Release|x64.ActiveCfg = Release|x64
 		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Release|x64.Build.0 = Release|x64
 		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Release|x64.Deploy.0 = Release|x64
+		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Release|x86.ActiveCfg = Release|Win32
+		{3B94F9AD-35E8-4E1E-B176-AAA091E4C3CC}.Release|x86.Build.0 = Release|Win32
+		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|x64.ActiveCfg = Debug|x64
+		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|x64.Build.0 = Debug|x64
+		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|x86.ActiveCfg = Debug|Win32
+		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Debug|x86.Build.0 = Debug|Win32
 		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|x64.ActiveCfg = Release|x64
 		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|x64.Build.0 = Release|x64
+		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|x86.ActiveCfg = Release|Win32
+		{99BEAFED-8DB9-4B7D-A0BE-5186158193FE}.Release|x86.Build.0 = Release|Win32
+		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|x64.ActiveCfg = Debug|x64
+		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|x64.Build.0 = Debug|x64
+		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|x86.ActiveCfg = Debug|Win32
+		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Debug|x86.Build.0 = Debug|Win32
 		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|x64.ActiveCfg = Release|x64
 		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|x64.Build.0 = Release|x64
+		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|x86.ActiveCfg = Release|Win32
+		{46F87D0E-C53A-4F95-AB48-A5DBA8014340}.Release|x86.Build.0 = Release|Win32
+		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Debug|x64.ActiveCfg = Debug|x64
+		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Debug|x64.Build.0 = Debug|x64
+		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Debug|x86.ActiveCfg = Debug|Win32
+		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Debug|x86.Build.0 = Debug|Win32
 		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Release|x64.ActiveCfg = Release|x64
 		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Release|x64.Build.0 = Release|x64
+		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Release|x86.ActiveCfg = Release|Win32
+		{742064B3-CEBB-4AF8-B43C-D213EB1C6D7A}.Release|x86.Build.0 = Release|Win32
+		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Debug|x64.ActiveCfg = Debug|x64
+		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Debug|x64.Build.0 = Debug|x64
+		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Debug|x86.ActiveCfg = Debug|Win32
+		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Debug|x86.Build.0 = Debug|Win32
 		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Release|x64.ActiveCfg = Release|x64
 		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Release|x64.Build.0 = Release|x64
+		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Release|x86.ActiveCfg = Release|Win32
+		{95B998F8-A73B-40AB-AF35-C38F2CB89480}.Release|x86.Build.0 = Release|Win32
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Debug|x64.ActiveCfg = Debug|x64
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Debug|x64.Build.0 = Debug|x64
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Debug|x86.ActiveCfg = Debug|Win32
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Debug|x86.Build.0 = Debug|Win32
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Release|x64.ActiveCfg = Release|x64
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Release|x64.Build.0 = Release|x64
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Release|x86.ActiveCfg = Release|Win32
+		{D5516BBE-9613-4AE3-BE4D-E754995F944C}.Release|x86.Build.0 = Release|Win32
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE
diff --git a/Makefile b/Makefile
index 839ebd8..00464dd 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,6 @@
 PR_TARGET=PoissonRecon
 SR_TARGET=SSDRecon
+PI_TARGET=PointInterpolant
 ST_TARGET=SurfaceTrimmer
 EH_TARGET=EDTInHeat
 IS_TARGET=ImageStitching
@@ -7,6 +8,7 @@ AV_TARGET=AdaptiveTreeVisualization
 CP_TARGET=ChunkPLY
 PR_SOURCE=PoissonRecon.cpp
 SR_SOURCE=SSDRecon.cpp
+PI_SOURCE=PointInterpolant.cpp
 ST_SOURCE=SurfaceTrimmer.cpp
 EH_SOURCE=EDTInHeat.cpp
 IS_SOURCE=ImageStitching.cpp
@@ -54,6 +56,7 @@ MD=mkdir
 
 PR_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(PR_SOURCE))))
 SR_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(SR_SOURCE))))
+PI_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(PI_SOURCE))))
 ST_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(ST_SOURCE))))
 EH_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(EH_SOURCE))))
 IS_OBJECTS=$(addprefix $(BIN), $(addsuffix .o, $(basename $(IS_SOURCE))))
@@ -66,6 +69,7 @@ all: LFLAGS += $(LFLAGS_RELEASE)
 all: make_dir
 all: $(BIN)$(PR_TARGET)
 all: $(BIN)$(SR_TARGET)
+all: $(BIN)$(PI_TARGET)
 all: $(BIN)$(ST_TARGET)
 all: $(BIN)$(EH_TARGET)
 all: $(BIN)$(IS_TARGET)
@@ -77,6 +81,7 @@ debug: LFLAGS += $(LFLAGS_DEBUG)
 debug: make_dir
 debug: $(BIN)$(PR_TARGET)
 debug: $(BIN)$(SR_TARGET)
+debug: $(BIN)$(PI_TARGET)
 debug: $(BIN)$(ST_TARGET)
 debug: $(BIN)$(EH_TARGET)
 debug: $(BIN)$(IS_TARGET)
@@ -93,6 +98,11 @@ ssdrecon: LFLAGS += $(LFLAGS_RELEASE)
 ssdrecon: make_dir
 ssdrecon: $(BIN)$(SR_TARGET)
 
+pointinterpolant: CFLAGS += $(CFLAGS_RELEASE)
+pointinterpolant: LFLAGS += $(LFLAGS_RELEASE)
+pointinterpolant: make_dir
+pointinterpolant: $(BIN)$(PI_TARGET)
+
 surfacetrimmer: CFLAGS += $(CFLAGS_RELEASE)
 surfacetrimmer: LFLAGS += $(LFLAGS_RELEASE)
 surfacetrimmer: make_dir
@@ -121,6 +131,7 @@ chunkply: $(BIN)$(CP_TARGET)
 clean:
 	rm -rf $(BIN)$(PR_TARGET)
 	rm -rf $(BIN)$(SR_TARGET)
+	rm -rf $(BIN)$(PI_TARGET)
 	rm -rf $(BIN)$(ST_TARGET)
 	rm -rf $(BIN)$(EH_TARGET)
 	rm -rf $(BIN)$(IS_TARGET)
@@ -128,6 +139,7 @@ clean:
 	rm -rf $(BIN)$(CP_TARGET)
 	rm -rf $(PR_OBJECTS)
 	rm -rf $(SR_OBJECTS)
+	rm -rf $(PI_OBJECTS)
 	rm -rf $(ST_OBJECTS)
 	rm -rf $(EH_OBJECTS)
 	rm -rf $(IS_OBJECTS)
@@ -147,6 +159,10 @@ $(BIN)$(SR_TARGET): $(SR_OBJECTS)
 	cd PNG  && make
 	$(CXX) -o $@ $(SR_OBJECTS) -L$(BIN) $(LFLAGS) -ljpeg -lmypng -lz
 
+$(BIN)$(PI_TARGET): $(PI_OBJECTS)
+	cd PNG  && make
+	$(CXX) -o $@ $(PI_OBJECTS) -L$(BIN) $(LFLAGS) -ljpeg -lmypng -lz
+
 $(BIN)$(ST_TARGET): $(ST_OBJECTS)
 	$(CXX) -o $@ $(ST_OBJECTS) $(LFLAGS)
 
diff --git a/PointInterpolant.vcxproj b/PointInterpolant.vcxproj
new file mode 100644
index 0000000..b4da247
--- /dev/null
+++ b/PointInterpolant.vcxproj
@@ -0,0 +1,167 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project DefaultTargets="Build" ToolsVersion="15.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup Label="ProjectConfigurations">
+    <ProjectConfiguration Include="Debug|Win32">
+      <Configuration>Debug</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|Win32">
+      <Configuration>Release</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Debug|x64">
+      <Configuration>Debug</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|x64">
+      <Configuration>Release</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="Src\PointInterpolant.cpp" />
+  </ItemGroup>
+  <PropertyGroup Label="Globals">
+    <VCProjectVersion>15.0</VCProjectVersion>
+    <ProjectGuid>{D5516BBE-9613-4AE3-BE4D-E754995F944C}</ProjectGuid>
+    <Keyword>Win32Proj</Keyword>
+    <RootNamespace>PointInterpolant</RootNamespace>
+    <WindowsTargetPlatformVersion>10.0.17763.0</WindowsTargetPlatformVersion>
+    <ProjectName>PointInterpolant</ProjectName>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v141</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v141</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v141</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v141</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>MultiByte</CharacterSet>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
+  <ImportGroup Label="ExtensionSettings">
+  </ImportGroup>
+  <ImportGroup Label="Shared">
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <PropertyGroup Label="UserMacros" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <LinkIncremental>false</LinkIncremental>
+    <OutDir>$(SolutionDir)\Bin\$(Platform)\$(Configuration)\</OutDir>
+    <IntDir>$(SolutionDir)\Obj\$(TargetName)\$(Platform)\$(Configuration)\</IntDir>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <LinkIncremental>true</LinkIncremental>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <LinkIncremental>true</LinkIncremental>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <LinkIncremental>false</LinkIncremental>
+  </PropertyGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <ClCompile>
+      <PrecompiledHeader>NotUsing</PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>MaxSpeed</Optimization>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>NOMINMAX;WIN32;NDEBUG;_CONSOLE;_CRT_SECURE_NO_DEPRECATE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>false</ConformanceMode>
+      <PrecompiledHeaderFile>pch.h</PrecompiledHeaderFile>
+      <AdditionalIncludeDirectories>.;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <EnableEnhancedInstructionSet>AdvancedVectorExtensions2</EnableEnhancedInstructionSet>
+      <OpenMPSupport>true</OpenMPSupport>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalLibraryDirectories>$(OutDir)</AdditionalLibraryDirectories>
+      <AdditionalDependencies>ZLIB.lib;JPEG.lib;PNG.lib;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)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <ClCompile>
+      <PrecompiledHeader>Use</PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>Disabled</Optimization>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+      <PrecompiledHeaderFile>pch.h</PrecompiledHeaderFile>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <ClCompile>
+      <PrecompiledHeader>Use</PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>Disabled</Optimization>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+      <PrecompiledHeaderFile>pch.h</PrecompiledHeaderFile>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <ClCompile>
+      <PrecompiledHeader>Use</PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>MaxSpeed</Optimization>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <SDLCheck>true</SDLCheck>
+      <PreprocessorDefinitions>WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <ConformanceMode>true</ConformanceMode>
+      <PrecompiledHeaderFile>pch.h</PrecompiledHeaderFile>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
+  <ImportGroup Label="ExtensionTargets">
+  </ImportGroup>
+</Project>
\ No newline at end of file
diff --git a/README.md b/README.md
index 05b0626..f31b85b 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-<center><h2>Adaptive Multigrid Solvers (Version 11.02)</h2></center>
+<center><h2>Adaptive Multigrid Solvers (Version 12.00)</h2></center>
 <center>
 <a href="#LINKS">links</a>
 <a href="#EXECUTABLES">executables</a>
@@ -27,10 +27,11 @@ This code-base was born from the Poisson Surface Reconstruction code. It has evo
 <a href="http://www.cs.jhu.edu/~misha/MyPapers/ToG13.pdf">[Kazhdan and Hoppe, 2013]</a>
 <br>
 <b>Executables: </b>
-<a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version11.02/AdaptiveSolvers.x64.zip">Win64</a><br>
+<a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version12.00/AdaptiveSolvers.x64.zip">Win64</a><br>
 <b>Source Code:</b>
-<a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version11.02/AdaptiveSolvers.zip">ZIP</a> <a href="https://github.com/mkazhdan/PoissonRecon">GitHub</a><br>
+<a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version12.00/AdaptiveSolvers.zip">ZIP</a> <a href="https://github.com/mkazhdan/PoissonRecon">GitHub</a><br>
 <b>Older Versions:</b>
+<a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version11.02/">V11.02</a>,
 <a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version11.01/">V11.01</a>,
 <a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version11.00/">V11.00</a>,
 <a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version10.07/">V10.07</a>,
@@ -336,6 +337,103 @@ individual components of the surface reconstructor.
 </ul>
 
 
+<ul>
+<dl>
+<DETAILS>
+<SUMMARY>
+<font size="+1"><b>PointInterpolant</b></font>:
+Fits a function to a set of sample values (and gradients)
+</SUMMARY>
+
+<dt><b>--in</b> &lt;<i>input sample positions and values</i>&gt;</dt>
+<dd> This string is the name of the file from which the positions and values will be read.<br>
+The file should be an ascii file with groups of <I>Dim</I>+1 or 2*<I>Dim</I>+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).<br>
+No information about the number of samples should be specified.</dd>
+
+<dt>[<b>--dim</b> &lt;<i>dimension of the samples</i>&gt;]</dt>
+<dd> This integerl value is the dimension of the samples.<BR>
+The default value is 2.<br></dd>
+
+<dt>[<b>--tree</b> &lt;<i>output octree and coefficients</i>&gt;]</dt>
+<dd> This string is the name of the file to which the the octree and function coefficients are to be written.</dd>
+
+<dt>[<b>--grid</b> &lt;<i>output grid</i>&gt;]</dt>
+<dd> 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^<i>d</i>,
+and the next 4 x 2^<i>d</i> x 2^<i>d</i> x ... bytes corresponding to the (single precision) floating point values
+of the implicit function.</dd>
+
+<dt>[<b>--degree</b> &lt;<i>B-spline degree</i>&gt;]</dt>
+<dd> 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).<br>
+The default value for this parameter is 2.</dt>
+
+<dt>[<b>--bType</b> &lt;<i>boundary type</i>&gt;]</dt>
+<dd> This integer specifies the boundary type for the finite elements. Valid values are:
+<ul>
+<li> <b>1</b>: Free boundary constraints</li>
+<li> <b>2</b>: Dirichlet boundary constraints</li>
+<li> <b>3</b>: Neumann boundary constraints</li>
+</ul>
+The default value for this parameter is 1 (free).
+
+<dt>[<b>--depth</b> &lt;<i>reconstruction depth</i>&gt;]</dt>
+<dd> This integer is the maximum depth of the tree that will be used for surface reconstruction.
+Running at depth <i>d</i> corresponds to solving on a grid whose resolution is no larger than
+2^<i>d</i> x 2^<i>d</i> x ... Note that since the reconstructor adapts the octree to the
+sampling density, the specified reconstruction depth is only an upper bound.<br>
+The default value for this parameter is 8.</dd>
+
+<dt>[<b>--width</b> &lt;<i>finest cell width</i>&gt;]</dt>
+<dd> This floating point value specifies the target width of the finest level octree cells.<br>
+This parameter is ignored if the <B>--depth</B> is also specified.</dd>
+
+<dt>[<b>--scale</b> &lt;<i>scale factor</i>&gt;]</dt>
+<dd> This floating point value specifies the ratio between the diameter of the cube used for reconstruction
+and the diameter of the samples' bounding cube.<br>
+The default value is 1.1.</dd>
+
+<dt>[<b>--valueWeight</b> &lt;<i>value interpolation weight</i>&gt;]</dt>
+<dd> This floating point value specifies the importance that interpolation of the samples' values
+is given in the fitting of the function.<br>
+The default value for this parameter is 1000.</dd>
+
+<dt>[<b>--gradientWeight</b> &lt;<i>gradient interpolation weight</i>&gt;]</dt>
+<dd> This floating point value specifies the importance that interpolation of the samples' gradients
+is given in the fitting of the function.<br>
+The default value for this parameter is 1.<BR>
+This value is ignored unless gradient interpolation is specified.</dd>
+
+<dt>[<b>--lapWeight</b> &lt;<i>Laplacian weight</i>&gt;]</dt>
+<dd> This floating point value specifies the importance that Laplacian regularization
+is given in the fitting of the function.<br>
+The default value for this parameter is 0.</dd>
+
+<dt>[<b>--biLapWeight</b> &lt;<i>bi-Laplacian weight</i>&gt;]</dt>
+<dd> This floating point value specifies the importance that bi-Laplacian regularization
+is given in the fitting of the function.<br>
+The default value for this parameter is 1.</dd>
+
+<dt>[<b>--iters</b> &lt;<i>GS iters</i>&gt;]</dt>
+<dd> This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.<br>
+The default value for this parameter is 8.</dd>
+
+<dt>[<b>--useGradients</b>]</dt>
+<dd> Enabling this flag indicates that the input file contains gradients as well as sample values.</dd>
+
+<dt>[<b>--performance</b>]</dt>
+<dd> Enabling this flag provides running time and peak memory usage at the end of the execution.</dd>
+
+<dt>[<b>--verbose</b>]</dt>
+<dd> Enabling this flag provides a more verbose description of the running times and memory usages of
+individual components of the surface reconstructor.</dd>
+
+</DETAILS>
+</dl>
+</ul>
+
+
 <ul>
 <dl>
 <DETAILS>
@@ -504,9 +602,15 @@ The default value for this parameter is equal to the numer of (virtual) processo
 Extracts iso-surfaces and a sampling on a regular grid from an implicit function represented over an adapted tree
 </SUMMARY>
 <dt><b>--in</b> &lt;<i>input tree and coefficients</i>&gt;
-</dt><dd> This string is the name of the file from which the tree and implicit functions coefficients are to be read. 
+</dt><dd> This string is the name of the file from which the tree and implicit functions coefficients are to be read.</dd>
 
-</dd><dt>[<b>--grid</b> &lt;<i>output value grid</i>&gt;]
+<dt>[<b>--samples</b> &lt;<i>input sample positions</i>&gt;]</dt>
+<dd> This string is the name of the file from which sampling positions are to be read.<BR>
+The file should be an ascii file with groups of <I>Dim</I> white space delimited, numbers giving the coordinates of the sampling  points' position.<br>
+No information about the number of samples should be specified.</dd>
+</dd>
+
+<dt>[<b>--grid</b> &lt;<i>output value grid</i>&gt;]
 </dt><dd> This string is the name of the file to which the sampling of the implicit along a regular grid will be written.<BR>
 The file is written out in binary, with the first 4 bytes corresponding to the (integer) sampling resolution, <i>R</i>,
 and the next 4 x <I>R</I>^<i>D</i> bytes corresponding to the (single precision) floating point values of the implicit function. (Here, <i>D</I> is the dimension.)
@@ -629,6 +733,36 @@ which partitions the reconstruction into 11 pieces.
 </dl>
 </ul>
 
+<ul>
+<dl>
+<DETAILS>
+<SUMMARY>
+<font size="+1"><b>PointInterpolant / AdaptiveTreeVisualization</b></font>
+</SUMMARY>
+For testing purposes, a pair of point-sets is provided:
+<ol>
+
+<li> <a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/quadratic.2D.fitting.samples"><b>fitting samples</b></a>:
+A set of 1000 random 2D samples from within the square [-1,1,]x[-1,1] along with the evaluation of the quadratic <i>f(x,y)=x*x+y*y</i> at each sample point (represented in ASCII format).
+<LI> <a href="http://www.cs.jhu.edu/~misha/Code/PoissonRecon/quadratic.2D.evaluation.samples"><b>evaluation samples</b></a>:
+A set of 4 2D positions at which the fit function is to be evaluated (represented in ASCII format).
+</ol>
+
+The function fitting the input samples can be by calling the point interpolant:
+<blockquote><code>% PointInterpolant --in quadratic.2D.fitting.samples --out quadratic.2D.tree --dim 2</code></blockquote>
+Then, the reconstructed function can be evaluated at the evaluation samples by calling the adaptive tree visualization:
+<blockquote><code>% AdaptiveTreeVisualization --in quadratic.2D.tree --samples quadratic.2D.evaluation.samples</code></blockquote>
+This will output the evaluation positions and values:
+<blockquote><CODE>
+0 0 1.33836e-05<br>
+0.5 0 0.25001<br>
+0.5 0.5 0.500006<br>
+2 2 nan
+</CODE></blockquote>
+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.
+</DETAILS>
+</dl>
+</ul>
 
 <ul>
 <dl>
@@ -871,4 +1005,3 @@ Similarly, to reduce compilation times, support for specific degrees can be remo
 <hr>
 <a name="SUPPORT"><b>SUPPORT</b></a><br>
 This work genersouly supported by NSF grants #0746039 and #1422325.
-
diff --git a/Src/AdaptiveTreeVisualization.cpp b/Src/AdaptiveTreeVisualization.cpp
index eabcafa..ac6a32d 100644
--- a/Src/AdaptiveTreeVisualization.cpp
+++ b/Src/AdaptiveTreeVisualization.cpp
@@ -42,6 +42,7 @@ DAMAGE.
 
 cmdLineParameter< char* >
 	In( "in" ) ,
+	Samples( "samples" ) ,
 	OutMesh( "mesh" ) ,
 	OutGrid( "grid" );
 
@@ -66,6 +67,7 @@ cmdLineParameter< float >
 cmdLineReadable* params[] =
 {
 	&In , 
+	&Samples ,
 	&OutMesh , &NonManifold , &PolygonMesh , &FlipOrientation , &ASCII , &NonLinearFit , &IsoValue ,
 	&OutGrid , &PrimalGrid ,
 	&Threads ,
@@ -81,6 +83,7 @@ void ShowUsage( char* ex )
 {
 	printf( "Usage: %s\n" , ex );
 	printf( "\t --%s <input tree>\n" , In.name );
+	printf( "\t[--%s sample positions>]\n" , Samples.name );
 	printf( "\t[--%s <ouput triangle mesh>]\n" , OutMesh.name );
 	printf( "\t[--%s <ouput grid>]\n" , OutGrid.name );
 	printf( "\t[--%s <num threads>=%d]\n" , Threads.name , Threads.value );
@@ -158,7 +161,7 @@ void WriteGrid( ConstPointer( Real ) values , int res , const char *fileName )
 }
 
 template< unsigned int Dim , class Real , unsigned int FEMSig >
-void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
+void _Execute( const FEMTree< Dim , Real >* tree , XForm< Real , Dim+1 > xForm , FILE* fp )
 {
 	ThreadPool::Init( (ThreadPool::ParallelType)ParallelType.value , Threads.value );
 	static const unsigned int Degree = FEMSignature< FEMSig >::Degree;
@@ -166,6 +169,37 @@ void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
 
 	coefficients.read( fp );
 
+	// Evaluate at the sample positions
+	if( Samples.set )
+	{
+		InputPointStream< Real , Dim > *pointStream;
+		char* ext = GetFileExtension( Samples.value );
+		if     ( !strcasecmp( ext , "bpts" ) ) pointStream = new BinaryInputPointStream< Real , Dim >( Samples.value );
+		else if( !strcasecmp( ext , "ply"  ) ) pointStream = new    PLYInputPointStream< Real , Dim >( Samples.value );
+		else                                   pointStream = new  ASCIIInputPointStream< Real , Dim >( Samples.value );
+		delete[] ext;
+		typename FEMTree< Dim , Real >::template MultiThreadedEvaluator< IsotropicUIntPack< Dim , FEMSig > , 0 > evaluator( tree , coefficients );
+		static const unsigned int CHUNK_SIZE = 1024;
+		Point< Real , Dim > points[ CHUNK_SIZE ];
+		Real values[ CHUNK_SIZE ];
+		size_t pointsRead;
+		while( ( pointsRead=pointStream->nextPoints( points , CHUNK_SIZE ) ) )
+		{
+			ThreadPool::Parallel_for( 0 , pointsRead , [&]( unsigned int thread , size_t j )
+			{
+				Point< Real , Dim > p = xForm * points[j];
+				bool inBounds = true;
+				for( int d=0 ; d<Dim ; d++ ) if( p[d]<0 || p[d]>1 ) inBounds = false;
+				if( inBounds ) values[j] = evaluator.values( xForm * points[j] , thread )[0];
+				else           values[j] = (Real)nan( "" );
+			}
+			);
+			for( int j=0 ; j<pointsRead ; j++ ) printf( "%g %g %g\n" , points[j][0] , points[j][1] , values[j] );
+		}
+
+		delete pointStream;
+	}
+
 	// Output the grid
 	if( OutGrid.set )
 	{
@@ -196,7 +230,7 @@ void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
 		if( Verbose.set ) printf( "Vertices / Polygons: %llu / %llu\n" , (unsigned long long)( mesh.outOfCorePointCount()+mesh.inCorePoints.size() ) , (unsigned long long)mesh.polygonCount() );
 
 		std::vector< std::string > comments;
-		if( !PlyWritePolygons< Vertex , node_index_type , Real , Dim >( OutMesh.value , &mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , comments , XForm< Real , Dim+1 >::Identity() ) )
+		if( !PlyWritePolygons< Vertex , node_index_type , Real , Dim >( OutMesh.value , &mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , comments , xForm.inverse() ) )
 			ERROR_OUT( "Could not write mesh to: " , OutMesh.value );
 	}
 }
@@ -205,7 +239,8 @@ void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
 template< unsigned int Dim , class Real >
 void Execute( FILE* fp , int degree , BoundaryType bType )
 {
-	FEMTree< Dim , Real > tree( fp , MEMORY_ALLOCATOR_BLOCK_SIZE );
+	XForm< Real , Dim+1 > xForm;
+	FEMTree< Dim , Real > tree( fp , xForm , MEMORY_ALLOCATOR_BLOCK_SIZE );
 
 	if( Verbose.set ) printf( "Leaf Nodes / Active Nodes / Ghost Nodes: %llu / %llu / %llu\n" , (unsigned long long)tree.leaves() , (unsigned long long)tree.nodes() , (unsigned long long)tree.ghostNodes() );
 
@@ -215,10 +250,10 @@ void Execute( FILE* fp , int degree , BoundaryType bType )
 	{
 		switch( degree )
 		{
-			case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
-			case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
-			case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
-			case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
+			case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
+			case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
+			case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
+			case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
 			default: ERROR_OUT( "Only B-Splines of degree 1 - 4 are supported" );
 		}
 	}
@@ -227,10 +262,10 @@ void Execute( FILE* fp , int degree , BoundaryType bType )
 	{
 		switch( degree )
 		{
-			case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
-			case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
-			case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
-			case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
+			case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
+			case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
+			case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
+			case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
 			default: ERROR_OUT( "Only B-Splines of degree 1 - 4 are supported" );
 		}
 	}
@@ -239,10 +274,10 @@ void Execute( FILE* fp , int degree , BoundaryType bType )
 	{
 		switch( degree )
 		{
-			case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
-			case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
-			case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
-			case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
+			case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
+			case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
+			case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
+			case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
 			default: ERROR_OUT( "Only B-Splines of degree 1 - 4 are supported" );
 		}
 	}
diff --git a/Src/EDTInHeat.cpp b/Src/EDTInHeat.cpp
index b5a2a73..40cb341 100644
--- a/Src/EDTInHeat.cpp
+++ b/Src/EDTInHeat.cpp
@@ -545,7 +545,7 @@ void _Execute( int argc , char* argv[] )
 			if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Out.value );
 			FEMTree< Dim , Real >::WriteParameter( fp );
 			DenseNodeData< Real , IsotropicUIntPack< Dim , FEMSig > >::WriteSignatures( fp );
-			tree.write( fp );
+			tree.write( fp , xForm );
 			edtSolution.write( fp );
 			fclose( fp );
 		}
diff --git a/Src/FEMTree.WeightedSamples.inl b/Src/FEMTree.WeightedSamples.inl
index a3e82b2..b263dda 100644
--- a/Src/FEMTree.WeightedSamples.inl
+++ b/Src/FEMTree.WeightedSamples.inl
@@ -116,7 +116,7 @@ void FEMTree< Dim , Real >::_getSampleDepthAndWeight( const DensityEstimator< We
 	{
 		Real oldWeight , newWeight;
 		oldWeight = newWeight = weight;
-		while( newWeight<(Real)1. && temp->parent )
+		while( newWeight<(Real)1. && _localDepth(temp) )
 		{
 			temp=temp->parent;
 			oldWeight = newWeight;
diff --git a/Src/FEMTree.h b/Src/FEMTree.h
index bc8d110..71702c8 100644
--- a/Src/FEMTree.h
+++ b/Src/FEMTree.h
@@ -2129,13 +2129,13 @@ public:
 	static void ResetLocalMemoryUsage( void ){ _LocalMemoryUsage = 0; }
 	static double MemoryUsage( void );
 	FEMTree( size_t blockSize );
-	FEMTree( FILE* fp , size_t blockSize );
+	FEMTree( FILE* fp , XForm< Real , Dim+1 > &xForm , size_t blockSize );
 	~FEMTree( void )
 	{
 		if( _tree ) for( int c=0 ; c<(1<<Dim) ; c++ ) _tree[c].cleanChildren( !nodeAllocators.size() );
 		for( size_t i=0 ; i<nodeAllocators.size() ; i++ ) delete nodeAllocators[i];
 	}
-	void write( FILE* fp ) const;
+	void write( FILE* fp , XForm< Real , Dim+1 > xForm ) const;
 	static void WriteParameter( FILE* fp )
 	{
 		FEMTreeRealType realType;
@@ -2153,17 +2153,19 @@ public:
 	template< unsigned int Radius , class IsThickenNode , class ... DenseOrSparseNodeData > void thicken( IsThickenNode F , DenseOrSparseNodeData* ... data ){ thicken< Radius , Radius >( F , data ... ); }
 	template< unsigned int DensityDegree >
 	typename FEMTree::template DensityEstimator< DensityDegree >* setDensityEstimator( const std::vector< PointSample >& samples , LocalDepth splatDepth , Real samplesPerNode , int coDimension );
-	template< unsigned int ... NormalSigs , unsigned int DensityDegree , class Data >
+	template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
+	SparseNodeData< OutData , UIntPack< DataSigs ... > > setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData & , Real & ) > ConversionAndBiasFunction );
+	template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
 #if defined(_WIN32) || defined(_WIN64)
-	SparseNodeData< Point< Real , Dim > , UIntPack< NormalSigs ... > > setNormalField( UIntPack< NormalSigs ... > , const std::vector< PointSample >& samples , const std::vector< Data >& normalData , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< Real ( Real ) > BiasFunction = []( Real ){ return 0.f; } );
+	SparseNodeData< OutData , UIntPack< DataSigs ... > > setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData& ) > ConversionFunction , std::function< Real ( InData ) > BiasFunction = []( InData ){ return 0.f; } );
 #else // !_WIN32 && !_WIN64
-	SparseNodeData< Point< Real , Dim > , UIntPack< NormalSigs ... > > setNormalField( UIntPack< NormalSigs ... > , const std::vector< PointSample >& samples , const std::vector< Data >& normalData , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< Real ( Real ) > BiasFunction = []( Real ){ return (Real)0; } );
+	SparseNodeData< OutData , UIntPack< DataSigs ... > > setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData& ) > ConversionFunction , std::function< Real ( InData ) > BiasFunction = []( InData ){ return (Real)0; } );
 #endif // _WIN32 || _WIN64
 
 	template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
 	SparseNodeData< Data , IsotropicUIntPack< Dim , DataSig > > setSingleDepthDataField( const std::vector< PointSample >& samples , const std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density );
 	template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
-	SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > setDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest=false );
+	SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > setMultiDepthDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest=false );
 	template< unsigned int MaxDegree , class HasDataFunctor , class ... DenseOrSparseNodeData > void finalizeForMultigrid( LocalDepth fullDepth , const HasDataFunctor F , DenseOrSparseNodeData* ... data );
 
 	template< unsigned int ... FEMSigs > DenseNodeData< Real , UIntPack< FEMSigs ... > > initDenseNodeData( UIntPack< FEMSigs ... > ) const;
diff --git a/Src/FEMTree.inl b/Src/FEMTree.inl
index e6cd90a..e4c2310 100644
--- a/Src/FEMTree.inl
+++ b/Src/FEMTree.inl
@@ -73,7 +73,7 @@ template< unsigned int Dim , class Real > FEMTree< Dim , Real >::FEMTree( size_t
 	memset( _refinableSigs , -1 , sizeof( _refinableSigs ) );
 }
 template< unsigned int Dim , class Real >
-FEMTree< Dim , Real >::FEMTree( FILE* fp , size_t blockSize ) : _nodeInitializer( *this )
+FEMTree< Dim , Real >::FEMTree( FILE* fp , XForm< Real , Dim+1 > &xForm , size_t blockSize ) : _nodeInitializer( *this )
 {
 	if( blockSize )
 	{
@@ -87,6 +87,7 @@ FEMTree< Dim , Real >::FEMTree( FILE* fp , size_t blockSize ) : _nodeInitializer
 	Allocator< FEMTreeNode > *nodeAllocator = nodeAllocators.size() ? nodeAllocators[0] : NULL;
 	if( fp )
 	{
+		if( fread( xForm.coords , sizeof( Real ) , (Dim+1)*(Dim+1) , fp )!=(Dim+1)*(Dim+1) ) ERROR_OUT( "Failed to read transform" );
 		if( fread( &_depthOffset , sizeof( int ) , 1 , fp )!=1 ) ERROR_OUT( "Failed to read depth offset" );
 		_tree = FEMTreeNode::NewBrood( nodeAllocator , _nodeInitializer );
 		_tree->read( fp , nodeAllocator , _nodeInitializer );
@@ -113,8 +114,9 @@ FEMTree< Dim , Real >::FEMTree( FILE* fp , size_t blockSize ) : _nodeInitializer
 		_depthOffset = 0;
 	}
 }
-template< unsigned int Dim , class Real > void FEMTree< Dim , Real >::write( FILE* fp ) const
+template< unsigned int Dim , class Real > void FEMTree< Dim , Real >::write( FILE* fp , XForm< Real , Dim+1 > xForm ) const
 {
+	fwrite( xForm.coords , sizeof( Real ) , (Dim+1)*(Dim+1) , fp );
 	fwrite( &_depthOffset , sizeof( int ) , 1 , fp );
 	_tree->write( fp );
 }
@@ -321,68 +323,81 @@ typename FEMTree< Dim , Real >::template DensityEstimator< DensityDegree >* FEMT
 	MemoryUsage();
 	return _density;
 }
+
 template< unsigned int Dim , class Real >
-template< unsigned int ... NormalSigs , unsigned int DensityDegree , class Data >
-SparseNodeData< Point< Real , Dim > , UIntPack< NormalSigs ... > > FEMTree< Dim , Real >::setNormalField( UIntPack< NormalSigs ... > , const std::vector< PointSample >& samples , const std::vector< Data >& normalData , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< Real ( Real ) > BiasFunction )
+template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
+SparseNodeData< OutData , UIntPack< DataSigs ... > > FEMTree< Dim , Real >::setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData& ) > ConversionFunction , std::function< Real ( InData ) > BiasFunction )
+{
+	std::function< bool ( InData , OutData & , Real & ) > ConversionAndBiasFunction = [&]( InData in , OutData &out , Real &bias )
+	{
+		if( ConversionFunction( in , out ) )
+		{
+			bias = BiasFunction( in );
+			return true;
+		}
+		else return false;
+	};
+	return setDataField( UIntPack< DataSigs ... >() , samples , data , density , pointWeightSum , ConversionAndBiasFunction );
+}
+template< unsigned int Dim , class Real >
+template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
+SparseNodeData< OutData , UIntPack< DataSigs ... > > FEMTree< Dim , Real >::setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData & , Real & ) > ConversionAndBiasFunction )
 {
 	LocalDepth maxDepth = _spaceRoot->maxDepth();
 	typedef PointSupportKey< IsotropicUIntPack< Dim , DensityDegree > > DensityKey;
-	typedef UIntPack< FEMSignature< NormalSigs >::Degree ... > NormalDegrees;
-	typedef PointSupportKey< UIntPack< FEMSignature< NormalSigs >::Degree ... > > NormalKey;
+	typedef UIntPack< FEMSignature< DataSigs >::Degree ... > DataDegrees;
+	typedef PointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > > DataKey;
 	std::vector< DensityKey > densityKeys( ThreadPool::NumThreads() );
-	std::vector<  NormalKey >  normalKeys( ThreadPool::NumThreads() );
-	bool oneKey = DensityDegree==NormalDegrees::Min() && DensityDegree==NormalDegrees::Max();
+	std::vector<    DataKey >    dataKeys( ThreadPool::NumThreads() );
+	bool oneKey = DensityDegree==DataDegrees::Min() && DensityDegree==DataDegrees::Max();
 	for( size_t i=0 ; i<densityKeys.size() ; i++ ) densityKeys[i].set( _localToGlobal( maxDepth ) );
-	if( !oneKey ) for( size_t i=0 ; i<normalKeys.size() ; i++ ) normalKeys[i].set( _localToGlobal( maxDepth ) );
+	if( !oneKey ) for( size_t i=0 ; i<dataKeys.size() ; i++ ) dataKeys[i].set( _localToGlobal( maxDepth ) );
 	Real weightSum = 0;
 	pointWeightSum = 0;
-	SparseNodeData< Point< Real , Dim > , UIntPack< NormalSigs ... > > normalField;
+	SparseNodeData< OutData , UIntPack< DataSigs ... > > dataField;
 	Real _pointWeightSum = 0;
 	ThreadPool::Parallel_for( 0 , samples.size() , [&]( unsigned int thread , size_t i )
 	{
 		DensityKey& densityKey = densityKeys[ thread ];
-		NormalKey& normalKey = normalKeys[ thread ];
+		DataKey& dataKey = dataKeys[ thread ];
 		const ProjectiveData< Point< Real , Dim > , Real >& sample = samples[i].sample;
 		if( sample.weight>0 )
 		{
-			Point< Real , Dim > p = sample.data / sample.weight , n = normalData[i].template data<0>();
-			Real l = (Real)Length( n );
-			// It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
-			if( !l ) return;
-			Real confidence = l / sample.weight;
-			n *= sample.weight / l;
-			Real depthBias = BiasFunction( confidence );
-			AddAtomic( weightSum , sample.weight );
-			if( !_InBounds(p) )
-			{
-				WARN( "Point sample is out of bounds" );
-				return;
-			}
-			Allocator< FEMTreeNode > *nodeAllocator = nodeAllocators.size() ? nodeAllocators[ thread ] : NULL;
+			Point< Real , Dim > p = sample.data / sample.weight;
+			InData in = data[i] / sample.weight;
+			OutData out;
 
-#if defined( __GNUC__ ) && __GNUC__ < 5
-#warning "you've got me gcc version<5"
-			if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , Point< Real , Dim > >( nodeAllocator , *density , p , n , normalField , densityKey , oneKey ? *( (NormalKey*)&densityKey ) : normalKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
-#else // !__GNUC__ || __GNUC__ >=5
-			if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , Point< Real , Dim > , NormalSigs ... >( nodeAllocator , *density , p , n , normalField , densityKey , oneKey ? *( (NormalKey*)&densityKey ) : normalKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
-#endif // __GNUC__ || __GNUC__ < 4
-			else
+			Real depthBias;
+			if( !_InBounds(p) ) WARN( "Point sample is out of bounds" );
+			else if( ConversionAndBiasFunction( in , out , depthBias ) )
 			{
-				Real width = (Real)( 1.0 / ( 1<<maxDepth ) );
+				AddAtomic( weightSum , sample.weight );
+				out *= sample.weight;
+				Allocator< FEMTreeNode > *nodeAllocator = nodeAllocators.size() ? nodeAllocators[ thread ] : NULL;
 #if defined( __GNUC__ ) && __GNUC__ < 5
-#warning "you've got me gcc version<5"
-				_splatPointData< true , true , Point< Real , Dim > >( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , n / (Real)pow( width , Dim ) , normalField , oneKey ? *( (NormalKey*)&densityKey ) : normalKey );
+				#warning "you've got me gcc version<5"
+					if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , OutData >( nodeAllocator , *density , p , out , dataField , densityKey , oneKey ? *( (DataKey*)&densityKey ) : dataKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
 #else // !__GNUC__ || __GNUC__ >=5
-				_splatPointData< true , true , Point< Real , Dim > , NormalSigs ... >( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , n / (Real)pow( width , Dim ) , normalField , oneKey ? *( (NormalKey*)&densityKey ) : normalKey );
+				if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , OutData , DataSigs ... >( nodeAllocator , *density , p , out , dataField , densityKey , oneKey ? *( (DataKey*)&densityKey ) : dataKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
 #endif // __GNUC__ || __GNUC__ < 4
-				AddAtomic( _pointWeightSum , sample.weight );
+				else
+				{
+					Real width = (Real)( 1.0 / ( 1<<maxDepth ) );
+#if defined( __GNUC__ ) && __GNUC__ < 5
+					#warning "you've got me gcc version<5"
+						_splatPointData< true , true , OutData >( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , out / (Real)pow( width , Dim ) , dataField , oneKey ? *( (DataKey*)&densityKey ) : dataKey );
+#else // !__GNUC__ || __GNUC__ >=5
+					_splatPointData< true , true , OutData , DataSigs ... >( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , out / (Real)pow( width , Dim ) , dataField , oneKey ? *( (DataKey*)&densityKey ) : dataKey );
+#endif // __GNUC__ || __GNUC__ < 4
+					AddAtomic( _pointWeightSum , sample.weight );
+				}
 			}
 		}
 	}
 	);
 	pointWeightSum = _pointWeightSum / weightSum;
 	MemoryUsage();
-	return normalField;
+	return dataField;
 }
 template< unsigned int Dim , class Real >
 template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
@@ -412,7 +427,7 @@ SparseNodeData< Data , IsotropicUIntPack< Dim , DataSig > > FEMTree< Dim , Real
 }
 template< unsigned int Dim , class Real >
 template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
-SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > FEMTree< Dim , Real >::setDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest )
+SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > FEMTree< Dim , Real >::setMultiDepthDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest )
 {
 	Allocator< FEMTreeNode > *nodeAllocator = nodeAllocators.size() ? nodeAllocators[0] : NULL;
 	LocalDepth maxDepth = _spaceRoot->maxDepth();
diff --git a/Src/Image.h b/Src/Image.h
index 563ff47..067147f 100644
--- a/Src/Image.h
+++ b/Src/Image.h
@@ -22,7 +22,7 @@ struct ImageReader
 	{
 		unsigned int channels;
 		ImageReader* reader = Get( fileName );
-		width = reader->width() , height = reader->height();
+		width = reader->width() , height = reader->height() , channels = reader->channels();
 		if( channels!=1 && channels!=3 ) ERROR_OUT( "Requres one- or three-channel input" );
 		unsigned char* pixels = new unsigned char[ width*height*3 ];
 		unsigned char* pixelRow = new unsigned char[ width*channels];
diff --git a/Src/Ply.h b/Src/Ply.h
index 1255977..c04fa1d 100644
--- a/Src/Ply.h
+++ b/Src/Ply.h
@@ -30,8 +30,8 @@ WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
 
 */
 
-#ifndef __PLY_H__
-#define __PLY_H__
+#ifndef PLY_INCLUDED
+#define PLY_INCLUDED
 
 #include <limits>
 #include <vector>
@@ -569,4 +569,4 @@ int PlyWritePolygons( const char* fileName , CoredMeshData< Vertex , Index >* me
 }
 inline int PlyDefaultFileType( void ){ return PLY_ASCII; }
 
-#endif /* !__PLY_H__ */
+#endif // PLY_INCLUDED
diff --git a/Src/PlyFile.h b/Src/PlyFile.h
index c851c65..42f8c35 100644
--- a/Src/PlyFile.h
+++ b/Src/PlyFile.h
@@ -29,8 +29,8 @@ EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
 WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.   
 */
 
-#ifndef __PLY_FILE_H__
-#define __PLY_FILE_H__
+#ifndef PLY_FILE_INCLUDED
+#define PLY_FILE_INCLUDED
 
 #include <string>
 #include <vector>
@@ -198,4 +198,4 @@ protected:
 };
 
 #include "PlyFile.inl"
-#endif /* !__PLY_FILE_H__ */
+#endif // PLY_FILE_INCLUDED
diff --git a/Src/PointInterpolant.cpp b/Src/PointInterpolant.cpp
new file mode 100644
index 0000000..bece832
--- /dev/null
+++ b/Src/PointInterpolant.cpp
@@ -0,0 +1,820 @@
+/*
+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 "PreProcessor.h"
+
+#undef USE_DOUBLE								// If enabled, double-precesion is used
+#define WEIGHT_DEGREE 2							// The order of the B-Spline used to splat in the weights for density estimation
+#define DEFAULT_FEM_DEGREE 2					// The default finite-element degree
+#define DEFAULT_FEM_BOUNDARY BOUNDARY_FREE		// The default finite-element boundary type
+#define DEFAULT_DIMENSION 2						// The dimension of the system
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <float.h>
+#include "MyMiscellany.h"
+#include "CmdLineParser.h"
+#include "PPolynomial.h"
+#include "FEMTree.h"
+#include "Ply.h"
+#include "PointStreamData.h"
+#include "Image.h"
+
+MessageWriter messageWriter;
+
+cmdLineParameter< char* >
+	In( "in" ) ,
+	Out( "out" ) ,
+	TempDir( "tempDir" ) ,
+	Grid( "grid" ) ,	
+	Tree( "tree" ) ,
+	Transform( "xForm" );
+
+cmdLineReadable
+	Performance( "performance" ) ,
+	ShowResidual( "showResidual" ) ,
+	PrimalGrid( "primalGrid" ) ,
+	ExactInterpolation( "exact" ) ,
+	InCore( "inCore" ) ,
+	NoValueConstraints( "noValues" ) ,
+	UseGradientConstraints( "useGradients" ) ,
+	NoComments( "noComments" ) ,
+	PolygonMesh( "polygonMesh" ) ,
+	NonManifold( "nonManifold" ) ,
+	NonLinearFit( "nonLinearFit" ) ,
+	ASCII( "ascii" ) ,
+	Verbose( "verbose" );
+
+cmdLineParameter< int >
+#ifndef FAST_COMPILE
+	Degree( "degree" , DEFAULT_FEM_DEGREE ) ,
+#endif // !FAST_COMPILE
+	Depth( "depth" , 8 ) ,
+	Iters( "iters" , 8 ) ,
+	FullDepth( "fullDepth" , 5 ) ,
+	BaseDepth( "baseDepth" , 5 ) ,
+	BaseVCycles( "baseVCycles" , 4 ) ,
+#ifndef FAST_COMPILE
+	BType( "bType" , DEFAULT_FEM_BOUNDARY+1 ) ,
+	Dimension( "dim" , DEFAULT_DIMENSION ) ,
+#endif // !FAST_COMPILE
+	MaxMemoryGB( "maxMemory" , 0 ) ,
+	ParallelType( "parallel" , (int)ThreadPool::OPEN_MP ) ,
+	ScheduleType( "schedule" , (int)ThreadPool::DefaultSchedule ) ,
+	ThreadChunkSize( "chunkSize" , (int)ThreadPool::DefaultChunkSize ) ,
+	Threads( "threads" , (int)std::thread::hardware_concurrency() );
+
+cmdLineParameter< float >
+	Scale( "scale" , 1.1f ) ,
+	Width( "width" , 0.f ) ,
+	CGSolverAccuracy( "cgAccuracy" , 1e-3f ) ,
+	IsoValue( "iso" , 0.f ) ,
+	ValueWeight   (    "valueWeight" , 1000.f ) ,
+	GradientWeight( "gradientWeight" , 1.f ) ,
+	LapWeight     (      "lapWeight" , 0.f ) ,
+	BiLapWeight   (    "biLapWeight" , 1.f );
+
+cmdLineReadable* params[] =
+{
+#ifndef FAST_COMPILE
+	&Degree , &BType , &Dimension ,
+#endif // !FAST_COMPILE
+	&In , &Out , &Depth , &Transform ,
+	&Width ,
+	&Scale , &Verbose , &CGSolverAccuracy , &NoComments ,
+	&NonManifold , &PolygonMesh , &ASCII , &ShowResidual ,
+	&ValueWeight , &GradientWeight ,
+	&LapWeight , &BiLapWeight ,
+	&Grid , &Threads ,
+	&Tree ,
+	&FullDepth ,
+	&BaseDepth , &BaseVCycles ,
+	&Iters ,
+	&IsoValue ,
+	&PrimalGrid ,
+	&ExactInterpolation ,
+	&Performance ,
+	&MaxMemoryGB ,
+	&InCore ,
+	&ParallelType ,
+	&ScheduleType ,
+	&ThreadChunkSize ,
+	&NoValueConstraints ,
+	&UseGradientConstraints ,
+	&NonLinearFit ,
+	NULL
+};
+
+void ShowUsage(char* ex)
+{
+	printf( "Usage: %s\n" , ex );
+	printf( "\t --%s <input points>\n" , In.name );
+	printf( "\t[--%s <ouput mesh>]\n" , Out.name );
+	printf( "\t[--%s <ouput grid>]\n" , Grid.name );
+	printf( "\t[--%s <ouput fem tree>]\n" , Tree.name );
+#ifndef FAST_COMPILE
+	printf( "\t[--%s <dimension>=%d]\n" , Dimension.name , Dimension.value );
+	printf( "\t[--%s <b-spline degree>=%d]\n" , Degree.name , Degree.value );
+	printf( "\t[--%s <boundary type>=%d]\n" , BType.name , BType.value );
+	for( int i=0 ; i<BOUNDARY_COUNT ; i++ ) printf( "\t\t%d] %s\n" , i+1 , BoundaryNames[i] );
+#endif // !FAST_COMPILE
+	printf( "\t[--%s <maximum reconstruction depth>=%d]\n" , Depth.name , Depth.value );
+	printf( "\t[--%s <grid width>]\n" , Width.name );
+	printf( "\t[--%s <full depth>=%d]\n" , FullDepth.name , FullDepth.value );
+	printf( "\t[--%s <coarse MG solver depth>=%d]\n" , BaseDepth.name , BaseDepth.value );
+	printf( "\t[--%s <coarse MG solver v-cycles>=%d]\n" , BaseVCycles.name , BaseVCycles.value );
+	printf( "\t[--%s <scale factor>=%f]\n" , Scale.name , Scale.value );
+	printf( "\t[--%s <zero-crossing weight>=%.3e]\n" , ValueWeight.name , ValueWeight.value );
+	printf( "\t[--%s <gradient weight>=%.3e]\n" , GradientWeight.name , GradientWeight.value );
+	printf( "\t[--%s <laplacian weight>=%.3e]\n" , LapWeight.name , LapWeight.value );
+	printf( "\t[--%s <bi-laplacian weight>=%.3e]\n" , BiLapWeight.name , BiLapWeight.value );
+	printf( "\t[--%s <iterations>=%d]\n" , Iters.name , Iters.value );
+	printf( "\t[--%s]\n" , ExactInterpolation.name );
+	printf( "\t[--%s <num threads>=%d]\n" , Threads.name , Threads.value );
+	printf( "\t[--%s <parallel type>=%d]\n" , ParallelType.name , ParallelType.value );
+	for( size_t i=0 ; i<ThreadPool::ParallelNames.size() ; i++ ) printf( "\t\t%d] %s\n" , (int)i , ThreadPool::ParallelNames[i].c_str() );
+	printf( "\t[--%s <schedue type>=%d]\n" , ScheduleType.name , ScheduleType.value );
+	for( size_t i=0 ; i<ThreadPool::ScheduleNames.size() ; i++ ) printf( "\t\t%d] %s\n" , (int)i , ThreadPool::ScheduleNames[i].c_str() );
+	printf( "\t[--%s <thread chunk size>=%d]\n" , ThreadChunkSize.name , ThreadChunkSize.value );
+	printf( "\t[--%s <cg solver accuracy>=%g]\n" , CGSolverAccuracy.name , CGSolverAccuracy.value );
+	printf( "\t[--%s <maximum memory (in GB)>=%d]\n" , MaxMemoryGB.name , MaxMemoryGB.value );
+	printf( "\t[--%s <iso-value>=%f]\n" , IsoValue.name , IsoValue.value );
+	printf( "\t[--%s]\n" , NoValueConstraints.name );
+	printf( "\t[--%s]\n" , UseGradientConstraints.name );
+	printf( "\t[--%s]\n" , Performance.name );
+	printf( "\t[--%s]\n" , PrimalGrid.name );
+	printf( "\t[--%s]\n" , NoComments.name );
+	printf( "\t[--%s]\n" , PolygonMesh.name );
+	printf( "\t[--%s]\n" , NonManifold.name );
+	printf( "\t[--%s]\n" , NonLinearFit.name );
+	printf( "\t[--%s]\n" , ASCII.name );
+	printf( "\t[--%s]\n" , InCore.name );
+	printf( "\t[--%s]\n" , Verbose.name );
+}
+
+template< unsigned int Dim , class Real >
+struct FEMTreeProfiler
+{
+	FEMTree< Dim , Real >& tree;
+	double t;
+
+	FEMTreeProfiler( FEMTree< Dim , Real >& t ) : tree(t) { ; }
+	void start( void ){ t = Time() , FEMTree< Dim , Real >::ResetLocalMemoryUsage(); }
+	void print( const char* header ) const
+	{
+		FEMTree< Dim , Real >::MemoryUsage();
+		if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+		else         printf(    "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" ,          Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+	}
+	void dumpOutput( const char* header ) const
+	{
+		FEMTree< Dim , Real >::MemoryUsage();
+		if( header ) messageWriter( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+		else         messageWriter(    "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" ,          Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+	}
+	void dumpOutput2( std::vector< std::string >& comments , const char* header ) const
+	{
+		FEMTree< Dim , Real >::MemoryUsage();
+		if( header ) messageWriter( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+		else         messageWriter( comments ,    "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" ,          Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+	}
+};
+
+template< class Real , unsigned int Dim >
+XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real scaleFactor )
+{
+	Point< Real , Dim > center = ( max + min ) / 2;
+	Real scale = max[0] - min[0];
+	for( int d=1 ; d<Dim ; d++ ) scale = std::max< Real >( scale , max[d]-min[d] );
+	scale *= scaleFactor;
+	for( int i=0 ; i<Dim ; i++ ) center[i] -= scale/2;
+	XForm< Real , Dim+1 > tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity();
+	for( int i=0 ; i<Dim ; i++ ) sXForm(i,i) = (Real)(1./scale ) , tXForm(Dim,i) = -center[i];
+	return sXForm * tXForm;
+}
+template< class Real , unsigned int Dim >
+XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real width , Real scaleFactor , int& depth )
+{
+	// Get the target resolution (along the largest dimension)
+	Real resolution = ( max[0]-min[0] ) / width;
+	for( int d=1 ; d<Dim ; d++ ) resolution = std::max< Real >( resolution , ( max[d]-min[d] ) / width );
+	resolution *= scaleFactor;
+	depth = 0;
+	while( (1<<depth)<resolution ) depth++;
+
+	Point< Real , Dim > center = ( max + min ) / 2;
+	Real scale = (1<<depth) * width;
+
+	for( int i=0 ; i<Dim ; i++ ) center[i] -= scale/2;
+	XForm< Real , Dim+1 > tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity();
+	for( int i=0 ; i<Dim ; i++ ) sXForm(i,i) = (Real)(1./scale ) , tXForm(Dim,i) = -center[i];
+	return sXForm * tXForm;
+}
+
+template< class Real , unsigned int Dim >
+XForm< Real , Dim+1 > GetPointXForm( InputPointStream< Real , Dim >& stream , Real width , Real scaleFactor , int& depth )
+{
+	Point< Real , Dim > min , max;
+	stream.boundingBox( min , max );
+	return GetBoundingBoxXForm( min , max , width , scaleFactor , depth );
+}
+template< class Real , unsigned int Dim >
+XForm< Real , Dim+1 > GetPointXForm( InputPointStream< Real , Dim >& stream , Real scaleFactor )
+{
+	Point< Real , Dim > min , max;
+	stream.boundingBox( min , max );
+	return GetBoundingBoxXForm( min , max , scaleFactor );
+}
+
+template< unsigned int Dim , typename Real , typename TotalPointSampleData > struct ValueAndGradientFromSample;
+
+template< unsigned int Dim , typename Real >
+struct ValueAndGradientFromSample< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > >
+{
+	typedef MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+	std::pair< Real , Point< Real , Dim > > operator()( TotalPointSampleData d ) const { return std::pair< Real , Point< Real , Dim > >( d.template data<0>() , d.template data<1>() ); }
+};
+
+template< unsigned int Dim , typename Real >
+struct ValueAndGradientFromSample< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > > >
+{
+	typedef MultiPointStreamData< Real , PointStreamValue< Real > > TotalPointSampleData;
+	std::pair< Real , Point< Real , Dim > > operator()( TotalPointSampleData d ) const { return std::pair< Real , Point< Real , Dim > >( d.template data<0>() , Point< Real , Dim >() ); }
+};
+
+template< unsigned int Dim , typename Real >
+struct ValueAndGradientFromSample< Dim , Real , MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > >
+{
+	typedef MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+	std::pair< Real , Point< Real , Dim > > operator()( TotalPointSampleData d ) const { return std::pair< Real , Point< Real , Dim > >( (Real)0 , d.template data<0>() ); }
+};
+
+
+template< unsigned int Dim , typename Real , typename TotalPointSampleData > struct ConstraintDual;
+
+template< unsigned int Dim , typename Real >
+struct ConstraintDual< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > >
+{
+	typedef MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+	Real vWeight , gWeight;
+	ConstraintDual( Real v , Real g ) : vWeight(v) , gWeight(g) { }
+	CumulativeDerivativeValues< Real , Dim , 1 > operator()( const Point< Real , Dim >& p , const TotalPointSampleData& data ) const 
+	{
+		Real value = data.template data<0>();
+		Point< Real , Dim > gradient = data.template data<1>();
+		CumulativeDerivativeValues< Real , Dim , 1 > cdv;
+		cdv[0] = value*vWeight;
+		for( int d=0 ; d<Dim ; d++ ) cdv[1+d] = gradient[d]*gWeight;
+		return cdv;
+	}
+};
+template< unsigned int Dim , typename Real >
+struct ConstraintDual< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > > >
+{
+	typedef MultiPointStreamData< Real , PointStreamValue< Real > > TotalPointSampleData;
+	Real vWeight , gWeight;
+	ConstraintDual( Real v , Real g ) : vWeight(v) , gWeight(g) { }
+	CumulativeDerivativeValues< Real , Dim , 1 > operator()( const Point< Real , Dim >& p , const TotalPointSampleData& data ) const 
+	{
+		Real value = data.template data<0>();
+		CumulativeDerivativeValues< Real , Dim , 1 > cdv;
+		cdv[0] = value*vWeight;
+		return cdv;
+	}
+};
+template< unsigned int Dim , typename Real >
+struct ConstraintDual< Dim , Real , MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > >
+{
+	typedef MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+	Real vWeight , gWeight;
+	ConstraintDual( Real v , Real g ) : vWeight(v) , gWeight(g) { }
+	CumulativeDerivativeValues< Real , Dim , 1 > operator()( const Point< Real , Dim >& p , const TotalPointSampleData& data ) const 
+	{
+		Point< Real , Dim > gradient = data.template data<0>();
+		CumulativeDerivativeValues< Real , Dim , 1 > cdv;
+		for( int d=0 ; d<Dim ; d++ ) cdv[1+d] = gradient[d]*gWeight;
+		return cdv;
+	}
+};
+template< unsigned int Dim , typename Real , typename TotalPointSampleData >
+struct SystemDual
+{
+	CumulativeDerivativeValues< Real , Dim , 1 > weight;
+	SystemDual( Real v , Real g )
+	{
+		weight[0] = v;
+		for( int d=0 ; d<Dim ; d++ ) weight[d+1] = g;
+	}
+	CumulativeDerivativeValues< Real , Dim , 1 > operator()( Point< Real , Dim > p , const TotalPointSampleData& data , const CumulativeDerivativeValues< Real , Dim , 1 >& dValues ) const
+	{
+		return dValues * weight;
+	}
+	CumulativeDerivativeValues< double , Dim , 1 > operator()( Point< Real , Dim > p , const TotalPointSampleData& data , const CumulativeDerivativeValues< double , Dim , 1 >& dValues ) const
+	{
+		return dValues * weight;
+	};
+};
+template< unsigned int Dim , class TotalPointSampleData >
+struct SystemDual< Dim , double , TotalPointSampleData >
+{
+	typedef double Real;
+	CumulativeDerivativeValues< Real , Dim , 1 > weight;
+	SystemDual( Real v , Real g ) : weight( v , g , g , g ) { }
+	CumulativeDerivativeValues< Real , Dim , 1 > operator()( Point< Real , Dim > p , const TotalPointSampleData& data , const CumulativeDerivativeValues< Real , Dim , 1 >& dValues ) const
+	{
+		return dValues * weight;
+	}
+};
+
+template< typename Vertex , typename Real , unsigned int ... FEMSigs , typename TotalPointSampleData >
+void ExtractMesh( UIntPack< FEMSigs ... > , FEMTree< sizeof ... ( FEMSigs ) , Real >& tree , const DenseNodeData< Real , UIntPack< FEMSigs ... > >& solution , Real isoValue , const std::vector< typename FEMTree< sizeof ... ( FEMSigs ) , Real >::PointSample >* samples , std::function< void ( Vertex& , Point< Real , sizeof ... ( FEMSigs ) > , Real , TotalPointSampleData ) > SetVertex , std::vector< std::string > &comments , XForm< Real , sizeof...(FEMSigs)+1 > iXForm )
+{
+	static const int Dim = sizeof ... ( FEMSigs );
+	typedef UIntPack< FEMSigs ... > Sigs;
+	static const unsigned int DataSig = FEMDegreeAndBType< WEIGHT_DEGREE , BOUNDARY_FREE >::Signature;
+
+	FEMTreeProfiler< Dim , Real > profiler( tree );
+
+	char tempHeader[1024];
+	{
+		char tempPath[1024];
+		tempPath[0] = 0;
+		if( TempDir.set ) strcpy( tempPath , TempDir.value );
+		else SetTempDirectory( tempPath , sizeof(tempPath) );
+		if( strlen(tempPath)==0 ) sprintf( tempPath , ".%c" , FileSeparator );
+		if( tempPath[ strlen( tempPath )-1 ]==FileSeparator ) sprintf( tempHeader , "%sPR_" , tempPath );
+		else                                                  sprintf( tempHeader , "%s%cPR_" , tempPath , FileSeparator );
+	}
+
+	CoredMeshData< Vertex , node_index_type > *mesh;
+	if( InCore.set ) mesh = new CoredVectorMeshData< Vertex , node_index_type >();
+	else             mesh = new CoredFileMeshData< Vertex , node_index_type >( tempHeader );
+	profiler.start();
+	typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats;
+#if defined( __GNUC__ ) && __GNUC__ < 5
+#warning "you've got me gcc version<5"
+	isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< TotalPointSampleData >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , (typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE >*)NULL , (SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > *)NULL , solution , isoValue , *mesh , SetVertex , NonLinearFit.set , !NonManifold.set , PolygonMesh.set , false );
+#else // !__GNUC__ || __GNUC__ >=5
+	isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< TotalPointSampleData >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , (typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE >*)NULL , NULL , solution , isoValue , *mesh , SetVertex , NonLinearFit.set , !NonManifold.set , PolygonMesh.set , false );
+#endif // __GNUC__ || __GNUC__ < 4
+	messageWriter( "Vertices / Polygons: %llu / %llu\n" , (unsigned long long)( mesh->outOfCorePointCount()+mesh->inCorePoints.size() ) , (unsigned long long)mesh->polygonCount() );
+	std::string isoStatsString = isoStats.toString() + std::string( "\n" );
+	messageWriter( isoStatsString.c_str() );
+	if( PolygonMesh.set ) profiler.dumpOutput2( comments , "#         Got polygons:" );
+	else                  profiler.dumpOutput2( comments , "#        Got triangles:" );
+
+	std::vector< std::string > noComments;
+	if( !PlyWritePolygons< Vertex , node_index_type , Real , Dim >( Out.value , mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , NoComments.set ? noComments : comments , iXForm ) )
+		ERROR_OUT( "Could not write mesh to: " , Out.value );
+	delete mesh;
+}
+
+template< typename Real , unsigned int Dim >
+void WriteGrid( ConstPointer( Real ) values , int res , const char *fileName )
+{
+	int resolution = 1;
+	for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+
+	char *ext = GetFileExtension( fileName );
+
+	if( Dim==2 && ImageWriter::ValidExtension( ext ) )
+	{
+		Real avg = 0;
+		std::vector< Real > avgs( ThreadPool::NumThreads() , 0 );
+		ThreadPool::Parallel_for( 0 , resolution , [&]( unsigned int thread , size_t i ){ avgs[thread] += values[i]; } );
+		for( unsigned int t=0 ; t<ThreadPool::NumThreads() ; t++ ) avg += avgs[t];
+		avg /= (Real)resolution;
+
+		Real std = 0;
+		std::vector< Real > stds( ThreadPool::NumThreads() , 0 );
+		ThreadPool::Parallel_for( 0 , resolution , [&]( unsigned int thread , size_t i ){ stds[thread] += ( values[i] - avg ) * ( values[i] - avg ); } );
+		for( unsigned int t=0 ; t<ThreadPool::NumThreads() ; t++ ) std += stds[t];
+		std = (Real)sqrt( std / resolution );
+
+		if( Verbose.set ) printf( "Grid to image: [%.2f,%.2f] -> [0,255]\n" , avg - 2*std , avg + 2*std );
+
+		unsigned char *pixels = new unsigned char[ resolution*3 ];
+		ThreadPool::Parallel_for( 0 , resolution , [&]( unsigned int , size_t i )
+		{
+			Real v = (Real)std::min< Real >( (Real)1. , std::max< Real >( (Real)-1. , ( values[i] - avg ) / (2*std ) ) );
+			v = (Real)( ( v + 1. ) / 2. * 256. );
+			unsigned char color = (unsigned char )std::min< Real >( (Real)255. , std::max< Real >( (Real)0. , v ) );
+			for( int c=0 ; c<3 ; c++ ) pixels[i*3+c ] = color;
+		}
+		);
+		ImageWriter::Write( fileName , pixels , res , res , 3 );
+		delete[] pixels;
+	}
+	else
+	{
+
+		FILE *fp = fopen( fileName , "wb" );
+		if( !fp ) ERROR_OUT( "Failed to open grid file for writing: " , fileName );
+		else
+		{
+			fwrite( &res , sizeof(int) , 1 , fp );
+			if( typeid(Real)==typeid(float) ) fwrite( values , sizeof(float) , resolution , fp );
+			else
+			{
+				float *fValues = new float[resolution];
+				for( int i=0 ; i<resolution ; i++ ) fValues[i] = float( values[i] );
+				fwrite( fValues , sizeof(float) , resolution , fp );
+				delete[] fValues;
+			}
+			fclose( fp );
+		}
+	}
+	delete[] ext;
+}
+
+
+template< class Real , typename TotalPointSampleData , unsigned int ... FEMSigs >
+void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
+{
+	static const int Dim = sizeof ... ( FEMSigs );
+	typedef UIntPack< FEMSigs ... > Sigs;
+	typedef UIntPack< FEMSignature< FEMSigs >::Degree ... > Degrees;
+	typedef UIntPack< FEMDegreeAndBType< WEIGHT_DEGREE , DerivativeBoundary< FEMSignature< FEMSigs >::BType , 1 >::BType >::Signature ... > DataSigs;
+	typedef typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator;
+	typedef typename FEMTree< Dim , Real >::template InterpolationInfo< Real , 1 > InterpolationInfo;
+	typedef InputPointStreamWithData< Real , Dim , TotalPointSampleData > InputPointStream;
+	typedef TransformedInputPointStreamWithData< Real , Dim , TotalPointSampleData > XInputPointStream;
+	std::vector< std::string > comments;
+	messageWriter( comments , "***********************************************\n" );
+	messageWriter( comments , "***********************************************\n" );
+	messageWriter( comments , "** Running Point Interpolant (Version %s) **\n" , VERSION );
+	messageWriter( comments , "***********************************************\n" );
+	messageWriter( comments , "***********************************************\n" );
+	if( !Threads.set ) messageWriter( comments , "Running with %d threads\n" , Threads.value );
+
+	ThreadPool::Init( (ThreadPool::ParallelType)ParallelType.value , Threads.value );
+
+	XForm< Real , Dim+1 > xForm , iXForm;
+	if( Transform.set )
+	{
+		FILE* fp = fopen( Transform.value , "r" );
+		if( !fp )
+		{
+			WARN( "Could not read x-form from: " , Transform.value );
+			xForm = XForm< Real , Dim+1 >::Identity();
+		}
+		else
+		{
+			for( int i=0 ; i<Dim+1 ; i++ ) for( int j=0 ; j<Dim+1 ; j++ )
+			{
+				float f;
+				if( fscanf( fp , " %f " , &f )!=1 ) ERROR_OUT( "Failed to read xform" );
+				xForm(i,j) = (Real)f;
+			}
+			fclose( fp );
+		}
+	}
+	else xForm = XForm< Real , Dim+1 >::Identity();
+
+	char str[1024];
+	for( int i=0 ; params[i] ; i++ )
+		if( params[i]->set )
+		{
+			params[i]->writeValue( str );
+			if( strlen( str ) ) messageWriter( comments , "\t--%s %s\n" , params[i]->name , str );
+			else                messageWriter( comments , "\t--%s\n" , params[i]->name );
+		}
+
+	double startTime = Time();
+
+	FEMTree< Dim , Real > tree( MEMORY_ALLOCATOR_BLOCK_SIZE );
+	FEMTreeProfiler< Dim , Real > profiler( tree );
+
+	if( Depth.set && Width.value>0 )
+	{
+		WARN( "Both --" , Depth.name , " and --" , Width.name , " set, ignoring --" , Width.name );
+		Width.value = 0;
+	}
+
+	size_t pointCount;
+
+	std::vector< typename FEMTree< Dim , Real >::PointSample >* samples = new std::vector< typename FEMTree< Dim , Real >::PointSample >();
+	std::vector< TotalPointSampleData >* sampleData = NULL;
+
+	// Read in the samples
+	{
+		profiler.start();
+		InputPointStream* pointStream;
+		char* ext = GetFileExtension( In.value );
+		sampleData = new std::vector< TotalPointSampleData >();
+		std::vector< std::pair< Point< Real , Dim > , TotalPointSampleData > > inCorePoints;
+		if( InCore.set )
+		{
+			InputPointStream *_pointStream;
+			if     ( !strcasecmp( ext , "bnpts" ) ) _pointStream = new BinaryInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadBinary );
+			else if( !strcasecmp( ext , "ply"   ) ) _pointStream = new    PLYInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::PlyReadProperties() , TotalPointSampleData::PlyReadNum , TotalPointSampleData::ValidPlyReadProperties );
+			else                                    _pointStream = new  ASCIIInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadASCII );
+			Point< Real , Dim > p;
+			TotalPointSampleData d;
+			while( _pointStream->nextPoint( p , d ) ) inCorePoints.push_back( std::pair< Point< Real , Dim > , TotalPointSampleData >( p , d ) );
+			delete _pointStream;
+
+			pointStream = new MemoryInputPointStreamWithData< Real , Dim , TotalPointSampleData >( inCorePoints.size() , &inCorePoints[0] );
+		}
+		else
+		{
+			if     ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadBinary );
+			else if( !strcasecmp( ext , "ply"   ) ) pointStream = new    PLYInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::PlyReadProperties() , TotalPointSampleData::PlyReadNum , TotalPointSampleData::ValidPlyReadProperties );
+			else                                    pointStream = new  ASCIIInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadASCII );
+		}
+		delete[] ext;
+		typename TotalPointSampleData::Transform _xForm( xForm );
+		XInputPointStream _pointStream( [&]( Point< Real , Dim >& p , TotalPointSampleData& d ){ p = xForm*p , d = _xForm(d); } , *pointStream );
+		if( Width.value>0 ) xForm = GetPointXForm< Real , Dim >( _pointStream , Width.value , (Real)( Scale.value>0 ? Scale.value : 1. ) , Depth.value ) * xForm;
+		else                xForm = Scale.value>0 ? GetPointXForm< Real , Dim >( _pointStream , (Real)Scale.value ) * xForm : xForm;
+		{
+			typename TotalPointSampleData::Transform _xForm( xForm );
+			XInputPointStream _pointStream( [&]( Point< Real , Dim >& p , TotalPointSampleData& d ){ p = xForm*p , d = _xForm(d); } , *pointStream );
+			auto ProcessData = []( const Point< Real , Dim >& p , TotalPointSampleData& d ){ return (Real)1.; };
+			pointCount = FEMTreeInitializer< Dim , Real >::template Initialize< TotalPointSampleData >( tree.spaceRoot() , _pointStream , Depth.value , *samples , *sampleData , true , tree.nodeAllocators.size() ? tree.nodeAllocators[0] : NULL , tree.initializer() , ProcessData );
+		}
+		iXForm = xForm.inverse();
+		delete pointStream;
+
+		messageWriter( "Input Points / Samples: %llu / %llu\n" , pointCount , (unsigned long long)samples->size() );
+		profiler.dumpOutput2( comments , "# Read input into tree:" );
+	}
+
+	DenseNodeData< Real , Sigs > solution;
+	{
+		DenseNodeData< Real , Sigs > constraints;
+		InterpolationInfo* iInfo = NULL;
+		int solveDepth = Depth.value;
+
+		tree.resetNodeIndices();
+
+		// Prepare for multigrid
+		{
+			profiler.start();
+			tree.template finalizeForMultigrid< Degrees::Max() >( FullDepth.value , []( const typename FEMTree< Dim , Real >::FEMTreeNode * ){ return true; } );
+			profiler.dumpOutput2( comments , "#       Finalized tree:" );
+		}
+
+		// Add the interpolation constraints
+		{
+			profiler.start();
+			if( ExactInterpolation.set ) iInfo = FEMTree< Dim , Real >::template       InitializeExactPointAndDataInterpolationInfo< Real , TotalPointSampleData , 1 >( tree , *samples , GetPointer( *sampleData ) , ConstraintDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , SystemDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , true , false );
+			else                         iInfo = FEMTree< Dim , Real >::template InitializeApproximatePointAndDataInterpolationInfo< Real , TotalPointSampleData , 1 >( tree , *samples , GetPointer( *sampleData ) , ConstraintDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , SystemDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , true , 1 );
+			constraints = tree.initDenseNodeData( Sigs() );
+			tree.addInterpolationConstraints( constraints , solveDepth , *iInfo );
+			profiler.dumpOutput2( comments , "#Set point constraints:" );
+		}
+
+		messageWriter( "Leaf Nodes / Active Nodes / Ghost Nodes: %llu / %llu / %llu\n" , (unsigned long long)tree.leaves() , (unsigned long long)tree.nodes() , (unsigned long long)tree.ghostNodes() );
+		messageWriter( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) );
+
+		// Solve the linear system
+		{
+			profiler.start();
+			typename FEMTree< Dim , Real >::SolverInfo sInfo;
+			sInfo.cgDepth = 0 , sInfo.cascadic = true , sInfo.vCycles = 1 , sInfo.iters = Iters.value , sInfo.cgAccuracy = CGSolverAccuracy.value , sInfo.verbose = Verbose.set , sInfo.showResidual = ShowResidual.set , sInfo.showGlobalResidual = SHOW_GLOBAL_RESIDUAL_NONE , sInfo.sliceBlockSize = 1;
+			sInfo.baseDepth = BaseDepth.value , sInfo.baseVCycles = BaseVCycles.value;
+			typename FEMIntegrator::template System< Sigs , IsotropicUIntPack< Dim , 2 > > F( { 0. , (double)LapWeight.value , (double)BiLapWeight.value } );
+			solution = tree.solveSystem( Sigs() , F , constraints , solveDepth , sInfo , iInfo );
+			profiler.dumpOutput2( comments , "# Linear system solved:" );
+			if( iInfo ) delete iInfo , iInfo = NULL;
+		}
+	}
+
+	if( Verbose.set )
+	{
+		typename FEMTree< Dim , Real >::template MultiThreadedEvaluator< Sigs , 1 > evaluator( &tree , solution );
+		std::pair< double , double > valueStat(0,0) , gradientStat(0,0);
+		std::vector< std::pair< double , double > > valueStats( ThreadPool::NumThreads() , std::pair< double , double >(0,0) ) , gradientStats( ThreadPool::NumThreads() , std::pair< double , double >(0,0) );
+		ValueAndGradientFromSample< Dim , Real , TotalPointSampleData > valueAndGradientFromSample;
+		ThreadPool::Parallel_for( 0 , samples->size() , [&]( unsigned int thread , size_t j )
+		{
+			ProjectiveData< Point< Real , Dim > , Real >& sample = (*samples)[j].sample;
+			Real w = sample.weight;
+			if( w>0 )
+			{
+				CumulativeDerivativeValues< Real , Dim , 1 > values = evaluator.values( sample.data / sample.weight , thread , (*samples)[j].node );
+				Real value = values[0];
+				Point< Real , Dim > gradient;
+				for( int d=0 ; d<Dim ; d++ ) gradient[d] = values[d+1];
+				std::pair< Real , Point< Real , Dim > > valueAndGradient = valueAndGradientFromSample( (*sampleData)[j] / w );
+				valueStats[ thread ].first += ( value - valueAndGradient.first ) * ( value - valueAndGradient.first ) * w;
+				valueStats[ thread ].second += ( value * value + valueAndGradient.first * valueAndGradient.first ) * w;
+				gradientStats[ thread ].first += Point< Real , Dim >::SquareNorm( gradient - valueAndGradient.second ) * w;
+				gradientStats[ thread ].second += ( Point< Real , Dim >::SquareNorm( gradient ) + Point< Real , Dim >::SquareNorm( valueAndGradient.second ) ) * w;
+			}
+		}
+		);
+		for( unsigned int t=0 ; t<ThreadPool::NumThreads() ; t++ ) valueStat.first += valueStats[t].first , valueStat.second += valueStats[t].second , gradientStat.first += gradientStats[t].first , gradientStat.second += gradientStats[t].second;
+		if( ValueWeight.value>0 && GradientWeight.value>0 ) messageWriter( "Value / Gradient Error: %g / %g\n" , (Real)sqrt( valueStat.first / valueStat.second ) , (Real)sqrt( gradientStat.first / gradientStat.second ) );
+		else if( ValueWeight.value>0 ) messageWriter( "Value Error: %g\n" , (Real)sqrt( valueStat.first / valueStat.second ) );
+		else if( GradientWeight.value>0 ) messageWriter( "Gradient Error: %g\n" , (Real)sqrt( gradientStat.first / gradientStat.second ) );
+	}
+
+	delete samples , samples = NULL;
+	delete sampleData , sampleData = NULL;
+
+
+	if( Tree.set )
+	{
+		FILE* fp = fopen( Tree.value , "wb" );
+		if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Tree.value );
+		FEMTree< Dim , Real >::WriteParameter( fp );
+		DenseNodeData< Real , Sigs >::WriteSignatures( fp );
+		tree.write( fp , xForm );
+		solution.write( fp );
+		fclose( fp );
+	}
+
+	if( Grid.set )
+	{
+		int res = 0;
+		profiler.start();
+		Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalGrid.set );
+		size_t resolution = 1;
+		for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+		profiler.dumpOutput( "Got grid:" );
+		WriteGrid< Real , Dim >( values , res , Grid.value );
+		DeletePointer( values );
+		if( Verbose.set )
+		{
+			printf( "Transform:\n" );
+			for( int i=0 ; i<Dim+1 ; i++ )
+			{
+				printf( "\t" );
+				for( int j=0 ; j<Dim+1 ; j++ ) printf( " %f" , iXForm(j,i) );
+				printf( "\n" );
+			}
+		}
+	}
+
+	if( Out.set )
+	{
+		typedef PlyVertex< Real , Dim > Vertex;
+		std::function< void ( Vertex& , Point< Real , Dim > , Real , TotalPointSampleData ) > SetVertex = []( Vertex& v , Point< Real , Dim > p , Real , TotalPointSampleData ){ v.point = p; };
+		ExtractMesh< Vertex >( UIntPack< FEMSigs ... >() , tree , solution , IsoValue.value , samples , SetVertex , comments , iXForm );
+	}
+
+	messageWriter( comments , "#          Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-startTime , FEMTree< Dim , Real >::MaxMemoryUsage() );
+}
+
+template< class Real , unsigned int ... FEMSigs >
+void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
+{
+	static const int Dim = sizeof ... ( FEMSigs );
+	if     ( !UseGradientConstraints.set ) Execute< Real , MultiPointStreamData< Real , PointStreamValue< Real >                                   > >( argc , argv , UIntPack< FEMSigs ... >() );
+	else if( NoValueConstraints.set )      Execute< Real , MultiPointStreamData< Real ,                            PointStreamNormal< Real , Dim > > >( argc , argv , UIntPack< FEMSigs ... >() );
+	else                                   Execute< Real , MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > >( argc , argv , UIntPack< FEMSigs ... >() );
+}
+
+#ifndef FAST_COMPILE
+template< unsigned int Dim , class Real >
+void Execute( int argc , char* argv[] )
+{
+	switch( BType.value )
+	{
+	case BOUNDARY_FREE+1:
+	{
+		switch( Degree.value )
+		{
+//			case 1: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BOUNDARY_FREE >::Signature >() );
+			case 2: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BOUNDARY_FREE >::Signature >() );
+			case 3: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BOUNDARY_FREE >::Signature >() );
+//			case 4: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BOUNDARY_FREE >::Signature >() );
+			default: ERROR_OUT( "Only B-Splines of degree 1 - 3 are supported" );
+		}
+	}
+	case BOUNDARY_NEUMANN+1:
+	{
+		switch( Degree.value )
+		{
+//			case 1: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BOUNDARY_NEUMANN >::Signature >() );
+			case 2: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BOUNDARY_NEUMANN >::Signature >() );
+			case 3: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BOUNDARY_NEUMANN >::Signature >() );
+//			case 4: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BOUNDARY_NEUMANN >::Signature >() );
+			default: ERROR_OUT( "Only B-Splines of degree 1 - 3 are supported" );
+		}
+	}
+	case BOUNDARY_DIRICHLET+1:
+	{
+		switch( Degree.value )
+		{
+//			case 1: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BOUNDARY_DIRICHLET >::Signature >() );
+			case 2: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BOUNDARY_DIRICHLET >::Signature >() );
+			case 3: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BOUNDARY_DIRICHLET >::Signature >() );
+//			case 4: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BOUNDARY_DIRICHLET >::Signature >() );
+			default: ERROR_OUT( "Only B-Splines of degree 1 - 3 are supported" );
+		}
+	}
+	default: ERROR_OUT( "Not a valid boundary type: " , BType.value );
+	}
+}
+#endif // !FAST_COMPILE
+
+int main( int argc , char* argv[] )
+{
+	Timer timer;
+#ifdef USE_SEG_FAULT_HANDLER
+	WARN( "using seg-fault handler" );
+	StackTracer::exec = argv[0];
+	signal( SIGSEGV , SignalHandler );
+#endif // USE_SEG_FAULT_HANDLER
+#ifdef ARRAY_DEBUG
+	WARN( "Array debugging enabled" );
+#endif // ARRAY_DEBUG
+
+	cmdLineParse( argc-1 , &argv[1] , params );
+#if 0
+	if( !In.set || !Out.set ) ERROR_OUT( "Need input and output" );
+	unsigned int width , height;
+	unsigned char *pixels = ImageReader::ReadColor( In.value , width , height );
+	FILE *fp = fopen( Out.value , "wb" );
+	if( !fp ) ERROR_OUT( "Failed to open file for reading: %s" , Out.value );
+	for( int i=0 ; i<10000 ; i++ )
+	{
+		int x = rand() % width , y = rand() % height;
+		double gray = (double)( pixels[ 3*(y*width+x) + 0 ] + pixels[ 3*(y*width+x) + 1 ] + pixels[ 3*(y*width+x) + 2 ] ) / ( 255. * 3 );
+		fprintf( fp , "%d %d  %f\n" , x , y , gray );
+	}
+	fclose( fp );
+#else
+
+	if( MaxMemoryGB.value>0 ) SetPeakMemoryMB( MaxMemoryGB.value<<10 );
+	ThreadPool::DefaultChunkSize = ThreadChunkSize.value;
+	ThreadPool::DefaultSchedule = (ThreadPool::ScheduleType)ScheduleType.value;
+	messageWriter.echoSTDOUT = Verbose.set;
+
+	if( !In.set )
+	{
+		ShowUsage( argv[0] );
+		return 0;
+	}
+	if( NoValueConstraints.set ) ValueWeight.value = 0;
+	if( !UseGradientConstraints.set ) GradientWeight.value = 0;
+
+	if( ValueWeight.value<0 ) ERROR_OUT( "Value weight must be non-negative: " , ValueWeight.value , "> 0" );
+	if( GradientWeight.value<0 ) ERROR_OUT( "Gradient weight must be non-negative: " , GradientWeight.value , "> 0" );
+	if( !ValueWeight.value && !GradientWeight.value ) ERROR_OUT( "Either value or gradient weight must be positive" );
+
+	if( LapWeight.value<0 ) ERROR_OUT( "Laplacian weight must be non-negative: " , LapWeight.value , " > 0" );
+	if( BiLapWeight.value<0 ) ERROR_OUT( "Bi-Laplacian weight must be non-negative: " , BiLapWeight.value , " > 0" );
+	if( !LapWeight.value && !BiLapWeight.value ) ERROR_OUT( "Eiter Laplacian or bi-Laplacian weight must be positive" );
+
+	if( BaseDepth.value>FullDepth.value )
+	{
+		if( BaseDepth.set ) WARN( "Base depth must be smaller than full depth: " , BaseDepth.value , " <= " , FullDepth.value );
+		BaseDepth.value = FullDepth.value;
+	}
+
+#ifdef USE_DOUBLE
+	typedef double Real;
+#else // !USE_DOUBLE
+	typedef float  Real;
+#endif // USE_DOUBLE
+
+#ifdef FAST_COMPILE
+	static const int Dimension = DIMENSION;
+	static const int Degree = DEFAULT_FEM_DEGREE;
+	static const BoundaryType BType = DEFAULT_FEM_BOUNDARY;
+	typedef IsotropicUIntPack< Dimension , FEMDegreeAndBType< Degree , BType >::Signature > FEMSigs;
+	WARN( "Compiled for degree-" , Degree , ", boundary-" , BoundaryNames[ BType ] , ", " , sizeof(Real)==4 ? "single" : "double" , "-precision _only_" );
+	Execute< Real >( argc , argv , FEMSigs() );
+#else // !FAST_COMPILE
+	if( Dimension.value==2 ) Execute< 2 , Real >( argc , argv );
+	else if( Dimension.value==3 ) Execute< 3 , Real >( argc , argv );
+	else ERROR_OUT( "Only Degrees 2 and 3 are supported" );
+#endif // FAST_COMPILE
+	if( Performance.set )
+	{
+		printf( "Time (Wall/CPU): %.2f / %.2f\n" , timer.wallTime() , timer.cpuTime() );
+		printf( "Peak Memory (MB): %d\n" , MemoryInfo::PeakMemoryUsageMB() );
+	}
+#endif
+	return EXIT_SUCCESS;
+}
diff --git a/Src/PointStream.inl b/Src/PointStream.inl
index a14e130..5cafe61 100644
--- a/Src/PointStream.inl
+++ b/Src/PointStream.inl
@@ -150,8 +150,8 @@ void PLYInputPointStream< Real , Dim >::reset( void )
 		{
 			foundVertices = true;
 			_pCount = num_elems , _pIdx = 0;
-			for( int i=0 ; i<PlyVertex< Real , Dim >::ReadComponents ; i++ ) 
-				if( !_ply->get_property( elem_name , &(PlyVertex< Real , Dim >::Properties()[i]) ) ) ERROR_OUT( "Failed to find property in ply file: %s" , PlyVertex< Real , Dim >::Properties()[i].name );
+			for( int i=0 ; i<PlyVertex< Real , Dim >::PlyReadNum ; i++ ) 
+				if( !_ply->get_property( elem_name , &(PlyVertex< Real , Dim >::PlyReadProperties()[i]) ) ) ERROR_OUT( "Failed to find property in ply file: %s" , PlyVertex< Real , Dim >::PlyReadProperties()[i].name );
 		}
 		for( int j=0 ; j<plist.size() ; j++ ) delete plist[j];
 		if( foundVertices ) break;
diff --git a/Src/PointStreamData.h b/Src/PointStreamData.h
index 23f2235..b27528e 100644
--- a/Src/PointStreamData.h
+++ b/Src/PointStreamData.h
@@ -89,7 +89,7 @@ struct PointStreamPosition : public PointStreamData< Real , Point< Real , Dim >
 	protected:
 		XForm< Real , Dim+1 > _xForm;
 	};
-	static void readASCII( FILE* fp , PointStreamPosition& p )
+	static void ReadASCII( FILE* fp , PointStreamPosition& p )
 	{
 		float f;
 		for( int i=0 ; i<Dim ; i++ )
@@ -271,8 +271,8 @@ struct PointStreamColor : public PointStreamData< Real , Color< Real > >
 	};
 	static void ReadASCII( FILE* fp , PointStreamColor& p )
 	{
-		unsigned char c[3];
-		if( fscanf( fp , " %c %c %c " , &c[0] , &c[1] , &c[2] )!=3 ) ERROR_OUT( "Failed to read color" );
+		unsigned int c[3];
+		if( fscanf( fp , " %d %d %d " , &c[0] , &c[1] , &c[2] )!=3 ) ERROR_OUT( "Failed to read color" );
 		p.psData[0] = (Real)c[0] , p.psData[1] = (Real)c[1] , p.psData[2] = (Real)c[2];
 	};
 	static void ReadBinary( FILE* fp , PointStreamColor& p )
diff --git a/Src/PoissonRecon.cpp b/Src/PoissonRecon.cpp
index cb2e8e9..a5659ed 100644
--- a/Src/PoissonRecon.cpp
+++ b/Src/PoissonRecon.cpp
@@ -338,7 +338,7 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
 	typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats;
 	if( sampleData )
 	{
-		SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
+		SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setMultiDepthDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
 		for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) )
 		{
 			ProjectiveData< TotalPointSampleData , Real >* clr = _sampleData( n );
@@ -582,9 +582,27 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 		{
 			profiler.start();
 			normalInfo = new SparseNodeData< Point< Real , Dim > , NormalSigs >();
-			if( ConfidenceBias.value>0 ) *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , [&]( Real conf ){ return (Real)( log( conf ) * ConfidenceBias.value / log( 1<<(Dim-1) ) ); } );
-			else                         *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum );
-			ThreadPool::Parallel_for( 0 , normalInfo->size() , [&]( unsigned int , size_t i ){ (*normalInfo)[i] *= (Real)-1.; } );
+			std::function< bool ( TotalPointSampleData , Point< Real , Dim >& ) > ConversionFunction = []( TotalPointSampleData in , Point< Real , Dim > &out )
+			{
+				Point< Real , Dim > n = in.template data<0>();
+				Real l = (Real)Length( n );
+				// It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+				if( !l ) return false;
+				out = n / l;
+				return true;
+			};
+			std::function< bool ( TotalPointSampleData , Point< Real , Dim >& , Real & ) > ConversionAndBiasFunction = []( TotalPointSampleData in , Point< Real , Dim > &out , Real &bias )
+			{
+				Point< Real , Dim > n = in.template data<0>();
+				Real l = (Real)Length( n );
+				// It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+				if( !l ) return false;
+				out = n / l;
+				bias = (Real)( log( l ) * ConfidenceBias.value / log( 1<<(Dim-1) ) );
+				return true;
+			};
+			if( ConfidenceBias.value>0 ) *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionAndBiasFunction );
+			else                         *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionFunction );			ThreadPool::Parallel_for( 0 , normalInfo->size() , [&]( unsigned int , size_t i ){ (*normalInfo)[i] *= (Real)-1.; } );
 			profiler.dumpOutput2( comments , "#     Got normal field:" );
 			messageWriter( "Point weight / Estimated Area: %g / %g\n" , pointWeightSum , pointCount*pointWeightSum );
 		}
@@ -672,7 +690,7 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 		if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Tree.value );
 		FEMTree< Dim , Real >::WriteParameter( fp );
 		DenseNodeData< Real , Sigs >::WriteSignatures( fp );
-		tree.write( fp );
+		tree.write( fp , xForm );
 		solution.write( fp );
 		fclose( fp );
 	}
diff --git a/Src/PreProcessor.h b/Src/PreProcessor.h
index b7d76b9..c07380d 100644
--- a/Src/PreProcessor.h
+++ b/Src/PreProcessor.h
@@ -44,7 +44,7 @@ DAMAGE.
 #define USE_DEEP_TREE_NODES						// Chances are that if you are using big data, you want to support a tree with depth>15.
 #endif // BIG_DATA
 
-#define VERSION "11.02"							// The version of the code
+#define VERSION "12.00"							// The version of the code
 #define MEMORY_ALLOCATOR_BLOCK_SIZE 1<<12		// The chunk size for memory allocation
 
 #endif // PRE_PROCESSOR_INCLUDED
\ No newline at end of file
diff --git a/Src/SSDRecon.cpp b/Src/SSDRecon.cpp
index 3068f40..3a0b55c 100644
--- a/Src/SSDRecon.cpp
+++ b/Src/SSDRecon.cpp
@@ -358,7 +358,7 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
 	typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats;
 	if( sampleData )
 	{
-		SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
+		SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setMultiDepthDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
 		for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) )
 		{
 			ProjectiveData< TotalPointSampleData , Real >* clr = _sampleData( n );
@@ -603,8 +603,27 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 		{
 			profiler.start();
 			normalInfo = new SparseNodeData< Point< Real , Dim > , NormalSigs >();
-			if( ConfidenceBias.value>0 ) *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , [&]( Real conf ){ return (Real)( log( conf ) * ConfidenceBias.value / log( 1<<(Dim-1) ) ); } );
-			else                         *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum );
+			std::function< bool ( TotalPointSampleData , Point< Real , Dim >& ) > ConversionFunction = []( TotalPointSampleData in , Point< Real , Dim > &out )
+			{
+				Point< Real , Dim > n = in.template data<0>();
+				Real l = (Real)Length( n );
+				// It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+				if( !l ) return false;
+				out = n / l;
+				return true;
+			};
+			std::function< bool ( TotalPointSampleData , Point< Real , Dim >& , Real & ) > ConversionAndBiasFunction = []( TotalPointSampleData in , Point< Real , Dim > &out , Real &bias )
+			{
+				Point< Real , Dim > n = in.template data<0>();
+				Real l = (Real)Length( n );
+				// It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+				if( !l ) return false;
+				out = n / l;
+				bias = (Real)( log( l ) * ConfidenceBias.value / log( 1<<(Dim-1) ) );
+				return true;
+			};
+			if( ConfidenceBias.value>0 ) *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionAndBiasFunction );
+			else                         *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionFunction );
 			profiler.dumpOutput2( comments , "#     Got normal field:" );
 			messageWriter( "Point weight / Estimated Area: %g / %g\n" , pointWeightSum , pointCount*pointWeightSum );
 		}
@@ -674,7 +693,7 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 		if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Tree.value );
 		FEMTree< Dim , Real >::WriteParameter( fp );
 		DenseNodeData< Real , Sigs >::WriteSignatures( fp );
-		tree.write( fp );
+		tree.write( fp , xForm );
 		solution.write( fp );
 		fclose( fp );
 	}
diff --git a/Src/SparseMatrix.h b/Src/SparseMatrix.h
index 44a9c2c..ee2487c 100644
--- a/Src/SparseMatrix.h
+++ b/Src/SparseMatrix.h
@@ -25,8 +25,8 @@ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 DAMAGE.
 */
-#ifndef __SPARSEMATRIX_HPP
-#define __SPARSEMATRIX_HPP
+#ifndef SPARSE_MATRIX_INCLUDED
+#define SPARSE_MATRIX_INCLUDED
 
 #include "SparseMatrixInterface.h"
 #include "Array.h"
@@ -151,4 +151,4 @@ public:
 };
 
 #include "SparseMatrix.inl"
-#endif /* __SPARSEMATRIX_HPP */
+#endif // SPARSE_MATRIX_INCLUDED