diff --git a/CMakeLists.txt b/CMakeLists.txt
index 011f5ab4..9e171f9c 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,7 +20,7 @@ ENDIF()
 # Set the project name
 PROJECT( LBPM-WIA )
 IF ( NOT CXX_STD )
-    SET( CXX_STD 98 )
+    SET( CXX_STD 11 )
 ENDIF()
 
 
diff --git a/common/Array.h b/common/Array.h
index a0c2f214..f926b101 100644
--- a/common/Array.h
+++ b/common/Array.h
@@ -269,6 +269,12 @@ public:
      */
     void scale( const TYPE &scale );
 
+    /*!
+     * Set the values of this array to pow(base, exp)
+     * @param base        Base array
+     * @param exp         Exponent value
+     */
+    void pow( const Array<TYPE> &baseArray, const TYPE &exp );
 
     //! Destructor
     ~Array();
@@ -283,7 +289,7 @@ public:
 
 
     //! Return the size of the Array
-    std::vector<size_t> size() const;
+    inline std::vector<size_t> size() const { return std::vector<size_t>( d_N, d_N + d_ndim ); }
 
 
     //! Return the size of the Array
@@ -356,6 +362,13 @@ public:
     template <class TYPE2>
     void copySubset( const std::vector<size_t> &index, const Array<TYPE2> &subset );
 
+    /*!
+     * Add data from an array into a subset of this array
+     * @param index         Index of the subset (imin,imax,jmin,jmax,kmin,kmax,...)
+     * @param subset        The subset array to add from
+     */
+    void addSubset( const std::vector<size_t> &index, const Array<TYPE> &subset );
+
 
     /*!
      * Access the desired element
@@ -526,6 +539,9 @@ public:
     //! Coarsen an array using the given filter
     Array<TYPE> coarsen( const Array<TYPE>& filter ) const;
 
+    //! Coarsen an array using the given filter
+    Array<TYPE> coarsen( const std::vector<size_t>& ratio, std::function<TYPE(const Array<TYPE>&)> filter ) const;
+
 private:
     int d_ndim;                  // Number of dimensions in array
     size_t d_N[ARRAY_NDIM_MAX];  // Size of each dimension
diff --git a/common/Array.hpp b/common/Array.hpp
index 6fa7b845..ca2482df 100644
--- a/common/Array.hpp
+++ b/common/Array.hpp
@@ -147,18 +147,6 @@ void Array<TYPE>::clear()
 }
 
 
-/********************************************************
-*  Return the size of the array                         *
-********************************************************/
-template <class TYPE>
-std::vector<size_t> Array<TYPE>::size() const
-{
-    if ( d_ndim==0 )
-        return std::vector<size_t>();
-    return std::vector<size_t>( d_N, d_N + d_ndim );
-}
-
-
 /********************************************************
 *  Check if the size of the array matches rhs           *
 ********************************************************/
@@ -416,6 +404,34 @@ void Array<TYPE>::copySubset( const std::vector<size_t> &index, const Array<TYPE
         }
     }
 }
+
+template <class TYPE>
+void Array<TYPE>::addSubset( const std::vector<size_t> &index, const Array<TYPE> &subset )
+{
+    // Get the subset indices
+    checkSubsetIndex( index );
+    std::array<size_t,5> first, last, N1;
+    getSubsetArrays( index, first, last, N1 );
+    std::array<size_t,5> N2 = getDimArray();
+    // add the sub-array
+    #if ARRAY_NDIM_MAX > 5
+        #error Function programmed for more than 5 dimensions
+    #endif
+    for (size_t i4=first[4]; i4<=last[4]; i4++) {
+        for (size_t i3=first[3]; i3<=last[3]; i3++) {
+            for (size_t i2=first[2]; i2<=last[2]; i2++) {
+                for (size_t i1=first[1]; i1<=last[1]; i1++) {
+                    for (size_t i0=first[0]; i0<=last[0]; i0++) {
+                        size_t k1 = GET_ARRAY_INDEX5D( N1, i0-first[0], 
+                            i1-first[1], i2-first[2], i3-first[3], i4-first[4] );
+                        size_t k2  = GET_ARRAY_INDEX5D( N2, i0, i1, i2, i3, i4 );
+                        d_data[k2] += subset.d_data[k1];
+                    }
+                }
+            }
+        }
+    }
+}
 // clang-format on
 
 
