From be731665f0f923af118e4e6a9661998c46216517 Mon Sep 17 00:00:00 2001
From: HiroIshida <h-ishida@jsk.imi.i.u-tokyo.ac.jp>
Date: Mon, 20 Jan 2025 02:20:37 +0900
Subject: [PATCH] chore: add SDFGen (MIT-licensed) source code under third

---
 third/SDFGen/.gitignore        |   2 +
 third/SDFGen/CMakeLists.txt    |  13 +
 third/SDFGen/README.md         |  25 ++
 third/SDFGen/array1.h          | 794 +++++++++++++++++++++++++++++++++
 third/SDFGen/array3.h          | 272 +++++++++++
 third/SDFGen/makelevelset3.cpp | 184 ++++++++
 third/SDFGen/makelevelset3.h   |  16 +
 third/SDFGen/util.h            | 470 +++++++++++++++++++
 third/SDFGen/vec.h             | 478 ++++++++++++++++++++
 9 files changed, 2254 insertions(+)
 create mode 100644 third/SDFGen/.gitignore
 create mode 100644 third/SDFGen/CMakeLists.txt
 create mode 100644 third/SDFGen/README.md
 create mode 100644 third/SDFGen/array1.h
 create mode 100644 third/SDFGen/array3.h
 create mode 100644 third/SDFGen/makelevelset3.cpp
 create mode 100644 third/SDFGen/makelevelset3.h
 create mode 100644 third/SDFGen/util.h
 create mode 100644 third/SDFGen/vec.h

