From 1e8ca14d857206c8aea62af646076c8f900fc84f Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Mon, 27 Jun 2016 11:15:42 -0400 Subject: [PATCH] Refactoring lbpm_uCT_pp.cpp --- IO/netcdf.cpp | 36 ++- analysis/eikonal.h | 6 +- analysis/eikonal.hpp | 1 + analysis/filters.cpp | 149 +++++++++++ analysis/filters.h | 28 ++ analysis/uCT.cpp | 395 ++++++++++++++++++++++++++++ analysis/uCT.h | 56 ++++ common/Domain.cpp | 37 +++ common/Domain.h | 8 + tests/lbpm_uCT_pp.cpp | 587 ++---------------------------------------- 10 files changed, 732 insertions(+), 571 deletions(-) create mode 100644 analysis/filters.cpp create mode 100644 analysis/filters.h create mode 100644 analysis/uCT.cpp create mode 100644 analysis/uCT.h diff --git a/IO/netcdf.cpp b/IO/netcdf.cpp index 54ef2433..2539edf8 100644 --- a/IO/netcdf.cpp +++ b/IO/netcdf.cpp @@ -405,6 +405,28 @@ std::vector defDim( int fid, const std::vector& names, const s return dimid; } template +int nc_put_vars_TYPE( int, int, const size_t*, const size_t*, const ptrdiff_t*, const TYPE* ); +template<> +int nc_put_vars_TYPE( int fid, int varid, const size_t* start, const size_t* count, const ptrdiff_t* stride, const short* data ) +{ + return nc_put_vars_short( fid, varid, start, count, stride, data ); +} +template<> +int nc_put_vars_TYPE( int fid, int varid, const size_t* start, const size_t* count, const ptrdiff_t* stride, const int* data ) +{ + return nc_put_vars_int( fid, varid, start, count, stride, data ); +} +template<> +int nc_put_vars_TYPE( int fid, int varid, const size_t* start, const size_t* count, const ptrdiff_t* stride, const float* data ) +{ + return nc_put_vars_float( fid, varid, start, count, stride, data ); +} +template<> +int nc_put_vars_TYPE( int fid, int varid, const size_t* start, const size_t* count, const ptrdiff_t* stride, const double* data ) +{ + return nc_put_vars_double( fid, varid, start, count, stride, data ); +} +template void write( int fid, const std::string& var, const std::vector& dimids, const Array& data, const std::vector& start, const std::vector& count, const std::vector& stride ) @@ -421,12 +443,20 @@ void write( int fid, const std::string& var, const std::vector& dimids, CHECK_NC_ERR( err ); // parallel write: each process writes its subarray to the file auto x = data.reverseDim(); - nc_put_vars_float( fid, varid, start.data(), count.data(), (const ptrdiff_t*) stride.data(), x.data() ); + nc_put_vars_TYPE( fid, varid, start.data(), count.data(), (const ptrdiff_t*) stride.data(), x.data() ); } -template -void write( int fid, const std::string& var, const std::vector& dimids, +template void write( int fid, const std::string& var, const std::vector& dimids, + const Array& data, const std::vector& start, + const std::vector& count, const std::vector& stride ); +template void write( int fid, const std::string& var, const std::vector& dimids, + const Array& data, const std::vector& start, + const std::vector& count, const std::vector& stride ); +template void write( int fid, const std::string& var, const std::vector& dimids, const Array& data, const std::vector& start, const std::vector& count, const std::vector& stride ); +template void write( int fid, const std::string& var, const std::vector& dimids, + const Array& data, const std::vector& start, + const std::vector& count, const std::vector& stride ); diff --git a/analysis/eikonal.h b/analysis/eikonal.h index 5959ab20..8f4545d8 100644 --- a/analysis/eikonal.h +++ b/analysis/eikonal.h @@ -21,7 +21,7 @@ * @param[in] timesteps Maximum number of timesteps to process * @return Returns the global variation */ -inline float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps); +float Eikonal3D( Array &Distance, const Array &ID, const Domain &Dm, const int timesteps); /*! @@ -34,7 +34,7 @@ inline float Eikonal3D( Array &Distance, const Array &ID, const Dom * @param[in] DM Domain information * @return Returns the global variation */ -inline void CalcDist3D( Array &Distance, const Array &ID, const Domain &Dm ); +void CalcDist3D( Array &Distance, const Array &ID, const Domain &Dm ); /*! @@ -47,7 +47,7 @@ inline void CalcDist3D( Array &Distance, const Array &ID, const Dom * @param[in] DM Domain information * @return Returns the global variation */ -inline void CalcDistMultiLevel( Array &Distance, const Array &ID, const Domain &Dm ); +void CalcDistMultiLevel( Array &Distance, const Array &ID, const Domain &Dm ); diff --git a/analysis/eikonal.hpp b/analysis/eikonal.hpp index 25369ebd..7bccbfa9 100644 --- a/analysis/eikonal.hpp +++ b/analysis/eikonal.hpp @@ -2,6 +2,7 @@ #define Eikonal_HPP_INC #include "analysis/eikonal.h" +#include "common/imfilter.h" diff --git a/analysis/filters.cpp b/analysis/filters.cpp new file mode 100644 index 00000000..2dee1def --- /dev/null +++ b/analysis/filters.cpp @@ -0,0 +1,149 @@ +#include "analysis/filters.h" + +#include "ProfilerApp.h" + + +void Med3D( const Array &Input, Array &Output ) +{ + PROFILE_START("Med3D"); + // Perform a 3D Median filter on Input array with specified width + int i,j,k,ii,jj,kk; + int imin,jmin,kmin,imax,jmax,kmax; + + float *List; + List=new float[27]; + + int Nx = int(Input.size(0)); + int Ny = int(Input.size(1)); + int Nz = int(Input.size(2)); + + for (k=1; k &Input, Array &Mean, + const Array &Distance, Array &Output, const int d, const float h) +{ + PROFILE_START("NLM3D"); + // Implemenation of 3D non-local means filter + // d determines the width of the search volume + // h is a free parameter for non-local means (i.e. 1/sigma^2) + // Distance is the signed distance function + // If Distance(i,j,k) > THRESHOLD_DIST then don't compute NLM + + float THRESHOLD_DIST = float(d); + float weight, sum; + int i,j,k,ii,jj,kk; + int imin,jmin,kmin,imax,jmax,kmax; + int returnCount=0; + + int Nx = int(Input.size(0)); + int Ny = int(Input.size(1)); + int Nz = int(Input.size(2)); + + // Compute the local means + for (k=1; k &Input, Array &Output ); + + +/*! + * @brief Filter image + * @details This routine performs a non-linear local means filter + * @param[in] Input Input image + * @param[in] Mean Mean value + * @param[out] Output Output image + */ +int NLM3D( const Array &Input, Array &Mean, + const Array &Distance, Array &Output, const int d, const float h); + + +#endif diff --git a/analysis/uCT.cpp b/analysis/uCT.cpp new file mode 100644 index 00000000..c96511f4 --- /dev/null +++ b/analysis/uCT.cpp @@ -0,0 +1,395 @@ +#include "analysis/uCT.h" +#include "analysis/analysis.h" +#include "analysis/eikonal.h" +#include "analysis/filters.h" +#include "analysis/uCT.h" +#include "common/imfilter.h" + + +template +inline int sign( T x ) +{ + if ( x==0 ) + return 0; + return x>0 ? 1:-1; +} + + +inline float trilinear( float dx, float dy, float dz, float f1, float f2, + float f3, float f4, float f5, float f6, float f7, float f8 ) +{ + double f, dx2, dy2, dz2, h0, h1; + dx2 = 1.0 - dx; + dy2 = 1.0 - dy; + dz2 = 1.0 - dz; + h0 = ( dx * f2 + dx2 * f1 ) * dy2 + ( dx * f4 + dx2 * f3 ) * dy; + h1 = ( dx * f6 + dx2 * f5 ) * dy2 + ( dx * f8 + dx2 * f7 ) * dy; + f = h0 * dz2 + h1 * dz; + return ( f ); +} + + +void InterpolateMesh( const Array &Coarse, Array &Fine ) +{ + PROFILE_START("InterpolateMesh"); + + // Interpolate values from a Coarse mesh to a fine one + // This routine assumes cell-centered meshes with 1 ghost cell + + // Fine mesh + int Nx = int(Fine.size(0))-2; + int Ny = int(Fine.size(1))-2; + int Nz = int(Fine.size(2))-2; + + // Coarse mesh + int nx = int(Coarse.size(0))-2; + int ny = int(Coarse.size(1))-2; + int nz = int(Coarse.size(2))-2; + + // compute the stride + int hx = Nx/nx; + int hy = Ny/ny; + int hz = Nz/nz; + ASSERT(nx*hx==Nx); + ASSERT(ny*hy==Ny); + ASSERT(nz*hz==Nz); + + // value to map distance between meshes (since distance is in voxels) + // usually hx=hy=hz (or something very close) + // the mapping is not exact + // however, it's assumed the coarse solution will be refined + // a good guess is the goal here! + float mapvalue = sqrt(hx*hx+hy*hy+hz*hz); + + // Interpolate to the fine mesh + for (int k=-1; k=-1&&k0=0&&dz<=1); + for (int j=-1; j=-1&&j0=0&&dy<=1); + for (int i=-1; i=-1&&i0=0&&dx<=1); + float val = trilinear( dx, dy, dz, + Coarse(i1,j1,k1), Coarse(i2,j1,k1), Coarse(i1,j2,k1), Coarse(i2,j2,k1), + Coarse(i1,j1,k2), Coarse(i2,j1,k2), Coarse(i1,j2,k2), Coarse(i2,j2,k2) ); + Fine(i+1,j+1,k+1) = mapvalue*val; + } + } + } + PROFILE_STOP("InterpolateMesh"); +} + + +// Smooth the data using the distance +void smooth( const Array& VOL, const Array& Dist, float sigma, Array& MultiScaleSmooth, fillHalo& fillFloat ) +{ + for (size_t i=0; i0 ? -1:1; + MultiScaleSmooth(i) = tmp*VOL(i) + (1-tmp)*value; + } + fillFloat.fill(MultiScaleSmooth); +} + + +// Segment the data +void segment( const Array& data, Array& ID, float tol ) +{ + ASSERT(data.size()==ID.size()); + for (size_t i=0; i tol ) + ID(i) = 0; + else + ID(i) = 1; + } +} + + +// Remove disconnected phases +void removeDisconnected( Array& ID, const Domain& Dm ) +{ + // Run blob identification to remove disconnected volumes + BlobIDArray GlobalBlobID; + DoubleArray SignDist(ID.size()); + DoubleArray Phase(ID.size()); + for (size_t i=0; i 0 ) + ID(i) = 0; + ID(i) = GlobalBlobID(i); + } +} + + +// Solve a level (without any coarse level information) +void solve( const Array& VOL, Array& Mean, Array& ID, + Array& Dist, Array& MultiScaleSmooth, Array& NonLocalMean, + fillHalo& fillFloat, const Domain& Dm, int nprocx ) +{ + PROFILE_SCOPED(timer,"solve"); + // Compute the median filter on the sparse array + Med3D( VOL, Mean ); + fillFloat.fill( Mean ); + segment( Mean, ID, 0.01 ); + // Compute the distance using the segmented volume + Eikonal3D( Dist, ID, Dm, ID.size(0)*nprocx ); + fillFloat.fill(Dist); + smooth( VOL, Dist, 2.0, MultiScaleSmooth, fillFloat ); + // Compute non-local mean + int depth = 5; + float sigsq=0.1; + int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq); + fillFloat.fill(NonLocalMean); +} + + +// Refine a solution from a coarse grid to a fine grid +void refine( const Array& Dist_coarse, + const Array& VOL, Array& Mean, Array& ID, + Array& Dist, Array& MultiScaleSmooth, Array& NonLocalMean, + fillHalo& 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)) }; + // Interpolate the distance from the coarse to fine grid + InterpolateMesh( Dist_coarse, Dist ); + // Compute the median filter on the array and segment + Med3D( VOL, Mean ); + fillFloat.fill( Mean ); + segment( Mean, ID, 0.01 ); + // If the ID has the wrong distance, set the distance to 0 and run a simple filter to set neighbors to 0 + for (size_t i=0; i0 ? 1:0; + if ( id != ID(i) ) + Dist(i) = 0; + } + fillFloat.fill( Dist ); + std::function filter_1D = []( int N, const float* data ) + { + bool zero = data[0]==0 || data[2]==0; + return zero ? data[1]*1e-12 : data[1]; + }; + std::vector BC(3,imfilter::BC::replicate); + std::vector> filter_set(3,filter_1D); + Dist = imfilter::imfilter_separable( Dist, {1,1,1}, filter_set, BC ); + fillFloat.fill( Dist ); + // Smooth the volume data + float lambda = 2*sqrt(double(ratio[0]*ratio[0]+ratio[1]*ratio[1]+ratio[2]*ratio[2])); + smooth( VOL, Dist, lambda, MultiScaleSmooth, fillFloat ); + // Compute non-local mean + int depth = 3; + float sigsq = 0.1; + int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq); + fillFloat.fill(NonLocalMean); + segment( NonLocalMean, ID, 0.001 ); + for (size_t i=0; i0 ? 1:0; + if ( id!=ID(i) || fabs(Dist(i))<1 ) + Dist(i) = 2.0*ID(i)-1.0; + } + // Remove disconnected domains + //removeDisconnected( ID, Dm ); + // Compute the distance using the segmented volume + 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& ID, Array& Dist, + fillHalo& fillFloat, const Domain& Dm, + Array& Mean, Array& Dist1, Array& 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 Dist0(Nx,Ny,Nz); + fillFloat.copy(Dist,Dist0); + float tmp = 0; + for (size_t i=0; i ID1 = ID; + Array ID2 = ID; + for (size_t i=0; i dx1 ? 1:0; + } + //Array Dist1 = Dist; + //Array 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; i0 && 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 fillDouble(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); + fillHalo 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 mean(N_blobs,0); + std::vector 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= 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); +} + + + +// Filter the original data +void filter_src( const Domain& Dm, Array& src ) +{ + PROFILE_START("Filter source data"); + int Nx = Dm.Nx-2; + int Ny = Dm.Ny-2; + int Nz = Dm.Nz-2; + fillHalo fillFloat(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); + // Perform a hot-spot filter on the data + std::vector BC = { imfilter::BC::replicate, imfilter::BC::replicate, imfilter::BC::replicate }; + std::function&)> filter_3D = []( const Array& data ) + { + float min1 = std::min(data(0,1,1),data(2,1,1)); + float min2 = std::min(data(1,0,1),data(1,2,1)); + float min3 = std::min(data(1,1,0),data(1,1,2)); + float max1 = std::max(data(0,1,1),data(2,1,1)); + float max2 = std::max(data(1,0,1),data(1,2,1)); + float max3 = std::max(data(1,1,0),data(1,1,2)); + float min = std::min(min1,std::min(min2,min3)); + float max = std::max(max1,std::max(max2,max3)); + return std::max(std::min(data(1,1,1),max),min); + }; + std::function&)> filter_1D = []( const Array& data ) + { + float min = std::min(data(0),data(2)); + float max = std::max(data(0),data(2)); + return std::max(std::min(data(1),max),min); + }; + //LOCVOL[0] = imfilter::imfilter( LOCVOL[0], {1,1,1}, filter_3D, BC ); + std::vector&)>> filter_set(3,filter_1D); + src = imfilter::imfilter_separable( src, {1,1,1}, filter_set, BC ); + fillFloat.fill( src ); + // Perform a gaussian filter on the data + int Nh[3] = { 2, 2, 2 }; + float sigma[3] = { 1.0, 1.0, 1.0 }; + std::vector> H(3); + H[0] = imfilter::create_filter( { Nh[0] }, "gaussian", &sigma[0] ); + H[1] = imfilter::create_filter( { Nh[1] }, "gaussian", &sigma[1] ); + H[2] = imfilter::create_filter( { Nh[2] }, "gaussian", &sigma[2] ); + src = imfilter::imfilter_separable( src, H, BC ); + fillFloat.fill( src ); + PROFILE_STOP("Filter source data"); +} diff --git a/analysis/uCT.h b/analysis/uCT.h new file mode 100644 index 00000000..32954d05 --- /dev/null +++ b/analysis/uCT.h @@ -0,0 +1,56 @@ +#ifndef uCT_H_INC +#define uCT_H_INC + +#include "common/Array.h" +#include "common/Domain.h" +#include "common/Communication.h" + + + +/*! + * @brief Interpolate between meshes + * @details This routine interpolates from a coarse to a fine mesh + * @param[in] Coarse Coarse mesh solution + * @param[out] Fine Fine mesh solution + */ +void InterpolateMesh( const Array &Coarse, Array &Fine ); + + +// Smooth the data using the distance +void smooth( const Array& VOL, const Array& Dist, float sigma, Array& MultiScaleSmooth, fillHalo& fillFloat ); + + +// Segment the data +void segment( const Array& data, Array& ID, float tol ); + + +// Remove disconnected phases +void removeDisconnected( Array& ID, const Domain& Dm ); + + +// Solve a level (without any coarse level information) +void solve( const Array& VOL, Array& Mean, Array& ID, + Array& Dist, Array& MultiScaleSmooth, Array& NonLocalMean, + fillHalo& fillFloat, const Domain& Dm, int nprocx ); + + +// Refine a solution from a coarse grid to a fine grid +void refine( const Array& Dist_coarse, + const Array& VOL, Array& Mean, Array& ID, + Array& Dist, Array& MultiScaleSmooth, Array& NonLocalMean, + fillHalo& fillFloat, const Domain& Dm, int nprocx, int level ); + + +// 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& ID, Array& Dist, + fillHalo& fillFloat, const Domain& Dm, + Array& Mean, Array& Dist1, Array& Dist2 ); + + +// Filter the original data +void filter_src( const Domain& Dm, Array& src ); + + +#endif diff --git a/common/Domain.cpp b/common/Domain.cpp index 7c99e6d9..c112620b 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -20,6 +20,43 @@ static int MAX_BLOB_COUNT=50; using namespace std; + +// Reading the domain information file +void read_domain( int rank, int nprocs, MPI_Comm comm, + int& nprocx, int& nprocy, int& nprocz, int& nx, int& ny, int& nz, + int& nspheres, double& Lx, double& Ly, double& Lz ) +{ + if (rank==0){ + ifstream domain("Domain.in"); + domain >> nprocx; + domain >> nprocy; + domain >> nprocz; + domain >> nx; + domain >> ny; + domain >> nz; + domain >> nspheres; + domain >> Lx; + domain >> Ly; + domain >> Lz; + + } + MPI_Barrier(comm); + // Computational domain + //................................................. + MPI_Bcast(&nx,1,MPI_INT,0,comm); + MPI_Bcast(&ny,1,MPI_INT,0,comm); + MPI_Bcast(&nz,1,MPI_INT,0,comm); + MPI_Bcast(&nprocx,1,MPI_INT,0,comm); + MPI_Bcast(&nprocy,1,MPI_INT,0,comm); + MPI_Bcast(&nprocz,1,MPI_INT,0,comm); + MPI_Bcast(&nspheres,1,MPI_INT,0,comm); + MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); + MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); + MPI_Barrier(comm); +} + + /******************************************************** * Constructor/Destructor * ********************************************************/ diff --git a/common/Domain.h b/common/Domain.h index ac194570..02148e5f 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -18,6 +18,14 @@ using namespace std; + +//! Read the domain information file +void read_domain( int rank, int nprocs, MPI_Comm comm, + int& nprocx, int& nprocy, int& nprocz, int& nx, int& ny, int& nz, + int& nspheres, double& Lx, double& Ly, double& Lz ); + + +//! Class to hold domain info struct Domain{ // Default constructor Domain(int nx, int ny, int nz, int rnk, int npx, int npy, int npz, diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp index c4065d73..db973c48 100644 --- a/tests/lbpm_uCT_pp.cpp +++ b/tests/lbpm_uCT_pp.cpp @@ -15,541 +15,18 @@ #include "common/Domain.h" #include "common/Communication.h" #include "common/MPI_Helpers.h" -#include "common/imfilter.h" #include "IO/MeshDatabase.h" #include "IO/Mesh.h" #include "IO/Writer.h" #include "IO/netcdf.h" #include "analysis/analysis.h" #include "analysis/eikonal.h" +#include "analysis/filters.h" +#include "analysis/uCT.h" #include "ProfilerApp.h" -template -inline int sign( T x ) -{ - if ( x==0 ) - return 0; - return x>0 ? 1:-1; -} - - -inline void Med3D( const Array &Input, Array &Output ) -{ - PROFILE_START("Med3D"); - // Perform a 3D Median filter on Input array with specified width - int i,j,k,ii,jj,kk; - int imin,jmin,kmin,imax,jmax,kmax; - - float *List; - List=new float[27]; - - int Nx = int(Input.size(0)); - int Ny = int(Input.size(1)); - int Nz = int(Input.size(2)); - - for (k=1; k &Coarse, Array &Fine ) -{ - PROFILE_START("InterpolateMesh"); - - // Interpolate values from a Coarse mesh to a fine one - // This routine assumes cell-centered meshes with 1 ghost cell - - // Fine mesh - int Nx = int(Fine.size(0))-2; - int Ny = int(Fine.size(1))-2; - int Nz = int(Fine.size(2))-2; - - // Coarse mesh - int nx = int(Coarse.size(0))-2; - int ny = int(Coarse.size(1))-2; - int nz = int(Coarse.size(2))-2; - - // compute the stride - int hx = Nx/nx; - int hy = Ny/ny; - int hz = Nz/nz; - ASSERT(nx*hx==Nx); - ASSERT(ny*hy==Ny); - ASSERT(nz*hz==Nz); - - // value to map distance between meshes (since distance is in voxels) - // usually hx=hy=hz (or something very close) - // the mapping is not exact - // however, it's assumed the coarse solution will be refined - // a good guess is the goal here! - float mapvalue = sqrt(hx*hx+hy*hy+hz*hz); - - // Interpolate to the fine mesh - for (int k=-1; k=-1&&k0=0&&dz<=1); - for (int j=-1; j=-1&&j0=0&&dy<=1); - for (int i=-1; i=-1&&i0=0&&dx<=1); - float val = trilinear( dx, dy, dz, - Coarse(i1,j1,k1), Coarse(i2,j1,k1), Coarse(i1,j2,k1), Coarse(i2,j2,k1), - Coarse(i1,j1,k2), Coarse(i2,j1,k2), Coarse(i1,j2,k2), Coarse(i2,j2,k2) ); - Fine(i+1,j+1,k+1) = mapvalue*val; - } - } - } - PROFILE_STOP("InterpolateMesh"); -} - - -inline int NLM3D( const Array &Input, Array &Mean, - const Array &Distance, Array &Output, const int d, const float h) -{ - PROFILE_START("NLM3D"); - // Implemenation of 3D non-local means filter - // d determines the width of the search volume - // h is a free parameter for non-local means (i.e. 1/sigma^2) - // Distance is the signed distance function - // If Distance(i,j,k) > THRESHOLD_DIST then don't compute NLM - - float THRESHOLD_DIST = float(d); - float weight, sum; - int i,j,k,ii,jj,kk; - int imin,jmin,kmin,imax,jmax,kmax; - int returnCount=0; - - int Nx = int(Input.size(0)); - int Ny = int(Input.size(1)); - int Nz = int(Input.size(2)); - - // Compute the local means - for (k=1; k> nprocx; - domain >> nprocy; - domain >> nprocz; - domain >> nx; - domain >> ny; - domain >> nz; - domain >> nspheres; - domain >> Lx; - domain >> Ly; - domain >> Lz; - - } - MPI_Barrier(comm); - // Computational domain - //................................................. - MPI_Bcast(&nx,1,MPI_INT,0,comm); - MPI_Bcast(&ny,1,MPI_INT,0,comm); - MPI_Bcast(&nz,1,MPI_INT,0,comm); - MPI_Bcast(&nprocx,1,MPI_INT,0,comm); - MPI_Bcast(&nprocy,1,MPI_INT,0,comm); - MPI_Bcast(&nprocz,1,MPI_INT,0,comm); - MPI_Bcast(&nspheres,1,MPI_INT,0,comm); - MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm); - MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm); - MPI_Barrier(comm); -} - - - -// Smooth the data using the distance -void smooth( const Array& VOL, const Array& Dist, float sigma, Array& MultiScaleSmooth, fillHalo& fillFloat ) -{ - for (size_t i=0; i0 ? -1:1; - MultiScaleSmooth(i) = tmp*VOL(i) + (1-tmp)*value; - } - fillFloat.fill(MultiScaleSmooth); -} - - -// Segment the data -void segment( const Array& data, Array& ID, float tol ) -{ - ASSERT(data.size()==ID.size()); - for (size_t i=0; i tol ) - ID(i) = 0; - else - ID(i) = 1; - } -} - - -// Remove disconnected phases -void removeDisconnected( Array& ID, const Domain& Dm ) -{ - // Run blob identification to remove disconnected volumes - BlobIDArray GlobalBlobID; - DoubleArray SignDist(ID.size()); - DoubleArray Phase(ID.size()); - for (size_t i=0; i 0 ) - ID(i) = 0; - ID(i) = GlobalBlobID(i); - } -} - - -// Solve a level (without any coarse level information) -void solve( const Array& VOL, Array& Mean, Array& ID, - Array& Dist, Array& MultiScaleSmooth, Array& NonLocalMean, - fillHalo& fillFloat, const Domain& Dm, int nprocx ) -{ - PROFILE_SCOPED(timer,"solve"); - // Compute the median filter on the sparse array - Med3D( VOL, Mean ); - fillFloat.fill( Mean ); - segment( Mean, ID, 0.01 ); - // Compute the distance using the segmented volume - Eikonal3D( Dist, ID, Dm, ID.size(0)*nprocx ); - fillFloat.fill(Dist); - smooth( VOL, Dist, 2.0, MultiScaleSmooth, fillFloat ); - // Compute non-local mean - int depth = 5; - float sigsq=0.1; - int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq); - fillFloat.fill(NonLocalMean); -} - - -// Refine a solution from a coarse grid to a fine grid -void refine( const Array& Dist_coarse, - const Array& VOL, Array& Mean, Array& ID, - Array& Dist, Array& MultiScaleSmooth, Array& NonLocalMean, - fillHalo& 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)) }; - // Interpolate the distance from the coarse to fine grid - InterpolateMesh( Dist_coarse, Dist ); - // Compute the median filter on the array and segment - Med3D( VOL, Mean ); - fillFloat.fill( Mean ); - segment( Mean, ID, 0.01 ); - // If the ID has the wrong distance, set the distance to 0 and run a simple filter to set neighbors to 0 - for (size_t i=0; i0 ? 1:0; - if ( id != ID(i) ) - Dist(i) = 0; - } - fillFloat.fill( Dist ); - std::function filter_1D = []( int N, const float* data ) - { - bool zero = data[0]==0 || data[2]==0; - return zero ? data[1]*1e-12 : data[1]; - }; - std::vector BC(3,imfilter::BC::replicate); - std::vector> filter_set(3,filter_1D); - Dist = imfilter::imfilter_separable( Dist, {1,1,1}, filter_set, BC ); - fillFloat.fill( Dist ); - // Smooth the volume data - float lambda = 2*sqrt(double(ratio[0]*ratio[0]+ratio[1]*ratio[1]+ratio[2]*ratio[2])); - smooth( VOL, Dist, lambda, MultiScaleSmooth, fillFloat ); - // Compute non-local mean - int depth = 3; - float sigsq = 0.1; - int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq); - fillFloat.fill(NonLocalMean); - segment( NonLocalMean, ID, 0.001 ); - for (size_t i=0; i0 ? 1:0; - if ( id!=ID(i) || fabs(Dist(i))<1 ) - Dist(i) = 2.0*ID(i)-1.0; - } - // Remove disconnected domains - //removeDisconnected( ID, Dm ); - // Compute the distance using the segmented volume - 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& ID, Array& Dist, - fillHalo& fillFloat, const Domain& Dm, - Array& Mean, Array& Dist1, Array& 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 Dist0(Nx,Ny,Nz); - fillFloat.copy(Dist,Dist0); - float tmp = 0; - for (size_t i=0; i ID1 = ID; - Array ID2 = ID; - for (size_t i=0; i dx1 ? 1:0; - } - //Array Dist1 = Dist; - //Array 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; i0 && 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 fillDouble(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1); - fillHalo 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 mean(N_blobs,0); - std::vector 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= 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) { @@ -615,7 +92,7 @@ int main(int argc, char **argv) } // array containing a distance mask - Array MASK(Nx[i]+2,Ny[i]+2,Nz[i]+2); + Array MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2); // Create the level data std::vector> ID(N_levels); @@ -675,50 +152,14 @@ int main(int argc, char **argv) // Filter the original data - PROFILE_START("Filter source data"); - { - // Perform a hot-spot filter on the data - std::vector BC = { imfilter::BC::replicate, imfilter::BC::replicate, imfilter::BC::replicate }; - std::function&)> filter_3D = []( const Array& data ) - { - float min1 = std::min(data(0,1,1),data(2,1,1)); - float min2 = std::min(data(1,0,1),data(1,2,1)); - float min3 = std::min(data(1,1,0),data(1,1,2)); - float max1 = std::max(data(0,1,1),data(2,1,1)); - float max2 = std::max(data(1,0,1),data(1,2,1)); - float max3 = std::max(data(1,1,0),data(1,1,2)); - float min = std::min(min1,std::min(min2,min3)); - float max = std::max(max1,std::max(max2,max3)); - return std::max(std::min(data(1,1,1),max),min); - }; - std::function&)> filter_1D = []( const Array& data ) - { - float min = std::min(data(0),data(2)); - float max = std::max(data(0),data(2)); - return std::max(std::min(data(1),max),min); - }; - //LOCVOL[0] = imfilter::imfilter( LOCVOL[0], {1,1,1}, filter_3D, BC ); - std::vector&)>> filter_set(3,filter_1D); - LOCVOL[0] = imfilter::imfilter_separable( LOCVOL[0], {1,1,1}, filter_set, BC ); - fillFloat[0]->fill( LOCVOL[0] ); - // Perform a gaussian filter on the data - int Nh[3] = { 2, 2, 2 }; - float sigma[3] = { 1.0, 1.0, 1.0 }; - std::vector> H(3); - H[0] = imfilter::create_filter( { Nh[0] }, "gaussian", &sigma[0] ); - H[1] = imfilter::create_filter( { Nh[1] }, "gaussian", &sigma[1] ); - H[2] = imfilter::create_filter( { Nh[2] }, "gaussian", &sigma[2] ); - LOCVOL[0] = imfilter::imfilter_separable( LOCVOL[0], H, BC ); - fillFloat[0]->fill( LOCVOL[0] ); - } - PROFILE_STOP("Filter source data"); + filter_src( *Dm[0], LOCVOL[0] ); // Set up the mask to be distance to cylinder (crop outside cylinder) float CylRad=900; for (int k=0;kiproc; int jproc = Dm[0]->jproc; @@ -824,7 +265,23 @@ int main(int argc, char **argv) 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] ); + //removeDisconnected( ID[0], *Dm[0] ); + + + // Write the distance function to a netcdf file + const char* netcdf_filename = "Distance.nc"; + { + RankInfoStruct info( rank, nprocx, nprocy, nprocz ); + size_t x = info.ix*nx; + size_t y = info.jy*ny; + size_t z = info.kz*nz; + std::vector dim = { Nx[0]*nprocx, Ny[0]*nprocy, Nz[0]*nprocz }; + int fid = netcdf::open( netcdf_filename, netcdf::CREATE, MPI_COMM_WORLD ); + auto dims = netcdf::defDim( fid, {"X", "Y", "Z"}, dim ); + netcdf::write( fid, "Distance", dims, Dist[0], {x,y,z}, Dist[0].size(), {1,1,1} ); + netcdf::close( fid ); + } + // Write the results to visit if (rank==0) printf("Writing output \n");