@@ -613,7 +629,17 @@ void Array<TYPE>::scale( const TYPE &value )
     for ( size_t i = 0; i < d_length; i++ )
         d_data[i] *= value;
 }
+template <class TYPE>
+    void Array<TYPE>::pow(const Array<TYPE> &baseArray, const TYPE &exp )
+{
+    // not insisting on the shapes being the same
+    // but insisting on the total size being the same
+    AMP_ASSERT(d_length==baseArray.length());
 
+    const auto base_data = baseArray.data();
+    for ( size_t i = 0; i < d_length; i++ )
+        d_data[i]  = std::pow(base_data[i], exp);
+}
 
 /********************************************************
 *  Simple math operations                               *
@@ -996,6 +1022,35 @@ Array<TYPE> Array<TYPE>::coarsen( const Array<TYPE>& filter ) const
     }
     return y;
 }
+template <class TYPE>
+Array<TYPE> Array<TYPE>::coarsen( const std::vector<size_t>& ratio, std::function<TYPE(const Array<TYPE>&)> filter ) const
+{
+    ASSERT((int)ratio.size()==d_ndim);
+    auto S2 = size();
+    for (size_t i=0; i<S2.size(); i++) {
+        S2[i] /= ratio[i];
+        INSIST(S2[i]*ratio[i]==size(i),"Array must be multiple of filter size");
+    }
+    Array<TYPE> tmp(ratio);
+    TYPE* tmp2 = tmp.data();
+    Array<TYPE> y( S2 );
+    INSIST(d_ndim<=3,"Function programmed for more than 3 dimensions");
+    for (size_t k1=0; k1<y.d_N[2]; k1++) {
+        for (size_t j1=0; j1<y.d_N[1]; j1++) {
+            for (size_t i1=0; i1<y.d_N[0]; i1++) {
+                for (size_t k2=0; k2<ratio[2]; k2++) {
+                    for (size_t j2=0; j2<ratio[1]; j2++) {
+                        for (size_t i2=0; i2<ratio[0]; i2++) {
+                            tmp2[GET_ARRAY_INDEX3D(tmp.d_N,i2,j2,k2)] = this->operator()(i1*ratio[0]+i2,j1*ratio[1]+j2,k1*ratio[2]+k2);
+                        }
+                    }
+                }
+                y(i1,j1,k1) = filter(tmp);
+            }
+        }
+    }
+    return y;
+}
 
 
 #endif
diff --git a/common/MPI_Helpers.h b/common/MPI_Helpers.h
index fde70618..5a02397c 100644
--- a/common/MPI_Helpers.h
+++ b/common/MPI_Helpers.h
@@ -169,6 +169,44 @@ void unpack( std::set<TYPE>& data, const char *buffer );
 
 
 
+// Helper functions
+inline float sumReduce( MPI_Comm comm, float x )
+{
+    float y = 0;
+	MPI_Allreduce(&x,&y,1,MPI_FLOAT,MPI_SUM,comm);
+    return y;
+}
+inline int sumReduce( MPI_Comm comm, int x )
+{
+    int y = 0;
+	MPI_Allreduce(&x,&y,1,MPI_INT,MPI_SUM,comm);
+    return y;
+}
+inline bool sumReduce( MPI_Comm comm, bool x )
+{
+    int y = sumReduce( comm, x?1:0 );
+    return y>0;
+}
+inline std::vector<float> sumReduce( MPI_Comm comm, const std::vector<float>& x )
+{
+    auto y = x;
+	MPI_Allreduce(x.data(),y.data(),x.size(),MPI_FLOAT,MPI_SUM,comm);
+    return y;
+}
+inline std::vector<int> sumReduce( MPI_Comm comm, const std::vector<int>& x )
+{
+    auto y = x;
+	MPI_Allreduce(x.data(),y.data(),x.size(),MPI_INT,MPI_SUM,comm);
+    return y;
+}
+inline float maxReduce( MPI_Comm comm, float x )
+{
+    float y = 0;
+	MPI_Allreduce(&x,&y,1,MPI_FLOAT,MPI_MAX,comm);
+    return y;
+}
+
+
 #endif
 
 
diff --git a/tests/lbpm_morph_pp.cpp b/tests/lbpm_morph_pp.cpp
index 503282ab..18da6292 100644
--- a/tests/lbpm_morph_pp.cpp
+++ b/tests/lbpm_morph_pp.cpp
@@ -133,7 +133,7 @@ int main(int argc, char **argv)
 	sprintf(LocalRankFilename,"SignDist.%05i",rank);
 	FILE *DIST = fopen(LocalRankFilename,"rb");
 	size_t ReadSignDist;
-	ReadSignDist=fread(SignDist.get(),8,N,DIST);
+	ReadSignDist=fread(SignDist.data(),8,N,DIST);
 	if (ReadSignDist != size_t(N)) printf("lbpm_morphdrain_pp: Error reading signed distance function (rank=%i)\n",rank);
 	fclose(DIST);
 
diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp
index 9365f5da..9db58f95 100644
--- a/tests/lbpm_uCT_pp.cpp
+++ b/tests/lbpm_uCT_pp.cpp
@@ -21,27 +21,17 @@
 #include "IO/Writer.h"
 #include "IO/netcdf.h"
 #include "analysis/analysis.h"
+#include "analysis/eikonal.h"
 
 #include "ProfilerApp.h"
 
 
-float sumReduce( MPI_Comm comm, float x )
+template<class T>
+inline int sign( T x )
 {
-    float y = 0;
-	MPI_Allreduce(&x,&y,1,MPI_FLOAT,MPI_SUM,comm);
-    return y;
-}
-int sumReduce( MPI_Comm comm, int x )
-{
-    int y = 0;
-	MPI_Allreduce(&x,&y,1,MPI_INT,MPI_SUM,comm);
-    return y;
-}
-float maxReduce( MPI_Comm comm, float x )
-{
-    float y = 0;
-	MPI_Allreduce(&x,&y,1,MPI_FLOAT,MPI_MAX,comm);
-    return y;
+    if ( x==0 )
+        return 0;
+    return x>0 ? 1:-1;
 }
 
 
@@ -172,160 +162,6 @@ inline void InterpolateMesh( const Array<float> &Coarse, Array<float> &Fine )
 	PROFILE_STOP("InterpolateMesh");
 }
 
-inline float minmod(float &a, float &b){
-
-	float value;
-
-	value = a;
-	if 	( a*b < 0.0)	    value=0.0;
-	else if (fabs(a) > fabs(b)) value = b;
-
-	return value;
-}
-
-
-inline float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps)
-{
-	PROFILE_START("Eikonal3D");
-
-	/*
-	 * This routine converts the data in the Distance array to a signed distance
-	 * by solving the equation df/dt = sign*(1-|grad f|), where Distance provides
-	 * the values of f on the mesh associated with domain Dm
-	 * It has been tested with segmented data initialized to values [-1,1]
-	 * and will converge toward the signed distance to the surface bounding the associated phases
-	 *
-	 * Reference:
-	 * Min C (2010) On reinitializing level set functions, Journal of Computational Physics	229
-	 */
-
-	int i,j,k;
-	float dt=0.1;
-	float Dx,Dy,Dz;
-	float Dxp,Dxm,Dyp,Dym,Dzp,Dzm;
-	float Dxxp,Dxxm,Dyyp,Dyym,Dzzp,Dzzm;
-	float sign,norm;
-	float LocalVar,GlobalVar,LocalMax,GlobalMax;
-
-	int xdim,ydim,zdim;
-	xdim=Dm.Nx-2;
-	ydim=Dm.Ny-2;
-	zdim=Dm.Nz-2;
-	fillHalo<float> fillData(Dm.Comm, Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
-
-	// Arrays to store the second derivatives
-	Array<float> Dxx(Dm.Nx,Dm.Ny,Dm.Nz);
-	Array<float> Dyy(Dm.Nx,Dm.Ny,Dm.Nz);
-	Array<float> Dzz(Dm.Nx,Dm.Ny,Dm.Nz);
-
-	int count = 0;
-	while (count < timesteps){
-
-		// Communicate the halo of values
-		fillData.fill(Distance);
-
-		// Compute second order derivatives
-		for (k=1;k<Dm.Nz-1;k++){
-			for (j=1;j<Dm.Ny-1;j++){
-				for (i=1;i<Dm.Nx-1;i++){
-					Dxx(i,j,k) = Distance(i+1,j,k) + Distance(i-1,j,k) - 2*Distance(i,j,k);
-					Dyy(i,j,k) = Distance(i,j+1,k) + Distance(i,j-1,k) - 2*Distance(i,j,k);
-					Dzz(i,j,k) = Distance(i,j,k+1) + Distance(i,j,k-1) - 2*Distance(i,j,k);
-				}
-			}
-		}
-		fillData.fill(Dxx);
-		fillData.fill(Dyy);
-		fillData.fill(Dzz);
-
-		LocalMax=LocalVar=0.0;
-		// Execute the next timestep
-        //  f(n+1) = f(n) + dt*sign(1-|grad f|)
-		for (k=1;k<Dm.Nz-1;k++){
-			for (j=1;j<Dm.Ny-1;j++){
-				for (i=1;i<Dm.Nx-1;i++){
-
-					int n = k*Dm.Nx*Dm.Ny + j*Dm.Nx + i;
-
-					sign = -1;
-					if (ID(i,j,k) == 1) sign = 1;
-
-					// local second derivative terms
-					Dxxp = minmod(Dxx(i,j,k),Dxx(i+1,j,k));
-					Dyyp = minmod(Dyy(i,j,k),Dyy(i,j+1,k));
-					Dzzp = minmod(Dzz(i,j,k),Dzz(i,j,k+1));
-					Dxxm = minmod(Dxx(i,j,k),Dxx(i-1,j,k));
-					Dyym = minmod(Dyy(i,j,k),Dyy(i,j-1,k));
-					Dzzm = minmod(Dzz(i,j,k),Dzz(i,j,k-1));
-
-					/* //............Compute upwind derivatives ...................
-                    Dxp = Distance(i+1,j,k) - Distance(i,j,k) + 0.5*Dxxp;
-                    Dyp = Distance(i,j+1,k) - Distance(i,j,k) + 0.5*Dyyp;
-                    Dzp = Distance(i,j,k+1) - Distance(i,j,k) + 0.5*Dzzp;
-
-                    Dxm = Distance(i,j,k) - Distance(i-1,j,k) + 0.5*Dxxm;
-                    Dym = Distance(i,j,k) - Distance(i,j-1,k) + 0.5*Dyym;
-                    Dzm = Distance(i,j,k) - Distance(i,j,k-1) + 0.5*Dzzm;
-					 */
-					Dxp = Distance(i+1,j,k);
-					Dyp = Distance(i,j+1,k);
-					Dzp = Distance(i,j,k+1);
-
-					Dxm = Distance(i-1,j,k);
-					Dym = Distance(i,j-1,k);
-					Dzm = Distance(i,j,k-1);
-
-					// Compute upwind derivatives for Godunov Hamiltonian
-					if (sign < 0.0){
-						if (Dxp > Dxm)  Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
-						else			Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
-
-						if (Dyp > Dym)  Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
-						else			Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
-
-						if (Dzp > Dzm)  Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
-						else			Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
-					}
-					else{
-						if (Dxp < Dxm)  Dx = Dxp - Distance(i,j,k) + 0.5*Dxxp;
-						else			Dx = Distance(i,j,k) - Dxm + 0.5*Dxxm;
-
-						if (Dyp < Dym)  Dy = Dyp - Distance(i,j,k) + 0.5*Dyyp;
-						else			Dy = Distance(i,j,k) - Dym + 0.5*Dyym;
-
-						if (Dzp < Dzm)  Dz = Dzp - Distance(i,j,k) + 0.5*Dzzp;
-						else			Dz = Distance(i,j,k) - Dzm + 0.5*Dzzm;
-					}
-
-					norm=sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
-					if (norm > 1.0) norm=1.0;
-
-					Distance(i,j,k) += dt*sign*(1.0 - norm);
-					LocalVar +=  dt*sign*(1.0 - norm);
-
-					if (fabs(dt*sign*(1.0 - norm)) > LocalMax)
-						LocalMax = fabs(dt*sign*(1.0 - norm));
-				}
-			}
-		}
-
-		MPI_Allreduce(&LocalVar,&GlobalVar,1,MPI_FLOAT,MPI_SUM,Dm.Comm);
-		MPI_Allreduce(&LocalMax,&GlobalMax,1,MPI_FLOAT,MPI_MAX,Dm.Comm);
-		GlobalVar /= (Dm.Nx-2)*(Dm.Ny-2)*(Dm.Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz;
-		count++;
-
-		if (count%50 == 0 && Dm.rank==0 )
-			printf("	Time=%i, Max variation=%f, Global variation=%f \n",count,GlobalMax,GlobalVar);
-
-		if (fabs(GlobalMax) < 1e-5){
-			if (Dm.rank==0) printf("	Exiting with max tolerance of 1e-5 \n");
-			count=timesteps;
-		}
-	}
-	PROFILE_STOP("Eikonal3D");
-	return GlobalVar;
-}
-
 
 inline int NLM3D( const Array<float> &Input, Array<float> &Mean, 
     const Array<float> &Distance, Array<float> &Output, const int d, const float h)
@@ -509,6 +345,7 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
     Array<float>& Dist, Array<float>& MultiScaleSmooth, Array<float>& NonLocalMean, 
     fillHalo<float>& fillFloat, const Domain& Dm, int nprocx )
 {
+    PROFILE_SCOPED(timer,"solve");
     // Compute the median filter on the sparse array
     Med3D( VOL, Mean );
     fillFloat.fill( Mean );
@@ -529,8 +366,9 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
 void refine( const Array<float>& Dist_coarse, 
     const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
     Array<float>& Dist, Array<float>& MultiScaleSmooth, Array<float>& NonLocalMean, 
-    fillHalo<float>& fillFloat, const Domain& Dm, int nprocx, bool solve_eikonal=true )
+    fillHalo<float>& fillFloat, const Domain& Dm, int nprocx, int level )
 {
+    PROFILE_SCOPED(timer,"refine");
     int ratio[3] = { int(Dist.size(0)/Dist_coarse.size(0)),
                      int(Dist.size(1)/Dist_coarse.size(1)),
                      int(Dist.size(2)/Dist_coarse.size(2)) };
@@ -573,13 +411,144 @@ void refine( const Array<float>& Dist_coarse,
     // Remove disconnected domains
     //removeDisconnected( ID, Dm );
     // Compute the distance using the segmented volume
-    if ( solve_eikonal ) {
-    	Eikonal3D( Dist, ID, Dm, ID.size(0)*nprocx );
+    if ( level > 0 ) {
+    	//Eikonal3D( Dist, ID, Dm, ID.size(0)*nprocx );
+        //CalcDist3D( Dist, ID, Dm );
+        CalcDistMultiLevel( Dist, ID, Dm );
 	    fillFloat.fill(Dist);
     }
 }
 
 
+// Remove regions that are likely noise by shrinking the volumes by dx,
+// removing all values that are more than dx+delta from the surface, and then
+// growing by dx+delta and intersecting with the original data
+void filter_final( Array<char>& ID, Array<float>& Dist,
+    fillHalo<float>& fillFloat, const Domain& Dm,
+    Array<float>& Mean, Array<float>& Dist1, Array<float>& Dist2 )
+{
+    PROFILE_SCOPED(timer,"filter_final");
+	int rank;
+	MPI_Comm_rank(Dm.Comm,&rank);
+    int Nx = Dm.Nx-2;
+    int Ny = Dm.Ny-2;
+    int Nz = Dm.Nz-2;
+    // Calculate the distance
+    CalcDistMultiLevel( Dist, ID, Dm );
+    fillFloat.fill(Dist);
+    // Compute the range to shrink the volume based on the L2 norm of the distance
+    Array<float> Dist0(Nx,Ny,Nz);
+    fillFloat.copy(Dist,Dist0);
+    float tmp = 0;
+    for (size_t i=0; i<Dist0.length(); i++)
+        tmp += Dist0(i)*Dist0(i);
+    tmp = sqrt( sumReduce(Dm.Comm,tmp) / sumReduce(Dm.Comm,(float)Dist0.length()) );
+    const float dx1 = 0.3*tmp;
+    const float dx2 = 1.05*dx1;
+    if (rank==0)
+        printf("   %0.1f %0.1f %0.1f\n",tmp,dx1,dx2);
+    // Update the IDs/Distance removing regions that are < dx of the range
+    Dist1 = Dist;
+    Dist2 = Dist;
+    Array<char> ID1 = ID;
+    Array<char> ID2 = ID;
+    for (size_t i=0; i<ID.length(); i++) {
+        ID1(i) = Dist(i)<-dx1 ? 1:0;
+        ID2(i) = Dist(i)> dx1 ? 1:0;
+    }
+    //Array<float> Dist1 = Dist;
+    //Array<float> Dist2 = Dist;
+    CalcDistMultiLevel( Dist1, ID1, Dm );
+    CalcDistMultiLevel( Dist2, ID2, Dm );
+    fillFloat.fill(Dist1);
+    fillFloat.fill(Dist2);
+    // Keep those regions that are within dx2 of the new volumes
+    Mean = Dist;
+    for (size_t i=0; i<ID.length(); i++) {
+        if ( Dist1(i)+dx2>0 && ID(i)<=0 ) {
+            Mean(i) = -1;
+        } else if ( Dist2(i)+dx2>0 && ID(i)>0 ) {
+            Mean(i) = 1;
+        } else {
+            Mean(i) = Dist(i)>0 ? 0.5:-0.5;
+        }
+    }
+    // Find regions of uncertainty that are entirely contained within another region
+    fillHalo<double> fillDouble(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
+    fillHalo<BlobIDType> fillInt(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
+    BlobIDArray GlobalBlobID;
+    DoubleArray SignDist(ID.size());
+    for (size_t i=0; i<ID.length(); i++)
+        SignDist(i) = fabs(Mean(i))==1 ? -1:1;
+    fillDouble.fill(SignDist);
+    DoubleArray Phase(ID.size());
+    Phase.fill(1);
+    ComputeGlobalBlobIDs( Nx, Ny, Nz, Dm.rank_info, Phase, SignDist, 0, 0, GlobalBlobID, Dm.Comm );
+    fillInt.fill(GlobalBlobID);
+    int N_blobs = maxReduce(Dm.Comm,GlobalBlobID.max()+1);
+    std::vector<float> mean(N_blobs,0);
+    std::vector<int> count(N_blobs,0);
+    for (int k=1; k<=Nz; k++) {
+        for (int j=1; j<=Ny; j++) {
+            for (int i=1; i<=Nx; i++) {
+                int id = GlobalBlobID(i,j,k);
+                if ( id >= 0 ) {
+                    if ( GlobalBlobID(i-1,j,k)<0 ) {
+                        mean[id] += Mean(i-1,j,k);
+                        count[id]++;
+                    }
+                    if ( GlobalBlobID(i+1,j,k)<0 ) {
+                        mean[id] += Mean(i+1,j,k);
+                        count[id]++;
+                    }
+                    if ( GlobalBlobID(i,j-1,k)<0 ) {
+                        mean[id] += Mean(i,j-1,k);
+                        count[id]++;
+                    }
+                    if ( GlobalBlobID(i,j+1,k)<0 ) {
+                        mean[id] += Mean(i,j+1,k);
+                        count[id]++;
+                    }
+                    if ( GlobalBlobID(i,j,k-1)<0 ) {
+                        mean[id] += Mean(i,j,k-1);
+                        count[id]++;
+                    }
+                    if ( GlobalBlobID(i,j,k+1)<0 ) {
+                        mean[id] += Mean(i,j,k+1);
+                        count[id]++;
+                    }
+                }
+            }
+        }
+    }
+    mean = sumReduce(Dm.Comm,mean);
+    count = sumReduce(Dm.Comm,count);
+    for (size_t i=0; i<mean.size(); i++)
+        mean[i] /= count[i];
+    /*if (rank==0) {
+        for (size_t i=0; i<mean.size(); i++)
+            printf("%i %0.4f\n",i,mean[i]);
+    }*/
+    for (size_t i=0; i<Mean.length(); i++) {
+        int id = GlobalBlobID(i);
+        if ( id >= 0 ) {
+            if ( fabs(mean[id]) > 0.95 ) {
+                // Isolated domain surrounded by one domain
+                GlobalBlobID(i) = -2;
+                Mean(i) = sign(mean[id]);
+            } else {
+                // Boarder volume, set to liquid
+                Mean(i) = 1;
+            }
+        }
+    }
+    // Perform the final segmentation and update the distance
+    fillFloat.fill(Mean);
+    segment( Mean, ID, 0.01 );
+    CalcDistMultiLevel( Dist, ID, Dm );
+    fillFloat.fill(Dist);
+}
+
 
 int main(int argc, char **argv)
 {
@@ -590,6 +559,7 @@ int main(int argc, char **argv)
 	MPI_Comm comm = MPI_COMM_WORLD;
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
+    Utilities::setErrorHandlers();
 	PROFILE_START("Main");
 
 	//std::vector<std::string> filenames;
@@ -806,14 +776,22 @@ int main(int argc, char **argv)
     if (rank==0)
         printf("Refine mesh\n");
     for (int i=int(Nx.size())-2; i>=0; i--) {
-        printf("   Refining to level %i\n",int(i));
+        if (rank==0)
+            printf("   Refining to level %i\n",int(i));
         refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i],
-            NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i>0 );
+            NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i );
     }
     PROFILE_STOP("Refine distance");
-removeDisconnected( ID[0], *Dm[0] );
 
+    // Perform a final filter
+    PROFILE_START("Filtering final domains");
+    if (rank==0)
+        printf("Filtering final domains\n");
+    Array<float> filter_Mean, filter_Dist1, filter_Dist2;
+    filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 );
+    PROFILE_STOP("Filtering final domains");
 
+//removeDisconnected( ID[0], *Dm[0] );
 
     // Write the results to visit
 	if (rank==0) printf("Writing output \n");
@@ -863,6 +841,29 @@ removeDisconnected( ID[0], *Dm[0] );
 	    meshData[i].vars.push_back(SmoothData);
         fillDouble[i]->copy( MultiScaleSmooth[i], SmoothData->data );
     }
