rm sdf

17 files changed, 18 insertions, 3180 deletions

diff --git a/extern/softbody/CMakeLists.txt b/extern/softbody/CMakeLists.txt
index 84b2de7fea7..78e4525322b 100644
--- a/extern/softbody/CMakeLists.txt
+++ b/extern/softbody/CMakeLists.txt
@@ -29,15 +29,15 @@ set(INC_SYS
 )
 
 set(SDFGEN_SRC
-	src/sdfgen/array1.h
-	src/sdfgen/array2.h
-	src/sdfgen/array3.h
-	src/sdfgen/hashgrid.h
-	src/sdfgen/hashtable.h
-	src/sdfgen/makelevelset3.cpp
-	src/sdfgen/makelevelset3.h
-	src/sdfgen/util.h
-	src/sdfgen/vec.h
+#	src/sdfgen/array1.h
+#	src/sdfgen/array2.h
+#	src/sdfgen/array3.h
+#	src/sdfgen/hashgrid.h
+#	src/sdfgen/hashtable.h
+#	src/sdfgen/makelevelset3.cpp
+#	src/sdfgen/makelevelset3.h
+#	src/sdfgen/util.h
+#	src/sdfgen/vec.h
 )
 
 set(SRC
diff --git a/extern/softbody/src/admmpd_bvh_traverse.cpp b/extern/softbody/src/admmpd_bvh_traverse.cpp
index 8fd4b1118f6..9ec4c281368 100644
--- a/extern/softbody/src/admmpd_bvh_traverse.cpp
+++ b/extern/softbody/src/admmpd_bvh_traverse.cpp
@@ -140,7 +140,6 @@ PointInTriangleMeshTraverse<T>::PointInTriangleMeshTraverse(
 	prim_inds(prim_inds_)
 {
 	//dir = VecType::Random();
-
 	BLI_assert(prim_verts->rows()>=0);
 	BLI_assert(prim_inds->rows()>=0);
 	BLI_assert(prim_inds->cols()==3);
diff --git a/extern/softbody/src/admmpd_sdf.cpp b/extern/softbody/src/admmpd_sdf.cpp
index 16bcc545903..f40615cdab6 100644
--- a/extern/softbody/src/admmpd_sdf.cpp
+++ b/extern/softbody/src/admmpd_sdf.cpp
@@ -1,6 +1,8 @@
 // Copyright Matt Overby 2020.
 // Distributed under the MIT License.
 
+#if 0
+
 #include "admmpd_sdf.h"
 #include "admmpd_geom.h"
 #include <vector>
@@ -344,4 +346,6 @@ void SDF<T>::get_cells(
 template class admmpd::SDF<double>;
 template class admmpd::SDF<float>;
 
-} // namespace admmpd
\ No newline at end of file
+} // namespace admmpd
+
+#endif
\ No newline at end of file
diff --git a/extern/softbody/src/admmpd_sdf.h b/extern/softbody/src/admmpd_sdf.h
index fce7689da88..898dbe3218c 100644
--- a/extern/softbody/src/admmpd_sdf.h
+++ b/extern/softbody/src/admmpd_sdf.h
@@ -1,6 +1,8 @@
 // Copyright Matt Overby 2020.
 // Distributed under the MIT License.
 
+#if 0
+
 #ifndef ADMMPD_SDF_H_
 #define ADMMPD_SDF_H_
 
@@ -84,3 +86,5 @@ public:
 } // namespace admmpd
 
 #endif // ADMMPD_SDF_H_
+
+#endif
\ No newline at end of file
diff --git a/extern/softbody/src/sdfgen/CMakeLists.txt b/extern/softbody/src/sdfgen/CMakeLists.txt
deleted file mode 100644
index 400c92beac7..00000000000
--- a/extern/softbody/src/sdfgen/CMakeLists.txt
+++ /dev/null
@@ -1,39 +0,0 @@
-cmake_minimum_required(VERSION 2.6)
-Project("SDFGen")
-
-# Set the build type.  Options are:
-#  Coverage       : w/ debug symbols, w/o optimization, w/ code-coverage
-#  Debug          : w/ debug symbols, w/o optimization
-#  Release        : w/o debug symbols, w/ optimization
-#  RelWithDebInfo : w/ debug symbols, w/ optimization
-#  MinSizeRel     : w/o debug symbols, w/ optimization, stripped binaries
-
-# SET(CMAKE_BUILD_TYPE Coverage)
-SET(CMAKE_BUILD_TYPE Release)
-
-#set the default path for built executables to the "bin" directory
-set(EXECUTABLE_OUTPUT_PATH ${CMAKE_CURRENT_BINARY_DIR}/bin)
-
-#These flags might not work on every system, especially the release flags, comment out as needed
-set(CMAKE_CXX_FLAGS "-O3")
-set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native")
-
-#checks if VTK is available
-find_package(VTK QUIET)
-if(VTK_FOUND)
-	set(HAVE_VTK 1)
-endif()
-
-#Generates config.h replacing expressions in config.h.in with actual values
-configure_file(
-  "${CMAKE_CURRENT_SOURCE_DIR}/config.h.in"
-  "${CMAKE_CURRENT_SOURCE_DIR}/config.h"
-)
-
-add_executable(${PROJECT_NAME} main.cpp  makelevelset3.cpp)
-
-if(VTK_FOUND)
-	include_directories(${VTK_INCLUDE_DIRS})
-	target_link_libraries(${PROJECT_NAME} ${VTK_LIBRARIES})
-endif()
\ No newline at end of file
diff --git a/extern/softbody/src/sdfgen/README.md b/extern/softbody/src/sdfgen/README.md
deleted file mode 100644
index ed85c981911..00000000000
--- a/extern/softbody/src/sdfgen/README.md
+++ /dev/null
@@ -1,25 +0,0 @@
-# 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/extern/softbody/src/sdfgen/array1.h b/extern/softbody/src/sdfgen/array1.h
deleted file mode 100644
index 08d39f2c284..00000000000
--- a/extern/softbody/src/sdfgen/array1.h
+++ /dev/null
@@ -1,798 +0,0 @@
-#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...
-
-namespace sdfgen {
-
-// 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/extern/softbody/src/sdfgen/array2.h b/extern/softbody/src/sdfgen/array2.h
deleted file mode 100644
index d08edd7c768..00000000000
--- a/extern/softbody/src/sdfgen/array2.h
+++ /dev/null
@@ -1,284 +0,0 @@
-#ifndef ARRAY2_H
-#define ARRAY2_H
-
-#include "array1.h"
-#include <algorithm>
-#include <cassert>
-#include <vector>
-
-namespace sdfgen {
-
-template<class T, class ArrayT=std::vector<T> >
-struct Array2
-{
-   // 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;
-   ArrayT a;
-
-   // the interface
-
-   Array2(void)
-      : ni(0), nj(0)
-   {}
-
-   Array2(int ni_, int nj_)
-      : ni(ni_), nj(nj_), a(ni_*nj_)
-   { assert(ni_>=0 && nj>=0); }
-
-   Array2(int ni_, int nj_, ArrayT& a_)
-      : ni(ni_), nj(nj_), a(a_)
-   { assert(ni_>=0 && nj>=0); }
-
-   Array2(int ni_, int nj_, const T& value)
-      : ni(ni_), nj(nj_), a(ni_*nj_, value)
-   { assert(ni_>=0 && nj>=0); }
-
-   Array2(int ni_, int nj_, const T& value, size_type max_n_)
-      : ni(ni_), nj(nj_), a(ni_*nj_, value, max_n_)
-   { assert(ni_>=0 && nj>=0); }
-
-   Array2(int ni_, int nj_, T* data_)
-      : ni(ni_), nj(nj_), a(ni_*nj_, data_)
-   { assert(ni_>=0 && nj>=0); }
-
-   Array2(int ni_, int nj_, T* data_, size_type max_n_)
-      : ni(ni_), nj(nj_), a(ni_*nj_, data_, max_n_)
-   { assert(ni_>=0 && nj>=0); }
-
-   template<class OtherArrayT>
-   Array2(Array2<T, OtherArrayT>& other)
-      : ni(other.ni), nj(other.nj), a(other.a)
-   {}
-
-   ~Array2(void)
-   {
-#ifndef NDEBUG
-      ni=nj=0;
-#endif
-   }
-
-   const T& operator()(int i, int j) const
-   {
-      assert(i>=0 && i<ni && j>=0 && j<nj);
-      return a[i+ni*j];
-   }
-
-   T& operator()(int i, int j)
-   {
-      assert(i>=0 && i<ni && j>=0 && j<nj);
-      return a[i+ni*j];
-   }
-
-   bool operator==(const Array2<T>& x) const
-   { return ni==x.ni && nj==x.nj && a==x.a; }
-
-   bool operator!=(const Array2<T>& x) const
-   { return ni!=x.ni || nj!=x.nj || a!=x.a; }
-
-   bool operator<(const Array2<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;
-      return a<x.a;
-   }
-
-   bool operator>(const Array2<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;
-      return a>x.a;
-   }
-
-   bool operator<=(const Array2<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;
-      return a<=x.a;
-   }
-
-   bool operator>=(const Array2<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;
-      return a>=x.a;
-   }
-
-   void assign(const T& value)
-   { a.assign(value); }
-
-   void assign(int ni_, int nj_, const T& value)
-   {
-      a.assign(ni_*nj_, value);
-      ni=ni_;
-      nj=nj_;
-   }
-    
-   void assign(int ni_, int nj_, const T* copydata)
-   {
-      a.assign(ni_*nj_, copydata);
-      ni=ni_;
-      nj=nj_;
-   }
-    
-   const T& at(int i, int j) const
-   {
-      assert(i>=0 && i<ni && j>=0 && j<nj);
-      return a[i+ni*j];
-   }
-
-   T& at(int i, int j)
-   {
-      assert(i>=0 && i<ni && j>=0 && j<nj);
-      return a[i+ni*j];
-   }
-
-   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=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_, const T& value)
-   {
-      a.fill(ni_*nj_, value);
-      ni=ni_;
-      nj=nj_;
-   }
-    
-   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)
-   { a.reserve(reserve_ni*reserve_nj); }
-
-   void resize(int ni_, int nj_)
-   {
-      assert(ni_>=0 && nj_>=0);
-      a.resize(ni_*nj_);
-      ni=ni_;
-      nj=nj_;
-   }
-
-   void resize(int ni_, int nj_, const T& value)
-   {
-      assert(ni_>=0 && nj_>=0);
-      a.resize(ni_*nj_, value);
-      ni=ni_;
-      nj=nj_;
-   }
-
-   void set_zero(void)
-   { a.set_zero(); }
-
-   size_type size(void) const
-   { return a.size(); }
-
-   void swap(Array2<T>& x)
-   {
-      std::swap(ni, x.ni);
-      std::swap(nj, x.nj);
-      a.swap(x.a);
-   }
-
-   void trim(void)
-   { a.trim(); }
-};
-
-// some common arrays
-
-typedef Array2<double, Array1<double> >                         Array2d;
-typedef Array2<float, Array1<float> >                           Array2f;
-typedef Array2<long long, Array1<long long> >                   Array2ll;
-typedef Array2<unsigned long long, Array1<unsigned long long> > Array2ull;
-typedef Array2<int, Array1<int> >                               Array2i;
-typedef Array2<unsigned int, Array1<unsigned int> >             Array2ui;
-typedef Array2<short, Array1<short> >                           Array2s;
-typedef Array2<unsigned short, Array1<unsigned short> >         Array2us;
-typedef Array2<char, Array1<char> >                             Array2c;
-typedef Array2<unsigned char, Array1<unsigned char> >           Array2uc;
-
-// and wrapped versions
-
-typedef Array2<double, WrapArray1<double> >                         WrapArray2d;
-typedef Array2<float, WrapArray1<float> >                           WrapArray2f;
-typedef Array2<long long, WrapArray1<long long> >                   WrapArray2ll;
-typedef Array2<unsigned long long, WrapArray1<unsigned long long> > WrapArray2ull;
-typedef Array2<int, WrapArray1<int> >                               WrapArray2i;
-typedef Array2<unsigned int, WrapArray1<unsigned int> >             WrapArray2ui;
-typedef Array2<short, WrapArray1<short> >                           WrapArray2s;
-typedef Array2<unsigned short, WrapArray1<unsigned short> >         WrapArray2us;
-typedef Array2<char, WrapArray1<char> >                             WrapArray2c;
-typedef Array2<unsigned char, WrapArray1<unsigned char> >           WrapArray2uc;
-
-}
-
-#endif
diff --git a/extern/softbody/src/sdfgen/array3.h b/extern/softbody/src/sdfgen/array3.h
deleted file mode 100644
index 6334d24df69..00000000000
--- a/extern/softbody/src/sdfgen/array3.h
+++ /dev/null
@@ -1,276 +0,0 @@
-#ifndef ARRAY3_H
-#define ARRAY3_H
-
-#include "array1.h"
-#include <algorithm>
-#include <cassert>
-#include <vector>
-
-namespace sdfgen {
-
-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/extern/softbody/src/sdfgen/config.h.in b/extern/softbody/src/sdfgen/config.h.in
deleted file mode 100644
index a1203207b6f..00000000000
--- a/extern/softbody/src/sdfgen/config.h.in
+++ /dev/null
@@ -1,11 +0,0 @@
-#ifndef SDFGEN_CONFIG_H
-#define SDFGEN_CONFIG_H
-
-/*
- * Configuration Header for SDFGEN
- */
- 
-/// Configured libraries
-#cmakedefine HAVE_VTK
-
-#endif //SDFGEN_CONFIG_H
diff --git a/extern/softbody/src/sdfgen/hashgrid.h b/extern/softbody/src/sdfgen/hashgrid.h
deleted file mode 100644
index a5eb9aaae8f..00000000000
--- a/extern/softbody/src/sdfgen/hashgrid.h
+++ /dev/null
@@ -1,141 +0,0 @@
-#ifndef HASHGRID_H
-#define HASHGRID_H
-
-#include "hashtable.h"
-#include "vec.h"
-
-//========================================================= first do 2D ============================
-
-namespace sdfgen {
-
-template<class DataType>
-struct HashGrid2
-{
-   double dx, overdx; // side-length of a grid cell and its reciprocal
-   HashTable<Vec2i,DataType> grid;
-
-   explicit HashGrid2(double dx_=1, int expected_size=512)
-      : dx(dx_), overdx(1/dx_), grid(expected_size)
-   {}
-
-   // only do this with an empty grid
-   void set_grid_size(double dx_)
-   { assert(size()==0); dx=dx_; overdx=1/dx; }
-
-   void add_point(const Vec2d &x, const DataType &datum)
-   { grid.add(round(x*overdx), datum); }
-
-   void delete_point(const Vec2d &x, const DataType &datum)
-   { grid.delete_entry(round(x*overdx), datum); }
-
-   void add_box(const Vec2d &xmin, const Vec2d &xmax, const DataType &datum)
-   {
-      Vec2i imin=round(xmin*overdx), imax=round(xmax*overdx);
-      for(int j=imin[1]; j<=imax[1]; ++j) for(int i=imin[0]; i<=imax[0]; ++i)
-         grid.add(Vec2i(i,j), datum);
-   }
-
-   void delete_box(const Vec2d &xmin, const Vec2d &xmax, const DataType &datum)
-   {
-      Vec2i imin=round(xmin*overdx), imax=round(xmax*overdx);
-      for(int j=imin[1]; j<=imax[1]; ++j) for(int i=imin[0]; i<=imax[0]; ++i)
-         grid.delete_entry(Vec2i(i,j), datum);
-   }
-
-   unsigned int size(void) const
-   { return grid.size(); }
-
-   void clear(void)
-   { grid.clear(); }
-
-   void reserve(unsigned int expected_size)
-   { grid.reserve(expected_size); }
-
-   bool find_first_point(const Vec2d &x, DataType &datum) const
-   { return grid.get_entry(round(x*overdx), datum); }
-
-   bool find_point(const Vec2d &x, std::vector<DataType> &data_list) const
-   {
-      data_list.resize(0);
-      grid.append_all_entries(round(x*overdx), data_list);
-      return data_list.size()>0;
-   }
-
-   bool find_box(const Vec2d &xmin, const Vec2d &xmax, std::vector<DataType> &data_list) const
-   {
-      data_list.resize(0);
-      Vec2i imin=round(xmin*overdx), imax=round(xmax*overdx);
-      for(int j=imin[1]; j<=imax[1]; ++j) for(int i=imin[0]; i<=imax[0]; ++i)
-         grid.append_all_entries(Vec2i(i,j), data_list);
-      return data_list.size()>0;
-   }
-};
-
-//==================================== and now in 3D =================================================
-
-template<class DataType>
-struct HashGrid3
-{
-   double dx, overdx; // side-length of a grid cell and its reciprocal
-   HashTable<Vec3i,DataType> grid;
-
-   explicit HashGrid3(double dx_=1, int expected_size=512)
-      : dx(dx_), overdx(1/dx_), grid(expected_size)
-   {}
-
-   // only do this with an empty grid
-   void set_grid_size(double dx_)
-   { assert(size()==0); dx=dx_; overdx=1/dx; }
-
-   void add_point(const Vec3d &x, const DataType &datum)
-   { grid.add(round(x*overdx), datum); }
-
-   void delete_point(const Vec3d &x, const DataType &datum)
-   { grid.delete_entry(round(x*overdx), datum); }
-
-   void add_box(const Vec3d &xmin, const Vec3d &xmax, const DataType &datum)
-   {
-      Vec3i imin=round(xmin*overdx), imax=round(xmax*overdx);
-      for(int k=imin[2]; k<=imax[2]; ++k) for(int j=imin[1]; j<=imax[1]; ++j) for(int i=imin[0]; i<=imax[0]; ++i)
-         grid.add(Vec3i(i,j,k), datum);
-   }
-
-   void delete_box(const Vec3d &xmin, const Vec3d &xmax, const DataType &datum)
-   {
-      Vec3i imin=round(xmin*overdx), imax=round(xmax*overdx);
-      for(int k=imin[2]; k<=imax[2]; ++k) for(int j=imin[1]; j<=imax[1]; ++j) for(int i=imin[0]; i<=imax[0]; ++i)
-         grid.delete_entry(Vec3i(i,j,k), datum);
-   }
-   
-   unsigned int size(void) const
-   { return grid.size(); }
-
-   void clear(void)
-   { grid.clear(); }
-
-   void reserve(unsigned int expected_size)
-   { grid.reserve(expected_size); }
-
-   bool find_first_point(const Vec3d &x, DataType &index) const
-   { return grid.get_entry(round(x*overdx), index); }
-
-   bool find_point(const Vec3d &x, std::vector<DataType> &data_list) const
-   {
-      data_list.resize(0);
-      grid.append_all_entries(round(x*overdx), data_list);
-      return data_list.size()>0;
-   }
-
-   bool find_box(const Vec3d &xmin, const Vec3d &xmax, std::vector<DataType> &data_list) const
-   {
-      data_list.resize(0);
-      Vec3i imin=round(xmin*overdx), imax=round(xmax*overdx);
-      for(int k=imin[2]; k<=imax[2]; ++k) for(int j=imin[1]; j<=imax[1]; ++j) for(int i=imin[0]; i<=imax[0]; ++i)
-         grid.append_all_entries(Vec3i(i, j, k), data_list);
-      return data_list.size()>0;
-   }
-};
-
-}
-
-#endif
diff --git a/extern/softbody/src/sdfgen/hashtable.h b/extern/softbody/src/sdfgen/hashtable.h
deleted file mode 100644
index c81fb59ea91..00000000000
--- a/extern/softbody/src/sdfgen/hashtable.h
+++ /dev/null
@@ -1,252 +0,0 @@
-#ifndef HASHTABLE_H
-#define HASHTABLE_H
-
-#include <functional>
-#include <iostream>
-#include <vector>
-
-namespace sdfgen {
-
-template<class Key, class Data>
-struct HashEntry
-{
-   Key key;
-   Data data;
-   int next;
-};
-
-// a useful core hash function
-inline unsigned int hash(unsigned int k)
-{ return k*2654435769u; }
-
-// default hash function object
-struct DefaultHashFunction
-{
-   template<typename Key>
-   unsigned int operator() (const Key &k) const { return hash(k); }
-};
-
-struct equal
-{
-   template<typename T>
-   bool operator() (const T &a, const T &b) const { return a==b; }
-};
-
-template<typename Key, typename Data, class HashFunction=DefaultHashFunction, class KeyEqual=equal>
-struct HashTable
-{
-   unsigned int table_rank;
-   unsigned int table_bits;
-   std::vector<int> table;
-   unsigned int num_entries;
-   std::vector<HashEntry<Key, Data> > pool;
-   int free_list;
-   const HashFunction hash_function;
-   const KeyEqual key_equal;
-
-   explicit HashTable(unsigned int expected_size=64)
-      : hash_function(HashFunction()), key_equal(KeyEqual())
-   { init(expected_size); }
-
-   explicit HashTable(const HashFunction &hf, unsigned int expected_size=64)
-      : hash_function(hf), key_equal(KeyEqual())
-   { init(expected_size); }
-
-   void init(unsigned int expected_size)
-   {
-      unsigned int i;
-      num_entries=0;
-      table_rank=4;
-      while(1u<<table_rank < expected_size)
-         ++table_rank;
-      ++table_rank; // give us some extra room
-      table_bits=(1u<<table_rank)-1;
-      table.resize(1u<<table_rank);
-      for(i=0; i<table.size(); ++i)
-         table[i]=-1; // empty list
-      pool.resize(1u<<table_rank);
-      free_list=0;
-      for(unsigned int i=0; i<pool.size()-1; ++i)
-         pool[i].next=i+1;
-      pool[pool.size()-1].next=-1; // end of free list
-   }
-
-   void add(const Key &k, const Data &d)
-   {
-      if(free_list==-1)
-         reserve(1u<<(table_rank+1));
-      int i=free_list; // where we're going to put the new entry
-      free_list=pool[i].next;
-      unsigned int t=hash_function(k)&table_bits; // index into table
-      pool[i].key=k;
-      pool[i].data=d;
-      pool[i].next=table[t]; // put the new entry at the start of table[t]'s list
-      table[t]=i;
-      ++num_entries;
-   }
-
-   void delete_entry(const Key &k, const Data &d) // delete first entry that matches both key and data
-   {
-      unsigned int t=hash_function(k)&table_bits;
-      int i=table[t], *p_i=&table[t];
-      while(i!=-1){
-         if(key_equal(k, pool[i].key) && d==pool[i].data){
-            *p_i=pool[i].next; // make list skip over this entry
-            pool[i].next=free_list; // and put it on the front of the free list
-            free_list=i;
-            return; // and we're done
-         }
-         p_i=&pool[i].next;
-         i=*p_i;
-      }
-   }
-
-   unsigned int size() const
-   { return num_entries; }
-
-   void clear()
-   {
-      unsigned int i=0;
-      num_entries=0;
-      for(i=0; i<table.size(); ++i)
-         table[i]=-1; // empty list
-      free_list=0;
-      for(i=0; i<pool.size()-1; ++i)
-         pool[i].next=i+1;
-      pool[pool.size()-1].next=-1;
-   }
-
-   void reserve(unsigned int expected_size)
-   {
-      if(expected_size<=pool.size())
-         return;
-      while(1u<<table_rank < expected_size)
-         ++table_rank;
-      table_bits=(1u<<table_rank)-1;
-      // increase room for new entries
-      unsigned int old_size=(unsigned int)pool.size(), i;
-      pool.resize(1u<<table_rank);
-      for(i=old_size; i<pool.size()-1; ++i)
-         pool[i].next=i+1;
-      pool[i].next=free_list;
-      free_list=old_size;
-      // And finally need to redo table (rehash entries)
-      old_size=(unsigned int)table.size();
-      table.resize(1u<<table_rank);
-      unsigned int t;
-      for(t=old_size; t<table.size(); ++t)
-         table[t]=-1; // initially make new lists empty
-      int j, *p_j;
-      for(t=0; t<old_size; ++t){
-         j=table[t]; 
-         p_j=&table[t];
-         while(j!=-1){
-            unsigned int new_t=hash_function(pool[j].key)&table_bits;
-            if(new_t!=t){ // j doesn't belong in this list anymore?
-               // delete from this list
-               *p_j=pool[j].next;
-               // add to correct list
-               pool[j].next=table[new_t];
-               table[new_t]=j;
-            }else
-               p_j=&(pool[j].next);
-            j=*p_j;
-         }
-      }
-   }
-
-   bool has_entry(const Key &k) const
-   {
-      unsigned int t=hash_function(k)&table_bits;
-      int i=table[t];
-      while(i!=-1){
-         if(key_equal(k, pool[i].key))
-            return true;
-         i=pool[i].next;
-      }
-      return false;
-   }
-
-   bool get_entry(const Key &k, Data &data_return) const
-   {
-      unsigned int t=hash_function(k)&table_bits;
-      int i=table[t];
-      while(i!=-1){
-         if(key_equal(k, pool[i].key)){
-            data_return=pool[i].data;
-            return true;
-         }
-         i=pool[i].next;
-      }
-      return false;
-   }
-
-   void append_all_entries(const Key& k, std::vector<Data>& data_return) const
-   {
-      unsigned int t=hash_function(k)&table_bits;
-      int i=table[t];
-      while(i!=-1){
-         if(key_equal(k, pool[i].key)) data_return.push_back(pool[i].data);
-         i=pool[i].next;
-      }
-   }
-
-   Data &operator() (const Key &k, const Data &missing_data)
-   {
-      unsigned int t=hash_function(k)&table_bits;
-      int i=table[t];
-      while(i!=-1){
-         if(key_equal(k, pool[i].key))
-            return pool[i].data;
-         i=pool[i].next;
-      }
-      add(k, missing_data); // note - this could cause the table to be resized, and t made out-of-date
-      return pool[table[hash_function(k)&table_bits]].data; // we know that add() puts it here!
-   }
-
-   const Data &operator() (const Key &k, const Data &missing_data) const
-   {
-      unsigned int t=hash_function(k)&table_bits;
-      int i=table[t];
-      while(i!=-1){
-         if(key_equal(k, pool[i].key))
-            return pool[i].data;
-         i=pool[i].next;
-      }
-      return missing_data;
-   }
-
-   void output_statistics() const
-   {
-      std::vector<int> lengthcount(table.size());
-      unsigned int t;
-      int total=0;
-      for(t=0; t<table.size(); ++t){
-         int i=table[t], length=0;
-         while(i!=-1){
-            ++length;
-            i=pool[i].next;
-         }
-         ++lengthcount[length];
-         ++total;
-      }
-      int subtotal=0;
-      int maxlength=0;
-      for(t=0; t<lengthcount.size() && t<10; ++t){
-         subtotal+=lengthcount[t];
-         if(lengthcount[t]>0){
-            std::cout<<"length "<<t<<": "<<lengthcount[t]<<"   ("<<lengthcount[t]/(float)total*100.0<<"%)"<<std::endl;
-            maxlength=t;
-         }
-      }
-      std::cout<<"rest: "<<total-subtotal<<"   ("<<100.0*(1.0-subtotal/(float)total)<<"%)"<<std::endl;
-      for(; t<lengthcount.size(); ++t)
-         if(lengthcount[t]>0)
-            maxlength=t;
-      std::cout<<"longest list: "<<maxlength<<std::endl;
-   }
-};
-
-}
-
-#endif
diff --git a/extern/softbody/src/sdfgen/main.cpp b/extern/softbody/src/sdfgen/main.cpp
deleted file mode 100644
index 4b5a6f12f59..00000000000
--- a/extern/softbody/src/sdfgen/main.cpp
+++ /dev/null
@@ -1,180 +0,0 @@
-//SDFGen - A simple grid-based signed distance field (level set) generator for triangle meshes.
-//Written by Christopher Batty (christopherbatty@yahoo.com, www.cs.columbia.edu/~batty)
-//...primarily using code from Robert Bridson's website (www.cs.ubc.ca/~rbridson)
-//This code is public domain. Feel free to mess with it, let me know if you like it.
-
-#include "makelevelset3.h"
-#include "config.h"
-
-#ifdef HAVE_VTK
-  #include <vtkImageData.h>
-  #include <vtkFloatArray.h>
-  #include <vtkXMLImageDataWriter.h>
-  #include <vtkPointData.h>
-  #include <vtkSmartPointer.h>
-#endif
-
-
-#include <fstream>
-#include <iostream>
-#include <sstream>
-#include <limits>
-
-int main(int argc, char* argv[]) {
-  
-  if(argc != 4) {
-    std::cout << "SDFGen - A utility for converting closed oriented triangle meshes into grid-based signed distance fields.\n";
-    std::cout << "\nThe output file format is:";
-    std::cout << "<ni> <nj> <nk>\n";
-    std::cout << "<origin_x> <origin_y> <origin_z>\n";
-    std::cout << "<dx>\n";
-    std::cout << "<value_1> <value_2> <value_3> [...]\n\n";
-    
-    std::cout << "(ni,nj,nk) are the integer dimensions of the resulting distance field.\n";
-    std::cout << "(origin_x,origin_y,origin_z) is the 3D position of the grid origin.\n";
-    std::cout << "<dx> is the grid spacing.\n\n";
-    std::cout << "<value_n> are the signed distance data values, in ascending order of i, then j, then k.\n";
-
-    std::cout << "The output filename will match that of the input, with the OBJ suffix replaced with SDF.\n\n";
-
-    std::cout << "Usage: SDFGen <filename> <dx> <padding>\n\n";
-    std::cout << "Where:\n";
-    std::cout << "\t<filename> specifies a Wavefront OBJ (text) file representing a *triangle* mesh (no quad or poly meshes allowed). File must use the suffix \".obj\".\n";
-    std::cout << "\t<dx> specifies the length of grid cell in the resulting distance field.\n";
-    std::cout << "\t<padding> specifies the number of cells worth of padding between the object bound box and the boundary of the distance field grid. Minimum is 1.\n\n";
-    
-    exit(-1);
-  }
-
-  std::string filename(argv[1]);
-  if(filename.size() < 5 || filename.substr(filename.size()-4) != std::string(".obj")) {
-    std::cerr << "Error: Expected OBJ file with filename of the form <name>.obj.\n";
-    exit(-1);
-  }
-
-  std::stringstream arg2(argv[2]);
-  float dx;
-  arg2 >> dx;
-  
-  std::stringstream arg3(argv[3]);
-  int padding;
-  arg3 >> padding;
-
-  if(padding < 1) padding = 1;
-  //start with a massive inside out bound box.
-  Vec3f min_box(std::numeric_limits<float>::max(),std::numeric_limits<float>::max(),std::numeric_limits<float>::max()), 
-    max_box(-std::numeric_limits<float>::max(),-std::numeric_limits<float>::max(),-std::numeric_limits<float>::max());
-  
-  std::cout << "Reading data.\n";
-
-  std::ifstream infile(argv[1]);
-  if(!infile) {
-    std::cerr << "Failed to open. Terminating.\n";
-    exit(-1);
-  }
-
-  int ignored_lines = 0;
-  std::string line;
-  std::vector<Vec3f> vertList;
-  std::vector<Vec3ui> faceList;
-  while(!infile.eof()) {
-    std::getline(infile, line);
-
-    //.obj files sometimes contain vertex normals indicated by "vn"
-    if(line.substr(0,1) == std::string("v") && line.substr(0,2) != std::string("vn")){
-      std::stringstream data(line);
-      char c;
-      Vec3f point;
-      data >> c >> point[0] >> point[1] >> point[2];
-      vertList.push_back(point);
-      update_minmax(point, min_box, max_box);
-    }
-    else if(line.substr(0,1) == std::string("f")) {
-      std::stringstream data(line);
-      char c;
-      int v0,v1,v2;
-      data >> c >> v0 >> v1 >> v2;
-      faceList.push_back(Vec3ui(v0-1,v1-1,v2-1));
-    }
-    else if( line.substr(0,2) == std::string("vn") ){
-      std::cerr << "Obj-loader is not able to parse vertex normals, please strip them from the input file. \n";
-      exit(-2); 
-    }
-    else {
-      ++ignored_lines; 
-    }
-  }
-  infile.close();
-  
-  if(ignored_lines > 0)
-    std::cout << "Warning: " << ignored_lines << " lines were ignored since they did not contain faces or vertices.\n";
-
-  std::cout << "Read in " << vertList.size() << " vertices and " << faceList.size() << " faces." << std::endl;
-
-  //Add padding around the box.
-  Vec3f unit(1,1,1);
-  min_box -= padding*dx*unit;
-  max_box += padding*dx*unit;
-  Vec3ui sizes = Vec3ui((max_box - min_box)/dx);
-  
-  std::cout << "Bound box size: (" << min_box << ") to (" << max_box << ") with dimensions " << sizes << "." << std::endl;
-
-  std::cout << "Computing signed distance field.\n";
-  Array3f phi_grid;
-  make_level_set3(faceList, vertList, min_box, dx, sizes[0], sizes[1], sizes[2], phi_grid);
-
-  std::string outname;
-
-  #ifdef HAVE_VTK
-    // If compiled with VTK, we can directly output a volumetric image format instead
-    //Very hackily strip off file suffix.
-    outname = filename.substr(0, filename.size()-4) + std::string(".vti");
-    std::cout << "Writing results to: " << outname << "\n";
-    vtkSmartPointer<vtkImageData> output_volume = vtkSmartPointer<vtkImageData>::New();
-
-    output_volume->SetDimensions(phi_grid.ni ,phi_grid.nj ,phi_grid.nk);
-    output_volume->SetOrigin( phi_grid.ni*dx/2, phi_grid.nj*dx/2,phi_grid.nk*dx/2);
-    output_volume->SetSpacing(dx,dx,dx);
-
-    vtkSmartPointer<vtkFloatArray> distance = vtkSmartPointer<vtkFloatArray>::New();
-    
-    distance->SetNumberOfTuples(phi_grid.a.size());
-    
-    output_volume->GetPointData()->AddArray(distance);
-    distance->SetName("Distance");
-
-    for(unsigned int i = 0; i < phi_grid.a.size(); ++i) {
-      distance->SetValue(i, phi_grid.a[i]);
-    }
-
-    vtkSmartPointer<vtkXMLImageDataWriter> writer =
-    vtkSmartPointer<vtkXMLImageDataWriter>::New();
-    writer->SetFileName(outname.c_str());
-
-    #if VTK_MAJOR_VERSION <= 5
-      writer->SetInput(output_volume);
-    #else
-      writer->SetInputData(output_volume);
-    #endif
-    writer->Write();
-
-  #else
-    // if VTK support is missing, default back to the original ascii file-dump.
-    //Very hackily strip off file suffix.
-    outname = filename.substr(0, filename.size()-4) + std::string(".sdf");
-    std::cout << "Writing results to: " << outname << "\n";
-    
-    std::ofstream outfile( outname.c_str());
-    outfile << phi_grid.ni << " " << phi_grid.nj << " " << phi_grid.nk << std::endl;
-    outfile << min_box[0] << " " << min_box[1] << " " << min_box[2] << std::endl;
-    outfile << dx << std::endl;
-    for(unsigned int i = 0; i < phi_grid.a.size(); ++i) {
-      outfile << phi_grid.a[i] << std::endl;
-    }
-    outfile.close();
-  #endif
-
-  std::cout << "Processing complete.\n";
-
-return 0;
-}
diff --git a/extern/softbody/src/sdfgen/makelevelset3.cpp b/extern/softbody/src/sdfgen/makelevelset3.cpp
deleted file mode 100644
index da2043ea6ae..00000000000
--- a/extern/softbody/src/sdfgen/makelevelset3.cpp
+++ /dev/null
@@ -1,187 +0,0 @@
-#include "makelevelset3.h"
-
-namespace sdfgen {
-
-// 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/extern/softbody/src/sdfgen/makelevelset3.h b/extern/softbody/src/sdfgen/makelevelset3.h
deleted file mode 100644
index efc3a67da7b..00000000000
--- a/extern/softbody/src/sdfgen/makelevelset3.h
+++ /dev/null
@@ -1,20 +0,0 @@
-#ifndef MAKELEVELSET3_H
-#define MAKELEVELSET3_H
-
-#include "array3.h"
-#include "vec.h"
-
-namespace sdfgen {
-
-// 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/extern/softbody/src/sdfgen/util.h b/extern/softbody/src/sdfgen/util.h
deleted file mode 100644
index 8fb9b23f344..00000000000
--- a/extern/softbody/src/sdfgen/util.h
+++ /dev/null
@@ -1,474 +0,0 @@
-#ifndef UTIL_H
-#define UTIL_H
-
-#include <algorithm>
-#include <vector>
-#include <cmath>
-#include <iostream>
-
-namespace sdfgen {
-
-#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/extern/softbody/src/sdfgen/vec.h b/extern/softbody/src/sdfgen/vec.h
deleted file mode 100644
index 14ed214a1e1..00000000000
--- a/extern/softbody/src/sdfgen/vec.h
+++ /dev/null
@@ -1,482 +0,0 @@
-#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.)
-
-namespace sdfgen {
-
-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