diff --git a/third/SDFGen/.gitignore b/third/SDFGen/.gitignore
new file mode 100644
index 0000000..6bace84
--- /dev/null
+++ b/third/SDFGen/.gitignore
@@ -0,0 +1,2 @@
+action.sh
+/build
diff --git a/third/SDFGen/CMakeLists.txt b/third/SDFGen/CMakeLists.txt
new file mode 100644
index 0000000..9e3c50c
--- /dev/null
+++ b/third/SDFGen/CMakeLists.txt
@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 2.6)
+Project("SDFGen")
+
+SET(CMAKE_BUILD_TYPE Release)
+
+set(EXECUTABLE_OUTPUT_PATH ${CMAKE_CURRENT_BINARY_DIR}/bin)
+
+set(CMAKE_CXX_FLAGS "-O3")
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native")
+# static
+add_library(sdfgen STATIC makelevelset3.cpp)
+# add_executable(${PROJECT_NAME} main.cpp  makelevelset3.cpp)
diff --git a/third/SDFGen/README.md b/third/SDFGen/README.md
new file mode 100644
index 0000000..ed85c98
--- /dev/null
+++ b/third/SDFGen/README.md
@@ -0,0 +1,25 @@
+# SDFGen
+A simple commandline utility to generate grid-based signed distance field (level set) generator from triangle meshes, using code from Robert Bridson's website.
+
+
+The MIT License (MIT)
+
+Copyright (c) 2015, Christopher Batty
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/third/SDFGen/array1.h b/third/SDFGen/array1.h
new file mode 100644
index 0000000..9c67b49
--- /dev/null
+++ b/third/SDFGen/array1.h
@@ -0,0 +1,794 @@
+#ifndef ARRAY1_H
+#define ARRAY1_H
+
+#include <algorithm>
+#include <cstring>
+#include <cassert>
+#include <climits>
+#include <cstdlib>
+#include <iostream>
+#include <stdexcept>
+#include <vector>
+
+// In this file:
+//   Array1<T>: a dynamic 1D array for plain-old-data (not objects)
+//   WrapArray1<T>: a 1D array wrapper around an existing array (perhaps objects, perhaps data)
+// For the most part std::vector operations are supported, though for the Wrap version
+// note that memory is never allocated/deleted and constructor/destructors are never called
+// from within the class, thus only shallow copies can be made and some operations such as
+// resize() and push_back() are limited.
+// Note: for the most part assertions are done with assert(), not exceptions...
+
+// gross template hacking to determine if a type is integral or not
+struct Array1True {};
+struct Array1False {};
+template<typename T> struct Array1IsIntegral{ typedef Array1False type; }; // default: no (specializations to yes follow)
+template<> struct Array1IsIntegral<bool>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<char>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<signed char>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned char>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<short>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned short>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<int>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned int>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<long>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned long>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<long long>{ typedef Array1True type; };
+template<> struct Array1IsIntegral<unsigned long long>{ typedef Array1True type; };
+
+//============================================================================
+template<typename T>
+struct Array1
+{
+   // STL-friendly typedefs
+
+   typedef T* iterator;
+   typedef const T* const_iterator;
+   typedef unsigned long size_type;
+   typedef long difference_type;
+   typedef T& reference;
+   typedef const T& const_reference;
+   typedef T value_type;
+   typedef T* pointer;
+   typedef const T* const_pointer;
+   typedef std::reverse_iterator<iterator> reverse_iterator;
+   typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
+
+   // the actual representation
+
+   unsigned long n;
+   unsigned long max_n;
+   T* data;
+
+   // STL vector's interface, with additions, but only valid when used with plain-old-data
+
+   Array1(void)
+      : n(0), max_n(0), data(0)
+   {}
+
+   // note: default initial values are zero
+   Array1(unsigned long n_)
+      : n(0), max_n(0), data(0)
+   {
+      if(n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+      data=(T*)std::calloc(n_, sizeof(T));
+      if(!data) throw std::bad_alloc();
+      n=n_;
+      max_n=n_;
+   }
+
+   Array1(unsigned long n_, const T& value)
+      : n(0), max_n(0), data(0)
+   {
+      if(n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+      data=(T*)std::calloc(n_, sizeof(T));
+      if(!data) throw std::bad_alloc();
+      n=n_;
+      max_n=n_;
+      for(unsigned long i=0; i<n; ++i) data[i]=value;
+   }
+
+   Array1(unsigned long n_, const T& value, unsigned long max_n_)
+      : n(0), max_n(0), data(0)
+   {
+      assert(n_<=max_n_);
+      if(max_n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+      data=(T*)std::calloc(max_n_, sizeof(T));
+      if(!data) throw std::bad_alloc();
+      n=n_;
+      max_n=max_n_;
+      for(unsigned long i=0; i<n; ++i) data[i]=value;
+   }
+
+   Array1(unsigned long n_, const T* data_)
+      : n(0), max_n(0), data(0)
+   {
+      if(n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+      data=(T*)std::calloc(n_, sizeof(T));
+      if(!data) throw std::bad_alloc();
+      n=n_;
+      max_n=n_;
+      assert(data_);
+      std::memcpy(data, data_, n*sizeof(T));
+   }
+
+   Array1(unsigned long n_, const T* data_, unsigned long max_n_)
+      : n(0), max_n(0), data(0)
+   {
+      assert(n_<=max_n_);
+      if(max_n_>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+      data=(T*)std::calloc(max_n_, sizeof(T));
+      if(!data) throw std::bad_alloc();
+      max_n=max_n_;
+      n=n_;
+      assert(data_);
+      std::memcpy(data, data_, n*sizeof(T));
+   }
+
+   Array1(const Array1<T> &x)
+      : n(0), max_n(0), data(0)
+   {
+      data=(T*)std::malloc(x.n*sizeof(T));
+      if(!data) throw std::bad_alloc();
+      n=x.n;
+      max_n=x.n;
+      std::memcpy(data, x.data, n*sizeof(T));
+   }
+
+   ~Array1(void)
+   {
+      std::free(data);
+#ifndef NDEBUG
+      data=0;
+      n=max_n=0;
+#endif
+   }
+
+   const T& operator[](unsigned long i) const
+   { return data[i]; }
+
+   T& operator[](unsigned long i)
+   { return data[i]; }
+
+   // these are range-checked (in debug mode) versions of operator[], like at()
+   const T& operator()(unsigned long i) const
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   T& operator()(unsigned long i)
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   Array1<T>& operator=(const Array1<T>& x)
+   {
+      if(max_n<x.n){
+         T* new_data=(T*)std::malloc(x.n*sizeof(T));
+         if(!new_data) throw std::bad_alloc();
+         std::free(data);
+         data=new_data;
+         max_n=x.n;
+      }
+      n=x.n;
+      std::memcpy(data, x.data, n*sizeof(T));
+      return *this;
+   }
+
+   bool operator==(const Array1<T>& x) const
+   {
+      if(n!=x.n) return false;
+      for(unsigned long i=0; i<n; ++i) if(!(data[i]==x.data[i])) return false;
+      return true;
+   }
+ 
+   bool operator!=(const Array1<T>& x) const
+   {
+      if(n!=x.n) return true;
+      for(unsigned long i=0; i<n; ++i) if(data[i]!=x.data[i]) return true;
+      return false;
+   }
+
+   bool operator<(const Array1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]<x[i]) return true;
+         else if(x[i]<data[i]) return false;
+      }
+      return n<x.n;
+   }
+
+   bool operator>(const Array1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]>x[i]) return true;
+         else if(x[i]>data[i]) return false;
+      }
+      return n>x.n;
+   }
+
+   bool operator<=(const Array1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]<x[i]) return true;
+         else if(x[i]<data[i]) return false;
+      }
+      return n<=x.n;
+   }
+
+   bool operator>=(const Array1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]>x[i]) return true;
+         else if(x[i]>data[i]) return false;
+      }
+      return n>=x.n;
+   }
+
+   void add_unique(const T& value)
+   {
+      for(unsigned long i=0; i<n; ++i) if(data[i]==value) return;
+      if(n==max_n) grow();
+      data[n++]=value;
+   }
+
+   void assign(const T& value)
+   { for(unsigned long i=0; i<n; ++i) data[i]=value; }
+
+   void assign(unsigned long num, const T& value)
+   { fill(num, value); } 
+
+   // note: copydata may not alias this array's data, and this should not be
+   // used when T is a full object (which defines its own copying operation)
+   void assign(unsigned long num, const T* copydata)
+   {
+      assert(num==0 || copydata);
+      if(num>max_n){
+         if(num>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+         std::free(data);
+         data=(T*)std::malloc(num*sizeof(T));
+         if(!data) throw std::bad_alloc();
+         max_n=num;
+      }
+      n=num;
+      std::memcpy(data, copydata, n*sizeof(T));
+   }
+
+   template<typename InputIterator>
+   void assign(InputIterator first, InputIterator last)
+   { assign_(first, last, typename Array1IsIntegral<InputIterator>::type()); }
+
+   template<typename InputIterator>
+   void assign_(InputIterator first, InputIterator last, Array1True check)
+   { fill(first, last); }
+
+   template<typename InputIterator>
+   void assign_(InputIterator first, InputIterator last, Array1False check)
+   {
+      unsigned long i=0;
+      InputIterator p=first;
+      for(; p!=last; ++p, ++i){
+         if(i==max_n) grow();
+         data[i]=*p;
+      }
+      n=i;
+   }
+
+   const T& at(unsigned long i) const
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   T& at(unsigned long i)
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   const T& back(void) const
+   { 
+      assert(data && n>0);
+      return data[n-1];
+   }
+
+   T& back(void)
+   {
+      assert(data && n>0);
+      return data[n-1];
+   }
+
+   const T* begin(void) const
+   { return data; }
+
+   T* begin(void)
+   { return data; }
+
+   unsigned long capacity(void) const
+   { return max_n; }
+
+   void clear(void)
+   {
+      std::free(data);
+      data=0;
+      max_n=0;
+      n=0;
+   }
+
+   bool empty(void) const
+   { return n==0; }
+
+   const T* end(void) const
+   { return data+n; }
+
+   T* end(void)
+   { return data+n; }
+
+   void erase(unsigned long index)
+   {
+      assert(index<n);
+      for(unsigned long i=index; i<n-1; ++i)
+         data[i]=data[i-1];
+      pop_back();
+   }
+
+   void fill(unsigned long num, const T& value)
+   {
+      if(num>max_n){
+         if(num>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+         std::free(data);
+         data=(T*)std::malloc(num*sizeof(T));
+         if(!data) throw std::bad_alloc();
+         max_n=num;
+      }
+      n=num;
+      for(unsigned long i=0; i<n; ++i) data[i]=value;
+   }
+
+   const T& front(void) const
+   {
+      assert(n>0);
+      return *data;
+   }
+
+   T& front(void)
+   {
+      assert(n>0);
+      return *data;
+   }
+
+   void grow(void)
+   {
+      unsigned long new_size=(max_n*sizeof(T)<ULONG_MAX/2 ? 2*max_n+1 : ULONG_MAX/sizeof(T));
+      T *new_data=(T*)std::realloc(data, new_size*sizeof(T));
+      if(!new_data) throw std::bad_alloc();
+      data=new_data;
+      max_n=new_size;
+   }
+
+   void insert(unsigned long index, const T& entry)
+   {
+      assert(index<=n);
+      push_back(back());
+      for(unsigned long i=n-1; i>index; --i)
+         data[i]=data[i-1];
+      data[index]=entry;
+   }
+
+   unsigned long max_size(void) const
+   { return ULONG_MAX/sizeof(T); }
+
+   void pop_back(void)
+   {
+      assert(n>0);
+      --n;
+   }
+
+   void push_back(const T& value)
+   {
+      if(n==max_n) grow();
+      data[n++]=value;
+   }
+
+   reverse_iterator rbegin(void)
+   { return reverse_iterator(end()); }
+
+   const_reverse_iterator rbegin(void) const
+   { return const_reverse_iterator(end()); }
+
+   reverse_iterator rend(void)
+   { return reverse_iterator(begin()); }
+
+   const_reverse_iterator rend(void) const
+   { return const_reverse_iterator(begin()); }
+
+   void reserve(unsigned long r)
+   {
+      if(r>ULONG_MAX/sizeof(T)) throw std::bad_alloc();
+      T *new_data=(T*)std::realloc(data, r*sizeof(T));
+      if(!new_data) throw std::bad_alloc();
+      data=new_data;
+      max_n=r;
+   }
+
+   void resize(unsigned long n_)
+   {
+      if(n_>max_n) reserve(n_);
+      n=n_;
+   }
+
+   void resize(unsigned long n_, const T& value)
+   {
+      if(n_>max_n) reserve(n_);
+      if(n<n_) for(unsigned long i=n; i<n_; ++i) data[i]=value;
+      n=n_;
+   }
+
+   void set_zero(void)
+   { std::memset(data, 0, n*sizeof(T)); }
+
+   unsigned long size(void) const
+   { return n; }
+
+   void swap(Array1<T>& x)
+   {
+      std::swap(n, x.n);
+      std::swap(max_n, x.max_n);
+      std::swap(data, x.data);
+   }
+
+   // resize the array to avoid wasted space, without changing contents
+   // (Note: realloc, at least on some platforms, will not do the trick)
+   void trim(void)
+   {
+      if(n==max_n) return;
+      T *new_data=(T*)std::malloc(n*sizeof(T));
+      if(!new_data) return;
+      std::memcpy(new_data, data, n*sizeof(T));
+      std::free(data);
+      data=new_data;
+      max_n=n;
+   }
+};
+
+// some common arrays
+
+typedef Array1<double>             Array1d;
+typedef Array1<float>              Array1f;
+typedef Array1<long long>          Array1ll;
+typedef Array1<unsigned long long> Array1ull;
+typedef Array1<int>                Array1i;
+typedef Array1<unsigned int>       Array1ui;
+typedef Array1<short>              Array1s;
+typedef Array1<unsigned short>     Array1us;
+typedef Array1<char>               Array1c;
+typedef Array1<unsigned char>      Array1uc;
+
+//============================================================================
+template<typename T>
+struct WrapArray1
+{
+   // STL-friendly typedefs
+
+   typedef T* iterator;
+   typedef const T* const_iterator;
+   typedef unsigned long size_type;
+   typedef long difference_type;
+   typedef T& reference;
+   typedef const T& const_reference;
+   typedef T value_type;
+   typedef T* pointer;
+   typedef const T* const_pointer;
+   typedef std::reverse_iterator<iterator> reverse_iterator;
+   typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
+
+   // the actual representation
+
+   unsigned long n;
+   unsigned long max_n;
+   T* data;
+
+   // most of STL vector's interface, with a few changes
+
+   WrapArray1(void)
+      : n(0), max_n(0), data(0)
+   {}
+
+   WrapArray1(unsigned long n_, T* data_)
+      : n(n_), max_n(n_), data(data_)
+   { assert(data || max_n==0); }
+
+   WrapArray1(unsigned long n_, T* data_, unsigned long max_n_)
+      : n(n_), max_n(max_n_), data(data_)
+   {
+      assert(n<=max_n);
+      assert(data || max_n==0);
+   }
+
+   // Allow for simple shallow copies of existing arrays
+   // Note that if the underlying arrays change where their data is, the WrapArray may be screwed up
+
+   WrapArray1(Array1<T>& a)
+      : n(a.n), max_n(a.max_n), data(a.data)
+   {}
+
+   WrapArray1(std::vector<T>& a)
+      : n(a.size()), max_n(a.capacity()), data(&a[0])
+   {}
+
+   void init(unsigned long n_, T* data_, unsigned long max_n_)
+   {
+      assert(n_<=max_n_);
+      assert(data_ || max_n_==0);
+      n=n_;
+      max_n=max_n_;
+      data=data_;
+   }
+
+   const T& operator[](unsigned long i) const
+   { return data[i]; }
+
+   T& operator[](unsigned long i)
+   { return data[i]; }
+
+   // these are range-checked (in debug mode) versions of operator[], like at()
+   const T& operator()(unsigned long i) const
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   T& operator()(unsigned long i)
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   bool operator==(const WrapArray1<T>& x) const
+   {
+      if(n!=x.n) return false;
+      for(unsigned long i=0; i<n; ++i) if(!(data[i]==x.data[i])) return false;
+      return true;
+   }
+ 
+   bool operator!=(const WrapArray1<T>& x) const
+   {
+      if(n!=x.n) return true;
+      for(unsigned long i=0; i<n; ++i) if(data[i]!=x.data[i]) return true;
+      return false;
+   }
+
+   bool operator<(const WrapArray1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]<x[i]) return true;
+         else if(x[i]<data[i]) return false;
+      }
+      return n<x.n;
+   }
+
+   bool operator>(const WrapArray1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]>x[i]) return true;
+         else if(x[i]>data[i]) return false;
+      }
+      return n>x.n;
+   }
+
+   bool operator<=(const WrapArray1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]<x[i]) return true;
+         else if(x[i]<data[i]) return false;
+      }
+      return n<=x.n;
+   }
+
+   bool operator>=(const WrapArray1<T>& x) const
+   {
+      for(unsigned long i=0; i<n && i<x.n; ++i){
+         if(data[i]>x[i]) return true;
+         else if(x[i]>data[i]) return false;
+      }
+      return n>=x.n;
+   }
+
+   void add_unique(const T& value)
+   {
+      for(unsigned long i=0; i<n; ++i) if(data[i]==value) return;
+      assert(n<max_n);
+      data[n++]=value;
+   }
+
+   void assign(const T& value)
+   { for(unsigned long i=0; i<n; ++i) data[i]=value; }
+
+   void assign(unsigned long num, const T& value)
+   { fill(num, value); } 
+
+   // note: copydata may not alias this array's data, and this should not be
+   // used when T is a full object (which defines its own copying operation)
+   void assign(unsigned long num, const T* copydata)
+   {
+      assert(num==0 || copydata);
+      assert(num<=max_n);
+      n=num;
+      std::memcpy(data, copydata, n*sizeof(T));
+   }
+
+   template<typename InputIterator>
+   void assign(InputIterator first, InputIterator last)
+   { assign_(first, last, typename Array1IsIntegral<InputIterator>::type()); }
+
+   template<typename InputIterator>
+   void assign_(InputIterator first, InputIterator last, Array1True check)
+   { fill(first, last); }
+
+   template<typename InputIterator>
+   void assign_(InputIterator first, InputIterator last, Array1False check)
+   {
+      unsigned long i=0;
+      InputIterator p=first;
+      for(; p!=last; ++p, ++i){
+         assert(i<max_n);
+         data[i]=*p;
+      }
+      n=i;
+   }
+
+   const T& at(unsigned long i) const
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   T& at(unsigned long i)
+   {
+      assert(i<n);
+      return data[i];
+   }
+
+   const T& back(void) const
+   { 
+      assert(data && n>0);
+      return data[n-1];
+   }
+
+   T& back(void)
+   {
+      assert(data && n>0);
+      return data[n-1];
+   }
+
+   const T* begin(void) const
+   { return data; }
+
+   T* begin(void)
+   { return data; }
+
+   unsigned long capacity(void) const
+   { return max_n; }
+
+   void clear(void)
+   { n=0; }
+
+   bool empty(void) const
+   { return n==0; }
+
+   const T* end(void) const
+   { return data+n; }
+
+   T* end(void)
+   { return data+n; }
+
+   void erase(unsigned long index)
+   {
+      assert(index<n);
+      for(unsigned long i=index; i<n-1; ++i)
+         data[i]=data[i-1];
+      pop_back();
+   }
+
+   void fill(unsigned long num, const T& value)
+   {
+      assert(num<=max_n);
+      n=num;
+      for(unsigned long i=0; i<n; ++i) data[i]=value;
+   }
+
+   const T& front(void) const
+   {
+      assert(n>0);
+      return *data;
+   }
+
+   T& front(void)
+   {
+      assert(n>0);
+      return *data;
+   }
+
+   void insert(unsigned long index, const T& entry)
+   {
+      assert(index<=n);
+      push_back(back());
+      for(unsigned long i=n-1; i>index; --i)
+         data[i]=data[i-1];
+      data[index]=entry;
+   }
+
+   unsigned long max_size(void) const
+   { return max_n; }
+
+   void pop_back(void)
+   {
+      assert(n>0);
+      --n;
+   }
+
+   void push_back(const T& value)
+   {
+      assert(n<max_n);
+      data[n++]=value;
+   }
+
+   reverse_iterator rbegin(void)
+   { return reverse_iterator(end()); }
+
+   const_reverse_iterator rbegin(void) const
+   { return const_reverse_iterator(end()); }
+
+   reverse_iterator rend(void)
+   { return reverse_iterator(begin()); }
+
+   const_reverse_iterator rend(void) const
+   { return const_reverse_iterator(begin()); }
+
+   void reserve(unsigned long r)
+   { assert(r<=max_n); }
+
+   void resize(unsigned long n_)
+   {
+      assert(n_<=max_n);
+      n=n_;
+   }
+
+   void resize(unsigned long n_, const T& value)
+   {
+      assert(n_<=max_n);
+      if(n<n_) for(unsigned long i=n; i<n_; ++i) data[i]=value;
+      n=n_;
+   }
+
+   // note: shouldn't be used when T is a full object (setting to zero may not make sense)
+   void set_zero(void)
+   { std::memset(data, 0, n*sizeof(T)); }
+
+   unsigned long size(void) const
+   { return n; }
+
+   void swap(WrapArray1<T>& x)
+   {
+      std::swap(n, x.n);
+      std::swap(max_n, x.max_n);
+      std::swap(data, x.data);
+   }
+};
+
+// some common arrays
+
+typedef WrapArray1<double>             WrapArray1d;
+typedef WrapArray1<float>              WrapArray1f;
+typedef WrapArray1<long long>          WrapArray1ll;
+typedef WrapArray1<unsigned long long> WrapArray1ull;
+typedef WrapArray1<int>                WrapArray1i;
+typedef WrapArray1<unsigned int>       WrapArray1ui;
+typedef WrapArray1<short>              WrapArray1s;
+typedef WrapArray1<unsigned short>     WrapArray1us;
+typedef WrapArray1<char>               WrapArray1c;
+typedef WrapArray1<unsigned char>      WrapArray1uc;
+
+#endif
diff --git a/third/SDFGen/array3.h b/third/SDFGen/array3.h
new file mode 100644
index 0000000..92028ae
--- /dev/null
+++ b/third/SDFGen/array3.h
@@ -0,0 +1,272 @@
+#ifndef ARRAY3_H
+#define ARRAY3_H
+
+#include "array1.h"
+#include <algorithm>
+#include <cassert>
+#include <vector>
+
+template<class T, class ArrayT=std::vector<T> >
+struct Array3
+{
+   // STL-friendly typedefs
+
+   typedef typename ArrayT::iterator iterator;
+   typedef typename ArrayT::const_iterator const_iterator;
+   typedef typename ArrayT::size_type size_type;
+   typedef long difference_type;
+   typedef T& reference;
+   typedef const T& const_reference;
+   typedef T value_type;
+   typedef T* pointer;
+   typedef const T* const_pointer;
+   typedef typename ArrayT::reverse_iterator reverse_iterator;
+   typedef typename ArrayT::const_reverse_iterator const_reverse_iterator;
+
+   // the actual representation
+
+   int ni, nj, nk;
+   ArrayT a;
+
+   // the interface
+
+   Array3(void)
+      : ni(0), nj(0), nk(0)
+   {}
+
+   Array3(int ni_, int nj_, int nk_)
+      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_)
+   { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+   Array3(int ni_, int nj_, int nk_, ArrayT& a_)
+      : ni(ni_), nj(nj_), nk(nk_), a(a_)
+   { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+   Array3(int ni_, int nj_, int nk_, const T& value)
+      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value)
+   { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+   Array3(int ni_, int nj_, int nk_, const T& value, size_type max_n_)
+      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, value, max_n_)
+   { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+   Array3(int ni_, int nj_, int nk_, T* data_)
+      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_)
+   { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+   Array3(int ni_, int nj_, int nk_, T* data_, size_type max_n_)
+      : ni(ni_), nj(nj_), nk(nk_), a(ni_*nj_*nk_, data_, max_n_)
+   { assert(ni_>=0 && nj_>=0 && nk_>=0); }
+
+   ~Array3(void)
+   {
+#ifndef NDEBUG
+      ni=nj=0;
+#endif
+   }
+
+   const T& operator()(int i, int j, int k) const
+   {
+      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+      return a[i+ni*(j+nj*k)];
+   }
+
+   T& operator()(int i, int j, int k)
+   {
+      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+      return a[i+ni*(j+nj*k)];
+   }
+
+   bool operator==(const Array3<T>& x) const
+   { return ni==x.ni && nj==x.nj && nk==x.nk && a==x.a; }
+
+   bool operator!=(const Array3<T>& x) const
+   { return ni!=x.ni || nj!=x.nj || nk!=x.nk || a!=x.a; }
+
+   bool operator<(const Array3<T>& x) const
+   {
+      if(ni<x.ni) return true; else if(ni>x.ni) return false;
+      if(nj<x.nj) return true; else if(nj>x.nj) return false;
+      if(nk<x.nk) return true; else if(nk>x.nk) return false;
+      return a<x.a;
+   }
+
+   bool operator>(const Array3<T>& x) const
+   {
+      if(ni>x.ni) return true; else if(ni<x.ni) return false;
+      if(nj>x.nj) return true; else if(nj<x.nj) return false;
+      if(nk>x.nk) return true; else if(nk<x.nk) return false;
+      return a>x.a;
+   }
+
+   bool operator<=(const Array3<T>& x) const
+   {
+      if(ni<x.ni) return true; else if(ni>x.ni) return false;
+      if(nj<x.nj) return true; else if(nj>x.nj) return false;
+      if(nk<x.nk) return true; else if(nk>x.nk) return false;
+      return a<=x.a;
+   }
+
+   bool operator>=(const Array3<T>& x) const
+   {
+      if(ni>x.ni) return true; else if(ni<x.ni) return false;
+      if(nj>x.nj) return true; else if(nj<x.nj) return false;
+      if(nk>x.nk) return true; else if(nk<x.nk) return false;
+      return a>=x.a;
+   }
+
+   void assign(const T& value)
+   { a.assign(value); }
+
+   void assign(int ni_, int nj_, int nk_, const T& value)
+   {
+      a.assign(ni_*nj_*nk_, value);
+      ni=ni_;
+      nj=nj_;
+      nk=nk_;
+   }
+    
+   void assign(int ni_, int nj_, int nk_, const T* copydata)
+   {
+      a.assign(ni_*nj_*nk_, copydata);
+      ni=ni_;
+      nj=nj_;
+      nk=nk_;
+   }
+    
+   const T& at(int i, int j, int k) const
+   {
+      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+      return a[i+ni*(j+nj*k)];
+   }
+
+   T& at(int i, int j, int k)
+   {
+      assert(i>=0 && i<ni && j>=0 && j<nj && k>=0 && k<nk);
+      return a[i+ni*(j+nj*k)];
+   }
+
+   const T& back(void) const
+   { 
+      assert(a.size());
+      return a.back();
+   }
+
+   T& back(void)
+   {
+      assert(a.size());
+      return a.back();
+   }
+
+   const_iterator begin(void) const
+   { return a.begin(); }
+
+   iterator begin(void)
+   { return a.begin(); }
+
+   size_type capacity(void) const
+   { return a.capacity(); }
+
+   void clear(void)
+   {
+      a.clear();
+      ni=nj=nk=0;
+   }
+
+   bool empty(void) const
+   { return a.empty(); }
+
+   const_iterator end(void) const
+   { return a.end(); }
+
+   iterator end(void)
+   { return a.end(); }
+
+   void fill(int ni_, int nj_, int nk_, const T& value)
+   {
+      a.fill(ni_*nj_*nk_, value);
+      ni=ni_;
+      nj=nj_;
+      nk=nk_;
+   }
+    
+   const T& front(void) const
+   {
+      assert(a.size());
+      return a.front();
+   }
+
+   T& front(void)
+   {
+      assert(a.size());
+      return a.front();
+   }
+
+   size_type max_size(void) const
+   { return a.max_size(); }
+
+   reverse_iterator rbegin(void)
+   { return reverse_iterator(end()); }
+
+   const_reverse_iterator rbegin(void) const
+   { return const_reverse_iterator(end()); }
+
+   reverse_iterator rend(void)
+   { return reverse_iterator(begin()); }
+
+   const_reverse_iterator rend(void) const
+   { return const_reverse_iterator(begin()); }
+
+   void reserve(int reserve_ni, int reserve_nj, int reserve_nk)
+   { a.reserve(reserve_ni*reserve_nj*reserve_nk); }
+
+   void resize(int ni_, int nj_, int nk_)
+   {
+      assert(ni_>=0 && nj_>=0 && nk_>=0);
+      a.resize(ni_*nj_*nk_);
+      ni=ni_;
+      nj=nj_;
+      nk=nk_;
+   }
+
+   void resize(int ni_, int nj_, int nk_, const T& value)
+   {
+      assert(ni_>=0 && nj_>=0 && nk_>=0);
+      a.resize(ni_*nj_*nk_, value);
+      ni=ni_;
+      nj=nj_;
+      nk=nk_;
+   }
+
+   void set_zero(void)
+   { a.set_zero(); }
+
+   size_type size(void) const
+   { return a.size(); }
+
+   void swap(Array3<T>& x)
+   {
+      std::swap(ni, x.ni);
+      std::swap(nj, x.nj);
+      std::swap(nk, x.nk);
+      a.swap(x.a);
+   }
+
+   void trim(void)
+   { a.trim(); }
+};
+
+// some common arrays
+
+typedef Array3<double, Array1<double> > Array3d;
+typedef Array3<float, Array1<float> > Array3f;
+typedef Array3<long long, Array1<long long> > Array3ll;
+typedef Array3<unsigned long long, Array1<unsigned long long> > Array3ull;
+typedef Array3<int, Array1<int> > Array3i;
+typedef Array3<unsigned int, Array1<unsigned int> > Array3ui;
+typedef Array3<short, Array1<short> > Array3s;
+typedef Array3<unsigned short, Array1<unsigned short> > Array3us;
+typedef Array3<char, Array1<char> > Array3c;
+typedef Array3<unsigned char, Array1<unsigned char> > Array3uc;
+
+#endif
diff --git a/third/SDFGen/makelevelset3.cpp b/third/SDFGen/makelevelset3.cpp
new file mode 100644
index 0000000..92f2139
--- /dev/null
+++ b/third/SDFGen/makelevelset3.cpp
@@ -0,0 +1,184 @@
+#include "makelevelset3.h"
+
+// find distance x0 is from segment x1-x2
+static float point_segment_distance(const Vec3f &x0, const Vec3f &x1, const Vec3f &x2)
+{
+   Vec3f dx(x2-x1);
+   double m2=mag2(dx);
+   // find parameter value of closest point on segment
+   float s12=(float)(dot(x2-x0, dx)/m2);
+   if(s12<0){
+      s12=0;
+   }else if(s12>1){
+      s12=1;
+   }
+   // and find the distance
+   return dist(x0, s12*x1+(1-s12)*x2);
+}
+
+// find distance x0 is from triangle x1-x2-x3
+static float point_triangle_distance(const Vec3f &x0, const Vec3f &x1, const Vec3f &x2, const Vec3f &x3)
+{
+   // first find barycentric coordinates of closest point on infinite plane
+   Vec3f x13(x1-x3), x23(x2-x3), x03(x0-x3);
+   float m13=mag2(x13), m23=mag2(x23), d=dot(x13,x23);
+   float invdet=1.f/max(m13*m23-d*d,1e-30f);
+   float a=dot(x13,x03), b=dot(x23,x03);
+   // the barycentric coordinates themselves
+   float w23=invdet*(m23*a-d*b);
+   float w31=invdet*(m13*b-d*a);
+   float w12=1-w23-w31;
+   if(w23>=0 && w31>=0 && w12>=0){ // if we're inside the triangle
+      return dist(x0, w23*x1+w31*x2+w12*x3); 
+   }else{ // we have to clamp to one of the edges
+      if(w23>0) // this rules out edge 2-3 for us
+         return min(point_segment_distance(x0,x1,x2), point_segment_distance(x0,x1,x3));
+      else if(w31>0) // this rules out edge 1-3
+         return min(point_segment_distance(x0,x1,x2), point_segment_distance(x0,x2,x3));
+      else // w12 must be >0, ruling out edge 1-2
+         return min(point_segment_distance(x0,x1,x3), point_segment_distance(x0,x2,x3));
+   }
+}
+
+static void check_neighbour(const std::vector<Vec3ui> &tri, const std::vector<Vec3f> &x,
+                            Array3f &phi, Array3i &closest_tri,
+                            const Vec3f &gx, int i0, int j0, int k0, int i1, int j1, int k1)
+{
+   if(closest_tri(i1,j1,k1)>=0){
+      unsigned int p, q, r; assign(tri[closest_tri(i1,j1,k1)], p, q, r);
+      float d=point_triangle_distance(gx, x[p], x[q], x[r]);
+      if(d<phi(i0,j0,k0)){
+         phi(i0,j0,k0)=d;
+         closest_tri(i0,j0,k0)=closest_tri(i1,j1,k1);
+      }
+   }
+}
+
+static void sweep(const std::vector<Vec3ui> &tri, const std::vector<Vec3f> &x,
+                  Array3f &phi, Array3i &closest_tri, const Vec3f &origin, float dx,
+                  int di, int dj, int dk)
+{
+   int i0, i1;
+   if(di>0){ i0=1; i1=phi.ni; }
+   else{ i0=phi.ni-2; i1=-1; }
+   int j0, j1;
+   if(dj>0){ j0=1; j1=phi.nj; }
+   else{ j0=phi.nj-2; j1=-1; }
+   int k0, k1;
+   if(dk>0){ k0=1; k1=phi.nk; }
+   else{ k0=phi.nk-2; k1=-1; }
+   for(int k=k0; k!=k1; k+=dk) for(int j=j0; j!=j1; j+=dj) for(int i=i0; i!=i1; i+=di){
+      Vec3f gx(i*dx+origin[0], j*dx+origin[1], k*dx+origin[2]);
+      check_neighbour(tri, x, phi, closest_tri, gx, i, j, k, i-di, j,    k);
+      check_neighbour(tri, x, phi, closest_tri, gx, i, j, k, i,    j-dj, k);
+      check_neighbour(tri, x, phi, closest_tri, gx, i, j, k, i-di, j-dj, k);
+      check_neighbour(tri, x, phi, closest_tri, gx, i, j, k, i,    j,    k-dk);
+      check_neighbour(tri, x, phi, closest_tri, gx, i, j, k, i-di, j,    k-dk);
+      check_neighbour(tri, x, phi, closest_tri, gx, i, j, k, i,    j-dj, k-dk);
+      check_neighbour(tri, x, phi, closest_tri, gx, i, j, k, i-di, j-dj, k-dk);
+   }
+}
+
+// calculate twice signed area of triangle (0,0)-(x1,y1)-(x2,y2)
+// return an SOS-determined sign (-1, +1, or 0 only if it's a truly degenerate triangle)
+static int orientation(double x1, double y1, double x2, double y2, double &twice_signed_area)
+{
+   twice_signed_area=y1*x2-x1*y2;
+   if(twice_signed_area>0) return 1;
+   else if(twice_signed_area<0) return -1;
+   else if(y2>y1) return 1;
+   else if(y2<y1) return -1;
+   else if(x1>x2) return 1;
+   else if(x1<x2) return -1;
+   else return 0; // only true when x1==x2 and y1==y2
+}
+
+// robust test of (x0,y0) in the triangle (x1,y1)-(x2,y2)-(x3,y3)
+// if true is returned, the barycentric coordinates are set in a,b,c.
+static bool point_in_triangle_2d(double x0, double y0, 
+                                 double x1, double y1, double x2, double y2, double x3, double y3,
+                                 double& a, double& b, double& c)
+{
+   x1-=x0; x2-=x0; x3-=x0;
+   y1-=y0; y2-=y0; y3-=y0;
+   int signa=orientation(x2, y2, x3, y3, a);
+   if(signa==0) return false;
+   int signb=orientation(x3, y3, x1, y1, b);
+   if(signb!=signa) return false;
+   int signc=orientation(x1, y1, x2, y2, c);
+   if(signc!=signa) return false;
+   double sum=a+b+c;
+   assert(sum!=0); // if the SOS signs match and are nonkero, there's no way all of a, b, and c are zero.
+   a/=sum;
+   b/=sum;
+   c/=sum;
+   return true;
+}
+
+void make_level_set3(const std::vector<Vec3ui> &tri, const std::vector<Vec3f> &x,
+                     const Vec3f &origin, float dx, int ni, int nj, int nk,
+                     Array3f &phi, const int exact_band)
+{
+   phi.resize(ni, nj, nk);
+   phi.assign((ni+nj+nk)*dx); // upper bound on distance
+   Array3i closest_tri(ni, nj, nk, -1);
+   Array3i intersection_count(ni, nj, nk, 0); // intersection_count(i,j,k) is # of tri intersections in (i-1,i]x{j}x{k}
+   // we begin by initializing distances near the mesh, and figuring out intersection counts
+   Vec3f ijkmin, ijkmax;
+   for(unsigned int t=0; t<tri.size(); ++t){
+     unsigned int p, q, r; assign(tri[t], p, q, r);
+     // coordinates in grid to high precision
+      double fip=((double)x[p][0]-origin[0])/dx, fjp=((double)x[p][1]-origin[1])/dx, fkp=((double)x[p][2]-origin[2])/dx;
+      double fiq=((double)x[q][0]-origin[0])/dx, fjq=((double)x[q][1]-origin[1])/dx, fkq=((double)x[q][2]-origin[2])/dx;
+      double fir=((double)x[r][0]-origin[0])/dx, fjr=((double)x[r][1]-origin[1])/dx, fkr=((double)x[r][2]-origin[2])/dx;
+      // do distances nearby
+      int i0=clamp(int(min(fip,fiq,fir))-exact_band, 0, ni-1), i1=clamp(int(max(fip,fiq,fir))+exact_band+1, 0, ni-1);
+      int j0=clamp(int(min(fjp,fjq,fjr))-exact_band, 0, nj-1), j1=clamp(int(max(fjp,fjq,fjr))+exact_band+1, 0, nj-1);
+      int k0=clamp(int(min(fkp,fkq,fkr))-exact_band, 0, nk-1), k1=clamp(int(max(fkp,fkq,fkr))+exact_band+1, 0, nk-1);
+      for(int k=k0; k<=k1; ++k) for(int j=j0; j<=j1; ++j) for(int i=i0; i<=i1; ++i){
+         Vec3f gx(i*dx+origin[0], j*dx+origin[1], k*dx+origin[2]);
+         float d=point_triangle_distance(gx, x[p], x[q], x[r]);
+         if(d<phi(i,j,k)){
+            phi(i,j,k)=d;
+            closest_tri(i,j,k)=t;
+         }
+      }
+      // and do intersection counts
+      j0=clamp((int)std::ceil(min(fjp,fjq,fjr)), 0, nj-1);
+      j1=clamp((int)std::floor(max(fjp,fjq,fjr)), 0, nj-1);
+      k0=clamp((int)std::ceil(min(fkp,fkq,fkr)), 0, nk-1);
+      k1=clamp((int)std::floor(max(fkp,fkq,fkr)), 0, nk-1);
+      for(int k=k0; k<=k1; ++k) for(int j=j0; j<=j1; ++j){
+         double a, b, c;
+         if(point_in_triangle_2d(j, k, fjp, fkp, fjq, fkq, fjr, fkr, a, b, c)){
+            double fi=a*fip+b*fiq+c*fir; // intersection i coordinate
+            int i_interval=int(std::ceil(fi)); // intersection is in (i_interval-1,i_interval]
+            if(i_interval<0) ++intersection_count(0, j, k); // we enlarge the first interval to include everything to the -x direction
+            else if(i_interval<ni) ++intersection_count(i_interval,j,k);
+            // we ignore intersections that are beyond the +x side of the grid
+         }
+      }
+   }
+   // and now we fill in the rest of the distances with fast sweeping
+   for(unsigned int pass=0; pass<2; ++pass){
+      sweep(tri, x, phi, closest_tri, origin, dx, +1, +1, +1);
+      sweep(tri, x, phi, closest_tri, origin, dx, -1, -1, -1);
+      sweep(tri, x, phi, closest_tri, origin, dx, +1, +1, -1);
+      sweep(tri, x, phi, closest_tri, origin, dx, -1, -1, +1);
+      sweep(tri, x, phi, closest_tri, origin, dx, +1, -1, +1);
+      sweep(tri, x, phi, closest_tri, origin, dx, -1, +1, -1);
+      sweep(tri, x, phi, closest_tri, origin, dx, +1, -1, -1);
+      sweep(tri, x, phi, closest_tri, origin, dx, -1, +1, +1);
+   }
+   // then figure out signs (inside/outside) from intersection counts
+   for(int k=0; k<nk; ++k) for(int j=0; j<nj; ++j){
+      int total_count=0;
+      for(int i=0; i<ni; ++i){
+         total_count+=intersection_count(i,j,k);
+         if(total_count%2==1){ // if parity of intersections so far is odd,
+            phi(i,j,k)=-phi(i,j,k); // we are inside the mesh
+         }
+      }
+   }
+}
+
diff --git a/third/SDFGen/makelevelset3.h b/third/SDFGen/makelevelset3.h
new file mode 100644
index 0000000..b1d7d1d
--- /dev/null
+++ b/third/SDFGen/makelevelset3.h
@@ -0,0 +1,16 @@
+#ifndef MAKELEVELSET3_H
+#define MAKELEVELSET3_H
+
+#include "array3.h"
+#include "vec.h"
+
+// tri is a list of triangles in the mesh, and x is the positions of the vertices
+// absolute distances will be nearly correct for triangle soup, but a closed mesh is
+// needed for accurate signs. Distances for all grid cells within exact_band cells of
+// a triangle should be exact; further away a distance is calculated but it might not
+// be to the closest triangle - just one nearby.
+void make_level_set3(const std::vector<Vec3ui> &tri, const std::vector<Vec3f> &x,
+                     const Vec3f &origin, float dx, int nx, int ny, int nz,
+                     Array3f &phi, const int exact_band=1);
+
+#endif
diff --git a/third/SDFGen/util.h b/third/SDFGen/util.h
new file mode 100644
index 0000000..c5758a1
--- /dev/null
+++ b/third/SDFGen/util.h
@@ -0,0 +1,470 @@
+#ifndef UTIL_H
+#define UTIL_H
+
+#include <algorithm>
+#include <vector>
+#include <cmath>
+#include <iostream>
+
+#ifndef M_PI
+const double M_PI = 3.1415926535897932384626433832795;
+#endif
+
+#ifdef WIN32
+#undef min
+#undef max
+#endif
+
+using std::min;
+using std::max;
+using std::swap;
+
+template<class T>
+inline T sqr(const T& x)
+{ return x*x; }
+
+template<class T>
+inline T cube(const T& x)
+{ return x*x*x; }
+
+template<class T>
+inline T min(T a1, T a2, T a3)
+{ return min(a1, min(a2, a3)); }
+
+template<class T>
+inline T min(T a1, T a2, T a3, T a4)
+{ return min(min(a1, a2), min(a3, a4)); }
+
+template<class T>
+inline T min(T a1, T a2, T a3, T a4, T a5)
+{ return min(min(a1, a2), min(a3, a4), a5); }
+
+template<class T>
+inline T min(T a1, T a2, T a3, T a4, T a5, T a6)
+{ return min(min(a1, a2), min(a3, a4), min(a5, a6)); }
+
+template<class T>
+inline T max(T a1, T a2, T a3)
+{ return max(a1, max(a2, a3)); }
+
+template<class T>
+inline T max(T a1, T a2, T a3, T a4)
+{ return max(max(a1, a2), max(a3, a4)); }
+
+template<class T>
+inline T max(T a1, T a2, T a3, T a4, T a5)
+{ return max(max(a1, a2), max(a3, a4),  a5); }
+
+template<class T>
+inline T max(T a1, T a2, T a3, T a4, T a5, T a6)
+{ return max(max(a1, a2), max(a3, a4),  max(a5, a6)); }
+
+template<class T>
+inline void minmax(T a1, T a2, T& amin, T& amax)
+{
+   if(a1<a2){
+      amin=a1;
+      amax=a2;
+   }else{
+      amin=a2;
+      amax=a1;
+   }
+}
+
+template<class T>
+inline void minmax(T a1, T a2, T a3, T& amin, T& amax)
+{
+   if(a1<a2){
+      if(a1<a3){
+         amin=a1;
+         if(a2<a3) amax=a3;
+         else amax=a2;
+      }else{
+         amin=a3;
+         if(a1<a2) amax=a2;
+         else amax=a1;
+      }
+   }else{
+      if(a2<a3){
+         amin=a2;
+         if(a1<a3) amax=a3;
+         else amax=a1;
+      }else{
+         amin=a3;
+         amax=a1;
+      }
+   }
+}
+
+template<class T>
+inline void minmax(T a1, T a2, T a3, T a4, T& amin, T& amax)
+{
+   if(a1<a2){
+      if(a3<a4){
+         amin=min(a1,a3);
+         amax=max(a2,a4);
+      }else{
+         amin=min(a1,a4);
+         amax=max(a2,a3);
+      }
+   }else{
+      if(a3<a4){
+         amin=min(a2,a3);
+         amax=max(a1,a4);
+      }else{
+         amin=min(a2,a4);
+         amax=max(a1,a3);
+      }
+   }
+}
+
+template<class T>
+inline void minmax(T a1, T a2, T a3, T a4, T a5, T& amin, T& amax)
+{
+   //@@@ the logic could be shortcircuited a lot!
+   amin=min(a1,a2,a3,a4,a5);
+   amax=max(a1,a2,a3,a4,a5);
+}
+
+template<class T>
+inline void minmax(T a1, T a2, T a3, T a4, T a5, T a6, T& amin, T& amax)
+{
+   //@@@ the logic could be shortcircuited a lot!
+   amin=min(a1,a2,a3,a4,a5,a6);
+   amax=max(a1,a2,a3,a4,a5,a6);
+}
+
+template<class T>
+inline void update_minmax(T a1, T& amin, T& amax)
+{
+   if(a1<amin) amin=a1;
+   else if(a1>amax) amax=a1;
+}
+
+template<class T>
+inline void sort(T &a, T &b, T &c)
+{
+   T temp;
+   if(a<b){
+      if(a<c){
+    if(c<b){ // a<c<b
+            temp=c;c=b;b=temp;
+    } // else: a<b<c
+      }else{ // c<a<b
+    temp=c;c=b;b=a;a=temp;
+      }
+   }else{
+      if(b<c){
+    if(a<c){ //b<a<c
+       temp=b;b=a;a=temp;
+    }else{ // b<c<a
+       temp=b;b=c;c=a;a=temp;
+    }
+      }else{ // c<b<a
+    temp=c;c=a;a=temp;
+      }
+   }
+}
+
+template<class T>
+inline T clamp(T a, T lower, T upper)
+{
+   if(a<lower) return lower;
+   else if(a>upper) return upper;
+   else return a;
+}
+
+// only makes sense with T=float or double
+template<class T>
+inline T smooth_step(T r)
+{
+   if(r<0) return 0;
+   else if(r>1) return 1;
+   return r*r*r*(10+r*(-15+r*6));
+}
+
+// only makes sense with T=float or double
+template<class T>
+inline T smooth_step(T r, T r_lower, T r_upper, T value_lower, T value_upper)
+{ return value_lower + smooth_step((r-r_lower)/(r_upper-r_lower)) * (value_upper-value_lower); }
+
+// only makes sense with T=float or double
+template<class T>
+inline T ramp(T r)
+{ return smooth_step((r+1)/2)*2-1; }
+
+#ifdef WIN32
+inline int lround(double x)
+{
+   if(x>0)
+      return (x-floor(x)<0.5) ? (int)floor(x) : (int)ceil(x);
+   else
+      return (x-floor(x)<=0.5) ? (int)floor(x) : (int)ceil(x);
+}
+
+inline double remainder(double x, double y)
+{
+   return x-std::floor(x/y+0.5)*y;
+}
+#endif
+
+inline unsigned int round_up_to_power_of_two(unsigned int n)
+{
+   int exponent=0;
+   --n;
+   while(n){
+      ++exponent;
+      n>>=1;
+   }
+   return 1<<exponent;
+}
+
+inline unsigned int round_down_to_power_of_two(unsigned int n)
+{
+   int exponent=0;
+   while(n>1){
+      ++exponent;
+      n>>=1;
+   }
+   return 1<<exponent;
+}
+
+// Transforms even the sequence 0,1,2,3,... into reasonably good random numbers 
+// Challenge: improve on this in speed and "randomness"!
+// This seems to pass several statistical tests, and is a bijective map (of 32-bit unsigned ints)
+inline unsigned int randhash(unsigned int seed)
+{
+   unsigned int i=(seed^0xA3C59AC3u)*2654435769u;
+   i^=(i>>16);
+   i*=2654435769u;
+   i^=(i>>16);
+   i*=2654435769u;
+   return i;
+}
+
+// the inverse of randhash
+inline unsigned int unhash(unsigned int h)
+{
+   h*=340573321u;
+   h^=(h>>16);
+   h*=340573321u;
+   h^=(h>>16);
+   h*=340573321u;
+   h^=0xA3C59AC3u;
+   return h;
+}
+
+// returns repeatable stateless pseudo-random number in [0,1]
+inline double randhashd(unsigned int seed)
+{ return randhash(seed)/(double)UINT_MAX; }
+inline float randhashf(unsigned int seed)
+{ return randhash(seed)/(float)UINT_MAX; }
+
+// returns repeatable stateless pseudo-random number in [a,b]
+inline double randhashd(unsigned int seed, double a, double b)
+{ return (b-a)*randhash(seed)/(double)UINT_MAX + a; }
+inline float randhashf(unsigned int seed, float a, float b)
+{ return ( (b-a)*randhash(seed)/(float)UINT_MAX + a); }
+
+inline int intlog2(int x)
+{
+   int exp=-1;
+   while(x){
+      x>>=1;
+      ++exp;
+   }
+   return exp;
+}
+
+template<class T>
+inline void get_barycentric(T x, int& i, T& f, int i_low, int i_high)
+{
+   T s=std::floor(x);
+   i=(int)s;
+   if(i<i_low){
+      i=i_low;
+      f=0;
+   }else if(i>i_high-2){
+      i=i_high-2;
+      f=1;
+   }else
+      f=(T)(x-s);
+}
+
+template<class S, class T>
+inline S lerp(const S& value0, const S& value1, T f)
+{ return (1-f)*value0 + f*value1; }
+
+template<class S, class T>
+inline S bilerp(const S& v00, const S& v10, 
+                const S& v01, const S& v11, 
+                T fx, T fy)
+{ 
+   return lerp(lerp(v00, v10, fx),
+               lerp(v01, v11, fx), 
+               fy);
+}
+
+template<class S, class T>
+inline S trilerp(const S& v000, const S& v100,
+                 const S& v010, const S& v110,
+                 const S& v001, const S& v101,  
+                 const S& v011, const S& v111,
+                 T fx, T fy, T fz) 
+{
+   return lerp(bilerp(v000, v100, v010, v110, fx, fy),
+               bilerp(v001, v101, v011, v111, fx, fy),
+               fz);
+}
+
+template<class S, class T>
+inline S quadlerp(const S& v0000, const S& v1000,
+                  const S& v0100, const S& v1100,
+                  const S& v0010, const S& v1010,  
+                  const S& v0110, const S& v1110,
+                  const S& v0001, const S& v1001,
+                  const S& v0101, const S& v1101,
+                  const S& v0011, const S& v1011,  
+                  const S& v0111, const S& v1111,
+                  T fx, T fy, T fz, T ft) 
+{
+   return lerp(trilerp(v0000, v1000, v0100, v1100, v0010, v1010, v0110, v1110, fx, fy, fz),
+               trilerp(v0001, v1001, v0101, v1101, v0011, v1011, v0111, v1111, fx, fy, fz),
+               ft);
+}
+
+// f should be between 0 and 1, with f=0.5 corresponding to balanced weighting between w0 and w2
+template<class T>
+inline void quadratic_bspline_weights(T f, T& w0, T& w1, T& w2)
+{
+   w0=T(0.5)*sqr(f-1);
+   w1=T(0.75)-sqr(f-T(0.5));;
+   w2=T(0.5)*sqr(f);
+}
+
+// f should be between 0 and 1
+template<class T>
+inline void cubic_interp_weights(T f, T& wneg1, T& w0, T& w1, T& w2)
+{
+   T f2(f*f), f3(f2*f);
+   wneg1=-T(1./3)*f+T(1./2)*f2-T(1./6)*f3;
+   w0=1-f2+T(1./2)*(f3-f);
+   w1=f+T(1./2)*(f2-f3);
+   w2=T(1./6)*(f3-f);
+}
+
+template<class S, class T>
+inline S cubic_interp(const S& value_neg1, const S& value0, const S& value1, const S& value2, T f)
+{
+   T wneg1, w0, w1, w2;
+   cubic_interp_weights(f, wneg1, w0, w1, w2);
+   return wneg1*value_neg1 + w0*value0 + w1*value1 + w2*value2;
+}
+
+template<class T>
+void zero(std::vector<T>& v)
+{ for(int i=(int)v.size()-1; i>=0; --i) v[i]=0; }
+
+template<class T>
+T abs_max(const std::vector<T>& v)
+{
+   T m=0;
+   for(int i=(int)v.size()-1; i>=0; --i){
+      if(std::fabs(v[i])>m)
+         m=std::fabs(v[i]);
+   }
+   return m;
+}
+
+template<class T>
+bool contains(const std::vector<T>& a, T e)
+{
+   for(unsigned int i=0; i<a.size(); ++i)
+      if(a[i]==e) return true;
+   return false;
+}
+
+template<class T>
+void add_unique(std::vector<T>& a, T e)
+{
+   for(unsigned int i=0; i<a.size(); ++i)
+      if(a[i]==e) return;
+   a.push_back(e);
+}
+
+template<class T>
+void insert(std::vector<T>& a, unsigned int index, T e)
+{
+   a.push_back(a.back());
+   for(unsigned int i=(unsigned int)a.size()-1; i>index; --i)
+      a[i]=a[i-1];
+   a[index]=e;
+}
+
+template<class T>
+void erase(std::vector<T>& a, unsigned int index)
+{
+   for(unsigned int i=index; i<a.size()-1; ++i)
+      a[i]=a[i+1];
+   a.pop_back();
+}
+
+template<class T>
+void erase_swap(std::vector<T>& a, unsigned int index)
+{
+   for(unsigned int i=index; i<a.size()-1; ++i)
+      swap(a[i], a[i+1]);
+   a.pop_back();
+}
+
+template<class T>
+void erase_unordered(std::vector<T>& a, unsigned int index)
+{
+   a[index]=a.back();
+   a.pop_back();
+}
+
+template<class T>
+void erase_unordered_swap(std::vector<T>& a, unsigned int index)
+{
+   swap(a[index], a.back());
+   a.pop_back();
+}
+
+template<class T>
+void find_and_erase_unordered(std::vector<T>& a, const T& doomed_element)
+{
+   for(unsigned int i=0; i<a.size(); ++i)
+      if(a[i]==doomed_element){
+         erase_unordered(a, i);
+         return;
+      }
+}
+
+template<class T>
+void replace_once(std::vector<T>& a, const T& old_element, const T& new_element)
+{
+   for(unsigned int i=0; i<a.size(); ++i)
+      if(a[i]==old_element){
+         a[i]=new_element;
+         return;
+      }
+}
+
+template<class T>
+void write_matlab(std::ostream& output, const std::vector<T>& a, const char *variable_name, bool column_vector=true, int significant_digits=18)
+{
+   output<<variable_name<<"=[";
+   std::streamsize old_precision=output.precision();
+   output.precision(significant_digits);
+   for(unsigned int i=0; i<a.size(); ++i){
+      output<<a[i]<<" ";
+   }
+   output<<"]";
+   if(column_vector)
+      output<<"'";
+   output<<";"<<std::endl;
+   output.precision(old_precision);
+}
+
+#endif
diff --git a/third/SDFGen/vec.h b/third/SDFGen/vec.h
new file mode 100644
index 0000000..d846af6
--- /dev/null
+++ b/third/SDFGen/vec.h
@@ -0,0 +1,478 @@
+#ifndef VEC_H
+#define VEC_H
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include "util.h"
+
+// Defines a thin wrapper around fixed size C-style arrays, using template parameters,
+// which is useful for dealing with vectors of different dimensions.
+// For example, float[3] is equivalent to Vec<3,float>.
+// Entries in the vector are accessed with the overloaded [] operator, so
+// for example if x is a Vec<3,float>, then the middle entry is x[1].
+// For convenience, there are a number of typedefs for abbreviation:
+//   Vec<3,float> -> Vec3f
+//   Vec<2,int>   -> Vec2i
+// and so on.
+// Arithmetic operators are appropriately overloaded, and functions are defined
+// for additional operations (such as dot-products, norms, cross-products, etc.)
+
+template<unsigned int N, class T>
+struct Vec
+{
+   T v[N];
+
+   Vec<N,T>(void)
+   {}
+
+   explicit Vec<N,T>(T value_for_all)
+   { for(unsigned int i=0; i<N; ++i) v[i]=value_for_all; }
+
+   template<class S>
+   explicit Vec<N,T>(const S *source)
+   { for(unsigned int i=0; i<N; ++i) v[i]=(T)source[i]; }
+
+   template <class S>
+   explicit Vec<N,T>(const Vec<N,S>& source)
+   { for(unsigned int i=0; i<N; ++i) v[i]=(T)source[i]; }
+
+   Vec<N,T>(T v0, T v1)
+   {
+      assert(N==2);
+      v[0]=v0; v[1]=v1;
+   }
+
+   Vec<N,T>(T v0, T v1, T v2)
+   {
+      assert(N==3);
+      v[0]=v0; v[1]=v1; v[2]=v2;
+   }
+
+   Vec<N,T>(T v0, T v1, T v2, T v3)
+   {
+      assert(N==4);
+      v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3;
+   }
+
+   Vec<N,T>(T v0, T v1, T v2, T v3, T v4)
+   {
+      assert(N==5);
+      v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4;
+   }
+
+     Vec<N,T>(T v0, T v1, T v2, T v3, T v4, T v5)
+   {
+      assert(N==6);
+      v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; v[4]=v4; v[5]=v5;
+   }
+
+   T &operator[](int index)
+   {
+      assert(0<=index && (unsigned int)index<N);
+      return v[index];
+   }
+
+   const T &operator[](int index) const
+   {
+      assert(0<=index && (unsigned int)index<N);
+      return v[index];
+   }
+
+   bool nonzero(void) const
+   {
+      for(unsigned int i=0; i<N; ++i) if(v[i]) return true;
+      return false;
+   }
+
+   Vec<N,T> operator+=(const Vec<N,T> &w)
+   {
+      for(unsigned int i=0; i<N; ++i) v[i]+=w[i];
+      return *this;
+   }
+
+   Vec<N,T> operator+(const Vec<N,T> &w) const
+   {
+      Vec<N,T> sum(*this);
+      sum+=w;
+      return sum;
+   }
+
+   Vec<N,T> operator-=(const Vec<N,T> &w)
+   {
+      for(unsigned int i=0; i<N; ++i) v[i]-=w[i];
+      return *this;
+   }
+
+   Vec<N,T> operator-(void) const // unary minus
+   {
+      Vec<N,T> negative;
+      for(unsigned int i=0; i<N; ++i) negative.v[i]=-v[i];
+      return negative;
+   }
+
+   Vec<N,T> operator-(const Vec<N,T> &w) const // (binary) subtraction
+   {
+      Vec<N,T> diff(*this);
+      diff-=w;
+      return diff;
+   }
+
+   Vec<N,T> operator*=(T a)
+   {
+      for(unsigned int i=0; i<N; ++i) v[i]*=a;
+      return *this;
+   }
+
+   Vec<N,T> operator*(T a) const
+   {
+      Vec<N,T> w(*this);
+      w*=a;
+      return w;
+   }
+
+   Vec<N,T> operator*=(const Vec<N,T> &w)
+   {
+      for(unsigned int i=0; i<N; ++i) v[i]*=w.v[i];
+      return *this;
+   }
+
+   Vec<N,T> operator*(const Vec<N,T> &w) const
+   {
+      Vec<N,T> componentwise_product;
+      for(unsigned int i=0; i<N; ++i) componentwise_product[i]=v[i]*w.v[i];
+      return componentwise_product;
+   }
+
+   Vec<N,T> operator/=(T a)
+   {
+      for(unsigned int i=0; i<N; ++i) v[i]/=a;
+      return *this;
+   }
+
+   Vec<N,T> operator/(T a) const
+   {
+      Vec<N,T> w(*this);
+      w/=a;
+      return w;
+   }
+};
+
+typedef Vec<2,double>         Vec2d;
+typedef Vec<2,float>          Vec2f;
+typedef Vec<2,int>            Vec2i;
+typedef Vec<2,unsigned int>   Vec2ui;
+typedef Vec<2,short>          Vec2s;
+typedef Vec<2,unsigned short> Vec2us;
+typedef Vec<2,char>           Vec2c;
+typedef Vec<2,unsigned char>  Vec2uc;
+
+typedef Vec<3,double>         Vec3d;
+typedef Vec<3,float>          Vec3f;
+typedef Vec<3,int>            Vec3i;
+typedef Vec<3,unsigned int>   Vec3ui;
+typedef Vec<3,short>          Vec3s;
+typedef Vec<3,unsigned short> Vec3us;
+typedef Vec<3,char>           Vec3c;
+typedef Vec<3,unsigned char>  Vec3uc;
+
+typedef Vec<4,double>         Vec4d;
+typedef Vec<4,float>          Vec4f;
+typedef Vec<4,int>            Vec4i;
+typedef Vec<4,unsigned int>   Vec4ui;
+typedef Vec<4,short>          Vec4s;
+typedef Vec<4,unsigned short> Vec4us;
+typedef Vec<4,char>           Vec4c;
+typedef Vec<4,unsigned char>  Vec4uc;
+
+typedef Vec<6,double>         Vec6d;
+typedef Vec<6,float>          Vec6f;
+typedef Vec<6,unsigned int>   Vec6ui;
+typedef Vec<6,int>            Vec6i;
+typedef Vec<6,short>          Vec6s;
+typedef Vec<6,unsigned short> Vec6us;
+typedef Vec<6,char>           Vec6c;
+typedef Vec<6,unsigned char>  Vec6uc;
+
+
+template<unsigned int N, class T>
+T mag2(const Vec<N,T> &a)
+{
+   T l=sqr(a.v[0]);
+   for(unsigned int i=1; i<N; ++i) l+=sqr(a.v[i]);
+   return l;
+}
+
+template<unsigned int N, class T>
+T mag(const Vec<N,T> &a)
+{ return sqrt(mag2(a)); }
+
+template<unsigned int N, class T> 
+inline T dist2(const Vec<N,T> &a, const Vec<N,T> &b)
+{ 
+   T d=sqr(a.v[0]-b.v[0]);
+   for(unsigned int i=1; i<N; ++i) d+=sqr(a.v[i]-b.v[i]);
+   return d;
+}
+
+template<unsigned int N, class T> 
+inline T dist(const Vec<N,T> &a, const Vec<N,T> &b)
+{ return std::sqrt(dist2(a,b)); }
+
+template<unsigned int N, class T> 
+inline void normalize(Vec<N,T> &a)
+{ a/=mag(a); }
+
+template<unsigned int N, class T> 
+inline Vec<N,T> normalized(const Vec<N,T> &a)
+{ return a/mag(a); }
+
+template<unsigned int N, class T> 
+inline T infnorm(const Vec<N,T> &a)
+{
+   T d=std::fabs(a.v[0]);
+   for(unsigned int i=1; i<N; ++i) d=max(std::fabs(a.v[i]),d);
+   return d;
+}
+
+template<unsigned int N, class T>
+void zero(Vec<N,T> &a)
+{ 
+   for(unsigned int i=0; i<N; ++i)
+      a.v[i] = 0;
+}
+
+template<unsigned int N, class T>
+std::ostream &operator<<(std::ostream &out, const Vec<N,T> &v)
+{
+   out<<v.v[0];
+   for(unsigned int i=1; i<N; ++i)
+      out<<' '<<v.v[i];
+   return out;
+}
+
+template<unsigned int N, class T>
+std::istream &operator>>(std::istream &in, Vec<N,T> &v)
+{
+   in>>v.v[0];
+   for(unsigned int i=1; i<N; ++i)
+      in>>v.v[i];
+   return in;
+}
+
+template<unsigned int N, class T> 
+inline bool operator==(const Vec<N,T> &a, const Vec<N,T> &b)
+{ 
+   bool t = (a.v[0] == b.v[0]);
+   unsigned int i=1;
+   while(i<N && t) {
+      t = t && (a.v[i]==b.v[i]); 
+      ++i;
+   }
+   return t;
+}
+
+template<unsigned int N, class T> 
+inline bool operator!=(const Vec<N,T> &a, const Vec<N,T> &b)
+{ 
+   bool t = (a.v[0] != b.v[0]);
+   unsigned int i=1;
+   while(i<N && !t) {
+      t = t || (a.v[i]!=b.v[i]); 
+      ++i;
+   }
+   return t;
+}
+
+template<unsigned int N, class T>
+inline Vec<N,T> operator*(T a, const Vec<N,T> &v)
+{
+   Vec<N,T> w(v);
+   w*=a;
+   return w;
+}
+
+template<unsigned int N, class T>
+inline T min(const Vec<N,T> &a)
+{
+   T m=a.v[0];
+   for(unsigned int i=1; i<N; ++i) if(a.v[i]<m) m=a.v[i];
+   return m;
+}
+
+template<unsigned int N, class T>
+inline Vec<N,T> min_union(const Vec<N,T> &a, const Vec<N,T> &b)
+{
+   Vec<N,T> m;
+   for(unsigned int i=0; i<N; ++i) (a.v[i] < b.v[i]) ? m.v[i]=a.v[i] : m.v[i]=b.v[i];
+   return m;
+}
+
+template<unsigned int N, class T>
+inline Vec<N,T> max_union(const Vec<N,T> &a, const Vec<N,T> &b)
+{
+   Vec<N,T> m;
+   for(unsigned int i=0; i<N; ++i) (a.v[i] > b.v[i]) ? m.v[i]=a.v[i] : m.v[i]=b.v[i];
+   return m;
+}
+
+template<unsigned int N, class T>
+inline T max(const Vec<N,T> &a)
+{
+   T m=a.v[0];
+   for(unsigned int i=1; i<N; ++i) if(a.v[i]>m) m=a.v[i];
+   return m;
+}
+
+template<unsigned int N, class T>
+inline T dot(const Vec<N,T> &a, const Vec<N,T> &b)
+{
+   T d=a.v[0]*b.v[0];
+   for(unsigned int i=1; i<N; ++i) d+=a.v[i]*b.v[i];
+   return d;
+}
+
+template<class T> 
+inline Vec<2,T> rotate(const Vec<2,T>& a, float angle) 
+{
+   T c = cos(angle);
+   T s = sin(angle);
+   return Vec<2,T>(c*a[0] - s*a[1],s*a[0] + c*a[1]); // counter-clockwise rotation
+}
+
+template<class T>
+inline Vec<2,T> perp(const Vec<2,T> &a)
+{ return Vec<2,T>(-a.v[1], a.v[0]); } // counter-clockwise rotation by 90 degrees
+
+template<class T>
+inline T cross(const Vec<2,T> &a, const Vec<2,T> &b)
+{ return a.v[0]*b.v[1]-a.v[1]*b.v[0]; }
+
+template<class T>
+inline Vec<3,T> cross(const Vec<3,T> &a, const Vec<3,T> &b)
+{ return Vec<3,T>(a.v[1]*b.v[2]-a.v[2]*b.v[1], a.v[2]*b.v[0]-a.v[0]*b.v[2], a.v[0]*b.v[1]-a.v[1]*b.v[0]); }
+
+template<class T>
+inline T triple(const Vec<3,T> &a, const Vec<3,T> &b, const Vec<3,T> &c)
+{ return a.v[0]*(b.v[1]*c.v[2]-b.v[2]*c.v[1])
+        +a.v[1]*(b.v[2]*c.v[0]-b.v[0]*c.v[2])
+        +a.v[2]*(b.v[0]*c.v[1]-b.v[1]*c.v[0]); }
+
+template<unsigned int N, class T>
+inline unsigned int hash(const Vec<N,T> &a)
+{
+   unsigned int h=a.v[0];
+   for(unsigned int i=1; i<N; ++i)
+      h=hash(h ^ a.v[i]);
+   return h;
+}
+
+template<unsigned int N, class T>
+inline void assign(const Vec<N,T> &a, T &a0, T &a1)
+{ 
+   assert(N==2);
+   a0=a.v[0]; a1=a.v[1];
+}
+
+template<unsigned int N, class T>
+inline void assign(const Vec<N,T> &a, T &a0, T &a1, T &a2)
+{ 
+   assert(N==3);
+   a0=a.v[0]; a1=a.v[1]; a2=a.v[2];
+}
+
+template<unsigned int N, class T>
+inline void assign(const Vec<N,T> &a, T &a0, T &a1, T &a2, T &a3)
+{ 
+   assert(N==4);
+   a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; a3=a.v[3];
+}
+
+template<unsigned int N, class T>
+inline void assign(const Vec<N,T> &a, T &a0, T &a1, T &a2, T &a3, T &a4, T &a5)
+{ 
+   assert(N==6);
+   a0=a.v[0]; a1=a.v[1]; a2=a.v[2]; a3=a.v[3]; a4=a.v[4]; a5=a.v[5];
+}
+
+template<unsigned int N, class T>
+inline Vec<N,int> round(const Vec<N,T> &a)
+{ 
+   Vec<N,int> rounded;
+   for(unsigned int i=0; i<N; ++i)
+      rounded.v[i]=lround(a.v[i]);
+   return rounded; 
+}
+
+template<unsigned int N, class T>
+inline Vec<N,int> floor(const Vec<N,T> &a)
+{ 
+   Vec<N,int> rounded;
+   for(unsigned int i=0; i<N; ++i)
+      rounded.v[i]=(int)floor(a.v[i]);
+   return rounded; 
+}
+
+template<unsigned int N, class T>
+inline Vec<N,int> ceil(const Vec<N,T> &a)
+{ 
+   Vec<N,int> rounded;
+   for(unsigned int i=0; i<N; ++i)
+      rounded.v[i]=(int)ceil(a.v[i]);
+   return rounded; 
+}
+
+template<unsigned int N, class T>
+inline Vec<N,T> fabs(const Vec<N,T> &a)
+{ 
+   Vec<N,T> result;
+   for(unsigned int i=0; i<N; ++i)
+      result.v[i]=fabs(a.v[i]);
+   return result; 
+}
+
+template<unsigned int N, class T>
+inline void minmax(const Vec<N,T> &x0, const Vec<N,T> &x1, Vec<N,T> &xmin, Vec<N,T> &xmax)
+{
+   for(unsigned int i=0; i<N; ++i)
+      minmax(x0.v[i], x1.v[i], xmin.v[i], xmax.v[i]);
+}
+
+template<unsigned int N, class T>
+inline void minmax(const Vec<N,T> &x0, const Vec<N,T> &x1, const Vec<N,T> &x2, Vec<N,T> &xmin, Vec<N,T> &xmax)
+{
+   for(unsigned int i=0; i<N; ++i)
+      minmax(x0.v[i], x1.v[i], x2.v[i], xmin.v[i], xmax.v[i]);
+}
+
+template<unsigned int N, class T>
+inline void minmax(const Vec<N,T> &x0, const Vec<N,T> &x1, const Vec<N,T> &x2, const Vec<N,T> &x3,
+                   Vec<N,T> &xmin, Vec<N,T> &xmax)
+{
+   for(unsigned int i=0; i<N; ++i)
+      minmax(x0.v[i], x1.v[i], x2.v[i], x3.v[i], xmin.v[i], xmax.v[i]);
+}
+
+template<unsigned int N, class T>
+inline void minmax(const Vec<N,T> &x0, const Vec<N,T> &x1, const Vec<N,T> &x2, const Vec<N,T> &x3, const Vec<N,T> &x4,
+                   Vec<N,T> &xmin, Vec<N,T> &xmax)
+{
+   for(unsigned int i=0; i<N; ++i)
+      minmax(x0.v[i], x1.v[i], x2.v[i], x3.v[i], x4.v[i], xmin.v[i], xmax.v[i]);
+}
+
+template<unsigned int N, class T>
+inline void minmax(const Vec<N,T> &x0, const Vec<N,T> &x1, const Vec<N,T> &x2, const Vec<N,T> &x3, const Vec<N,T> &x4,
+                   const Vec<N,T> &x5, Vec<N,T> &xmin, Vec<N,T> &xmax)
+{
+   for(unsigned int i=0; i<N; ++i)
+      minmax(x0.v[i], x1.v[i], x2.v[i], x3.v[i], x4.v[i], x5.v[i], xmin.v[i], xmax.v[i]);
+}
+
+template<unsigned int N, class T>
+inline void update_minmax(const Vec<N,T> &x, Vec<N,T> &xmin, Vec<N,T> &xmax)
+{
+   for(unsigned int i=0; i<N; ++i) update_minmax(x[i], xmin[i], xmax[i]);
+}
+
+#endif