+    #if 0
+        std::shared_ptr<IO::Variable> filter_Mean_var( new IO::Variable() );
+	    filter_Mean_var->name = "Mean";
+	    filter_Mean_var->type = IO::VolumeVariable;
+	    filter_Mean_var->dim = 1;
+	    filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]);
+	    meshData[0].vars.push_back(filter_Mean_var);
+        fillDouble[0]->copy( filter_Mean, filter_Mean_var->data );
+        std::shared_ptr<IO::Variable> filter_Dist1_var( new IO::Variable() );
+	    filter_Dist1_var->name = "Dist1";
+	    filter_Dist1_var->type = IO::VolumeVariable;
+	    filter_Dist1_var->dim = 1;
+	    filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]);
+	    meshData[0].vars.push_back(filter_Dist1_var);
+        fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data );
+        std::shared_ptr<IO::Variable> filter_Dist2_var( new IO::Variable() );
+	    filter_Dist2_var->name = "Dist2";
+	    filter_Dist2_var->type = IO::VolumeVariable;
+	    filter_Dist2_var->dim = 1;
+	    filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]);
+	    meshData[0].vars.push_back(filter_Dist2_var);
+        fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data );
+    #endif
 
 	// Write visulization data
 	IO::writeData( 0, meshData, 2, comm );