OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

committing lbpm_uCT_pp due to impending merge conflicts. Will resolve on next commit

Browse Source
This commit is contained in:
James E McClure 2016-06-27 13:58:59 -04:00
commit dc607718e6
13 changed files with 947 additions and 615 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

166
IO/netcdf.cpp
View File

@ -9,6 +9,7 @@
#include <netcdf.h>
#include <netcdf_par.h>
#define CHECK_NC_ERR( ERR ) \
@ -24,6 +25,45 @@
namespace netcdf {
// Convert nc_type to VariableType
static inline VariableType convertType( nc_type type )
{
VariableType type2 = UNKNOWN;
if ( type == NC_BYTE )
type2 = BYTE;
else if ( type == NC_CHAR )
type2 = STRING;
else if ( type == NC_SHORT )
type2 = SHORT;
else if ( type == NC_USHORT )
type2 = USHORT;
else if ( type == NC_INT )
type2 = INT;
else if ( type == NC_UINT )
type2 = UINT;
else if ( type == NC_INT64 )
type2 = INT64;
else if ( type == NC_UINT64 )
type2 = UINT64;
else if ( type == NC_FLOAT )
type2 = FLOAT;
else if ( type == NC_DOUBLE )
type2 = DOUBLE;
else
ERROR("Unknown type");
return type2;
}
// Get nc_type from the template
template<class T> inline nc_type getType();
template<> inline nc_type getType<char>() { return NC_CHAR; }
template<> inline nc_type getType<short>() { return NC_SHORT; }
template<> inline nc_type getType<int>() { return NC_INT; }
template<> inline nc_type getType<float>() { return NC_FLOAT; }
template<> inline nc_type getType<double>() { return NC_DOUBLE; }
// Function to reverse an array
template<class TYPE>
inline std::vector<TYPE> reverse( const std::vector<TYPE>& x )
@ -76,11 +116,36 @@ std::string VariableTypeName( VariableType type )
/****************************************************
* Open/close a file *
****************************************************/
int open( const std::string& filename )
int open( const std::string& filename, FileMode mode, MPI_Comm comm )
{
int fid = 0;
int err = nc_open( filename.c_str(), NC_NOWRITE, &fid );
CHECK_NC_ERR( err );
if ( comm == MPI_COMM_NULL ) {
if ( mode == READ ) {
int err = nc_open( filename.c_str(), NC_NOWRITE, &fid );
CHECK_NC_ERR( err );
} else if ( mode == WRITE ) {
int err = nc_open( filename.c_str(), NC_WRITE, &fid );
CHECK_NC_ERR( err );
} else if ( mode == CREATE ) {
int err = nc_create( filename.c_str(), NC_SHARE|NC_64BIT_OFFSET, &fid );
CHECK_NC_ERR( err );
} else {
ERROR("Unknown file mode");
}
} else {
if ( mode == READ ) {
int err = nc_open_par( filename.c_str(), NC_MPIPOSIX, comm, MPI_INFO_NULL, &fid );
CHECK_NC_ERR( err );
} else if ( mode == WRITE ) {
int err = nc_open_par( filename.c_str(), NC_WRITE|NC_MPIPOSIX, comm, MPI_INFO_NULL, &fid );
CHECK_NC_ERR( err );
} else if ( mode == CREATE ) {
int err = nc_create_par( filename.c_str(), NC_MPIPOSIX, comm, MPI_INFO_NULL, &fid );
CHECK_NC_ERR( err );
} else {
ERROR("Unknown file mode");
}
}
return fid;
}
void close( int fid )
@ -154,33 +219,6 @@ std::vector<std::string> getAttNames( int fid )
}
return att;
}
static inline VariableType convertType( nc_type type )
{
VariableType type2 = UNKNOWN;
if ( type == NC_BYTE )
type2 = BYTE;
else if ( type == NC_CHAR )
type2 = STRING;
else if ( type == NC_SHORT )
type2 = SHORT;
else if ( type == NC_USHORT )
type2 = USHORT;
else if ( type == NC_INT )
type2 = INT;
else if ( type == NC_UINT )
type2 = UINT;
else if ( type == NC_INT64 )
type2 = INT64;
else if ( type == NC_UINT64 )
type2 = UINT64;
else if ( type == NC_FLOAT )
type2 = FLOAT;
else if ( type == NC_DOUBLE )
type2 = DOUBLE;
else
ERROR("Unknown type");
return type2;
}
VariableType getVarType( int fid, const std::string& var )
{
int varid = -1;
@ -354,6 +392,74 @@ Array<std::string> getAtt<std::string>( int fid, const std::string& att )
}
/****************************************************
* Write an array to a file *
****************************************************/
std::vector<int> defDim( int fid, const std::vector<std::string>& names, const std::vector<int>& dims )
{
std::vector<int> dimid(names.size(),0);
for (size_t i=0; i<names.size(); i++) {
int err = nc_def_dim( fid, names[i].c_str(), dims[i], &dimid[i]);
CHECK_NC_ERR( err );
}
return dimid;
}
template<class TYPE>
int nc_put_vars_TYPE( int, int, const size_t*, const size_t*, const ptrdiff_t*, const TYPE* );
template<>
int nc_put_vars_TYPE<short>( 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>( 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<float>( 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<double>( 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<class TYPE>
void write( int fid, const std::string& var, const std::vector<int>& dimids,
const Array<TYPE>& data, const std::vector<size_t>& start,
const std::vector<size_t>& count, const std::vector<size_t>& stride )
{
// Define the variable
int varid = 0;
int err = nc_def_var( fid, var.c_str(), getType<TYPE>(), data.ndim(), dimids.data(), &varid );
CHECK_NC_ERR( err );
// exit define mode
err = nc_enddef( fid );
CHECK_NC_ERR( err );
// set the access method to use MPI/PnetCDF collective I/O
err = nc_var_par_access( fid, varid, NC_COLLECTIVE );
CHECK_NC_ERR( err );
// parallel write: each process writes its subarray to the file
auto x = data.reverseDim();
nc_put_vars_TYPE<TYPE>( fid, varid, start.data(), count.data(), (const ptrdiff_t*) stride.data(), x.data() );
}
template void write<short>( int fid, const std::string& var, const std::vector<int>& dimids,
const Array<short>& data, const std::vector<size_t>& start,
const std::vector<size_t>& count, const std::vector<size_t>& stride );
template void write<int>( int fid, const std::string& var, const std::vector<int>& dimids,
const Array<int>& data, const std::vector<size_t>& start,
const std::vector<size_t>& count, const std::vector<size_t>& stride );
template void write<float>( int fid, const std::string& var, const std::vector<int>& dimids,
const Array<float>& data, const std::vector<size_t>& start,
const std::vector<size_t>& count, const std::vector<size_t>& stride );
template void write<double>( int fid, const std::string& var, const std::vector<int>& dimids,
const Array<double>& data, const std::vector<size_t>& start,
const std::vector<size_t>& count, const std::vector<size_t>& stride );
}; // netcdf namespace
#else

28
IO/netcdf.h
View File

@ -5,6 +5,8 @@
#include <vector>
#include "common/Array.h"
#include "common/MPI_Helpers.h"
namespace netcdf {
@ -13,6 +15,10 @@ namespace netcdf {
//! Enum to hold variable type
enum VariableType { BYTE, SHORT, USHORT, INT, UINT, INT64, UINT64, FLOAT, DOUBLE, STRING, UNKNOWN };
//! Enum to hold variable type
enum FileMode { READ, WRITE, CREATE };
//! Convert the VariableType to a string
std::string VariableTypeName( VariableType type );
@ -22,8 +28,10 @@ std::string VariableTypeName( VariableType type );
* @detailed This function opens a netcdf file
* @return This function returns a handle to the file
* @param filename File to open
* @param mode Open the file for reading or writing
* @param comm MPI communicator to use (MPI_COMM_WORLD: don't use parallel netcdf)
*/
int open( const std::string& filename );
int open( const std::string& filename, FileMode mode, MPI_Comm comm=MPI_COMM_NULL );
/*!
@ -111,5 +119,23 @@ template<class TYPE>
Array<TYPE> getAtt( int fid, const std::string& att );
/*!
* @brief Write the dimensions
* @detailed This function writes the grid dimensions to netcdf.
* @param fid Handle to the open file
*/
std::vector<int> defDim( int fid, const std::vector<std::string>& names, const std::vector<int>& dims );
/*!
* @brief Write a variable
* @detailed This function writes a variable to netcdf.
* @param fid Handle to the open file
*/
template<class TYPE>
void write( int fid, const std::string& var, const std::vector<int>& dimids, const Array<TYPE>& data,
const std::vector<size_t>& start, const std::vector<size_t>& count, const std::vector<size_t>& stride );
}; // netcdf namespace
#endif

6
analysis/eikonal.h
View File

@ -21,7 +21,7 @@
* @param[in] timesteps Maximum number of timesteps to process
* @return Returns the global variation
*/
inline float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps);
float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm, const int timesteps);
/*!
@ -34,7 +34,7 @@ inline float Eikonal3D( Array<float> &Distance, const Array<char> &ID, const Dom
* @param[in] DM Domain information
* @return Returns the global variation
*/
inline void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm );
void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Domain &Dm );
/*!
@ -47,7 +47,7 @@ inline void CalcDist3D( Array<float> &Distance, const Array<char> &ID, const Dom
* @param[in] DM Domain information
* @return Returns the global variation
*/
inline void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm );
void CalcDistMultiLevel( Array<float> &Distance, const Array<char> &ID, const Domain &Dm );

1
analysis/eikonal.hpp
View File

@ -2,6 +2,7 @@
#define Eikonal_HPP_INC
#include "analysis/eikonal.h"
#include "common/imfilter.h"

149
analysis/filters.cpp Normal file
View File

@ -0,0 +1,149 @@
#include "analysis/filters.h"
#include "ProfilerApp.h"
void Med3D( const Array<float> &Input, Array<float> &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<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
// Just use a 3x3x3 window (hit recursively if needed)
imin = i-1;
jmin = j-1;
kmin = k-1;
imax = i+2;
jmax = j+2;
kmax = k+2;
// Populate the list with values in the window
int Number=0;
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
List[Number++] = Input(ii,jj,kk);
}
}
}
// Sort the first 5 entries and return the median
for (ii=0; ii<14; ii++){
for (jj=ii+1; jj<27; jj++){
if (List[jj] < List[ii]){
float tmp = List[ii];
List[ii] = List[jj];
List[jj] = tmp;
}
}
}
// Return the median
Output(i,j,k) = List[13];
}
}
}
PROFILE_STOP("Med3D");
}
int NLM3D( const Array<float> &Input, Array<float> &Mean,
const Array<float> &Distance, Array<float> &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<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
imin = std::max(0,i-d);
jmin = std::max(0,j-d);
kmin = std::max(0,k-d);
imax = std::min(Nx-1,i+d);
jmax = std::min(Ny-1,j+d);
kmax = std::min(Nz-1,k+d);
// Populate the list with values in the window
sum = 0; weight=0;
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
sum += Input(ii,jj,kk);
weight++;
}
}
}
Mean(i,j,k) = sum / weight;
}
}
}
// Compute the non-local means
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
if (fabs(Distance(i,j,k)) < THRESHOLD_DIST){
// compute the expensive non-local means
sum = 0; weight=0;
imin = std::max(0,i-d);
jmin = std::max(0,j-d);
kmin = std::max(0,k-d);
imax = std::min(Nx-1,i+d);
jmax = std::min(Ny-1,j+d);
kmax = std::min(Nz-1,k+d);
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
float tmp = Mean(i,j,k) - Mean(ii,jj,kk);
sum += exp(-tmp*tmp*h)*Input(ii,jj,kk);
weight += exp(-tmp*tmp*h);
}
}
}
returnCount++;
//Output(i,j,k) = Mean(i,j,k);
Output(i,j,k) = sum / weight;
}
else{
// Just return the mean
Output(i,j,k) = Mean(i,j,k);
}
}
}
}
// Return the number of sites where NLM was applied
PROFILE_STOP("NLM3D");
return returnCount;
}

28
analysis/filters.h Normal file
View File

@ -0,0 +1,28 @@
#ifndef Filters_H_INC
#define Filters_H_INC
#include "common/Array.h"
/*!
* @brief Filter image
* @details This routine performs a median filter
* @param[in] Input Input image
* @param[out] Output Output image
*/
void Med3D( const Array<float> &Input, Array<float> &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<float> &Input, Array<float> &Mean,
const Array<float> &Distance, Array<float> &Output, const int d, const float h);
#endif

395
analysis/uCT.cpp Normal file
View File

@ -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<class T>
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<float> &Coarse, Array<float> &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<Nz+1; k++){
int k0 = floor((k-0.5*hz)/hz);
int k1 = k0+1;
int k2 = k0+2;
float dz = ( (k+0.5) - (k0+0.5)*hz ) / hz;
ASSERT(k0>=-1&&k0<nz+1&&dz>=0&&dz<=1);
for (int j=-1; j<Ny+1; j++){
int j0 = floor((j-0.5*hy)/hy);
int j1 = j0+1;
int j2 = j0+2;
float dy = ( (j+0.5) - (j0+0.5)*hy ) / hy;
ASSERT(j0>=-1&&j0<ny+1&&dy>=0&&dy<=1);
for (int i=-1; i<Nx+1; i++){
int i0 = floor((i-0.5*hx)/hx);
int i1 = i0+1;
int i2 = i0+2;
float dx = ( (i+0.5) - (i0+0.5)*hx ) / hx;
ASSERT(i0>=-1&&i0<nx+1&&dx>=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<float>& VOL, const Array<float>& Dist, float sigma, Array<float>& MultiScaleSmooth, fillHalo<float>& fillFloat )
{
for (size_t i=0; i<VOL.length(); i++) {
// use exponential weight based on the distance
float dst = Dist(i);
float tmp = exp(-(dst*dst)/(sigma*sigma));
float value = dst>0 ? -1:1;
MultiScaleSmooth(i) = tmp*VOL(i) + (1-tmp)*value;
}
fillFloat.fill(MultiScaleSmooth);
}
// Segment the data
void segment( const Array<float>& data, Array<char>& ID, float tol )
{
ASSERT(data.size()==ID.size());
for (size_t i=0; i<data.length(); i++) {
if ( data(i) > tol )
ID(i) = 0;
else
ID(i) = 1;
}
}
// Remove disconnected phases
void removeDisconnected( Array<char>& 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<ID.length(); i++) {
SignDist(i) = (2*ID(i)-1);
Phase(i) = 1;
}
ComputeGlobalBlobIDs( ID.size(0)-2, ID.size(1)-2, ID.size(2)-2,
Dm.rank_info, Phase, SignDist, 0, 0, GlobalBlobID, Dm.Comm );
for (size_t i=0; i<ID.length(); i++) {
if ( GlobalBlobID(i) > 0 )
ID(i) = 0;
ID(i) = GlobalBlobID(i);
}
}
// Solve a level (without any coarse level information)
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 );
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<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, 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; i<ID.length(); i++) {
char id = Dist(i)>0 ? 1:0;
if ( id != ID(i) )
Dist(i) = 0;
}
fillFloat.fill( Dist );
std::function<float(int,const float*)> 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<imfilter::BC> BC(3,imfilter::BC::replicate);
std::vector<std::function<float(int,const float*)>> filter_set(3,filter_1D);
Dist = imfilter::imfilter_separable<float>( 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; i<ID.length(); i++) {
char id = Dist(i)>0 ? 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<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);
}
// Filter the original data
void filter_src( const Domain& Dm, Array<float>& src )
{
PROFILE_START("Filter source data");
int Nx = Dm.Nx-2;
int Ny = Dm.Ny-2;
int Nz = Dm.Nz-2;
fillHalo<float> fillFloat(Dm.Comm,Dm.rank_info,Nx,Ny,Nz,1,1,1,0,1);
// Perform a hot-spot filter on the data
std::vector<imfilter::BC> BC = { imfilter::BC::replicate, imfilter::BC::replicate, imfilter::BC::replicate };
std::function<float(const Array<float>&)> filter_3D = []( const Array<float>& 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<float(const Array<float>&)> filter_1D = []( const Array<float>& 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<float>( LOCVOL[0], {1,1,1}, filter_3D, BC );
std::vector<std::function<float(const Array<float>&)>> filter_set(3,filter_1D);
src = imfilter::imfilter_separable<float>( 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<Array<float>> H(3);
H[0] = imfilter::create_filter<float>( { Nh[0] }, "gaussian", &sigma[0] );
H[1] = imfilter::create_filter<float>( { Nh[1] }, "gaussian", &sigma[1] );
H[2] = imfilter::create_filter<float>( { Nh[2] }, "gaussian", &sigma[2] );
src = imfilter::imfilter_separable( src, H, BC );
fillFloat.fill( src );
PROFILE_STOP("Filter source data");
}

56
analysis/uCT.h Normal file
View File

@ -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<float> &Coarse, Array<float> &Fine );
// Smooth the data using the distance
void smooth( const Array<float>& VOL, const Array<float>& Dist, float sigma, Array<float>& MultiScaleSmooth, fillHalo<float>& fillFloat );
// Segment the data
void segment( const Array<float>& data, Array<char>& ID, float tol );
// Remove disconnected phases
void removeDisconnected( Array<char>& ID, const Domain& Dm );
// Solve a level (without any coarse level information)
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 );
// Refine a solution from a coarse grid to a fine grid
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, 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<char>& ID, Array<float>& Dist,
fillHalo<float>& fillFloat, const Domain& Dm,
Array<float>& Mean, Array<float>& Dist1, Array<float>& Dist2 );
// Filter the original data
void filter_src( const Domain& Dm, Array<float>& src );
#endif

37
common/Domain.cpp
View File

@ -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 *
********************************************************/

8
common/Domain.h
View File

@ -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,

2
tests/CMakeLists.txt
View File

@ -38,7 +38,7 @@ ADD_LBPM_TEST_1_2_4( testCommunication )
ADD_LBPM_TEST_1_2_4( testUtilities )
ADD_LBPM_TEST( TestWriter )
IF ( USE_NETCDF )
ADD_LBPM_PROVISIONAL_TEST( TestNetcdf )
ADD_LBPM_TEST_PARALLEL( TestNetcdf 8 )
ADD_LBPM_EXECUTABLE( lbpm_uCT_pp )
ENDIF()

97
tests/TestNetcdf.cpp
View File

@ -1,21 +1,76 @@
// Sequential blob analysis
// Reads parallel simulation data and performs connectivity analysis
// and averaging on a blob-by-blob basis
// James E. McClure 2014
// Test reading/writing netcdf files
#include "IO/netcdf.h"
#include "common/MPI_Helpers.h"
#include "common/Communication.h"
#include "common/UnitTest.h"
#include "ProfilerApp.h"
void load( const std::string filename )
void load( const std::string& );
void test_NETCDF( UnitTest& ut )
{
int fid = netcdf::open( filename );
const int rank = comm_rank( MPI_COMM_WORLD );
const int size = comm_size( MPI_COMM_WORLD );
int nprocx = 2;
int nprocy = 2;
int nprocz = 2;
RankInfoStruct info( rank, nprocx, nprocy, nprocz );
Array<float> data( 4, 5, 6 );
for (size_t i=0; i<data.length(); i++)
data(i) = i;
size_t x = info.ix*data.size(0);
size_t y = info.jy*data.size(1);
size_t z = info.kz*data.size(2);
const char* filename = "test.nc";
std::vector<int> dim = { 4*nprocx, 5*nprocy, 6*nprocz };
int fid = netcdf::open( filename, netcdf::CREATE, MPI_COMM_WORLD );
auto dims = netcdf::defDim( fid, {"X", "Y", "Z"}, dim );
netcdf::write( fid, "tmp", dims, data, {x,y,z}, data.size(), {1,1,1} );
netcdf::close( fid );
MPI_Barrier( MPI_COMM_WORLD );
// Read the contents of the file we created
fid = netcdf::open( filename, netcdf::READ );
Array<float> tmp = netcdf::getVar<float>( fid, "tmp" );
if ( (int)tmp.size(0)!=dim[0] || (int)tmp.size(1)!=dim[1] || (int)tmp.size(2)!=dim[2] ) {
ut.failure("Array sizes do not match");
return;
}
bool pass = true;
for (size_t i=0; i<data.size(0); i++) {
for (size_t j=0; j<data.size(1); j++) {
for (size_t k=0; k<data.size(2); k++) {
pass = pass && data(i,j,k) == tmp(i+x,j+y,k+z);
}
}
}
if ( pass ) {
ut.passes("write/read simple parallel file");
} else {
ut.failure("write/read simple parallel file");
}
}
inline void print( const std::string& name, const std::vector<size_t> size )
{
printf(" Reading variable %s (%i",name.c_str(),(int)size[0]);
for (size_t i=1; i<size.size(); i++)
printf(",%i",(int)size[i]);
printf(")\n");
}
void load( const std::string& filename )
{
printf("Reading %s\n",filename.c_str());
int fid = netcdf::open( filename, netcdf::READ );
std::vector<std::string> vars = netcdf::getVarNames( fid );
for (size_t i=0; i<vars.size(); i++) {
printf("Reading variable %s\n",vars[i].c_str());
netcdf::VariableType type = netcdf::getVarType( fid, vars[i] );
print( vars[i], netcdf::getVarDim(fid,vars[i]) );
if ( type == netcdf::STRING )
Array<std::string> tmp = netcdf::getVar<std::string>( fid, vars[i] );
else if ( type == netcdf::SHORT )
@ -26,7 +81,7 @@ void load( const std::string filename )
std::vector<std::string> attr = netcdf::getAttNames( fid );
for (size_t i=0; i<attr.size(); i++) {
printf("Reading attribute %s\n",attr[i].c_str());
printf(" Reading attribute %s\n",attr[i].c_str());
netcdf::VariableType type = netcdf::getAttType( fid, attr[i] );
if ( type == netcdf::STRING )
Array<std::string> tmp = netcdf::getAtt<std::string>( fid, attr[i] );
@ -39,19 +94,29 @@ void load( const std::string filename )
int main(int argc, char **argv)
{
// Initialize MPI
MPI_Init(&argc,&argv);
int rank = comm_rank(MPI_COMM_WORLD);
UnitTest ut;
PROFILE_START("Main");
std::vector<std::string> filenames;
if ( argc==0 ) {
printf("At least one filename must be specified\n");
return 1;
// Test reading existing netcdf files
if ( rank==0 ) {
for (int i=1; i<argc; i++)
load( argv[i] );
}
for (int i=1; i<argc; i++)
load( argv[i] );
// Test writing/reading netcdf file
test_NETCDF( ut );
// Print the results
ut.report();
int N_errors = ut.NumFailGlobal();
PROFILE_SAVE("TestNetcdf");
return 0;
// Close MPI
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
return N_errors;
}

589
tests/lbpm_uCT_pp.cpp
View File

@ -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<class T>
inline int sign( T x )
{
if ( x==0 )
return 0;
return x>0 ? 1:-1;
}
inline void Med3D( const Array<float> &Input, Array<float> &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<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
// Just use a 3x3x3 window (hit recursively if needed)
imin = i-1;
jmin = j-1;
kmin = k-1;
imax = i+2;
jmax = j+2;
kmax = k+2;
// Populate the list with values in the window
int Number=0;
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
List[Number++] = Input(ii,jj,kk);
}
}
}
// Sort the first 5 entries and return the median
for (ii=0; ii<14; ii++){
for (jj=ii+1; jj<27; jj++){
if (List[jj] < List[ii]){
float tmp = List[ii];
List[ii] = List[jj];
List[jj] = tmp;
}
}
}
// Return the median
Output(i,j,k) = List[13];
}
}
}
PROFILE_STOP("Med3D");
}
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 );
}
inline void InterpolateMesh( const Array<float> &Coarse, Array<float> &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<Nz+1; k++){
int k0 = floor((k-0.5*hz)/hz);
int k1 = k0+1;
int k2 = k0+2;
float dz = ( (k+0.5) - (k0+0.5)*hz ) / hz;
ASSERT(k0>=-1&&k0<nz+1&&dz>=0&&dz<=1);
for (int j=-1; j<Ny+1; j++){
int j0 = floor((j-0.5*hy)/hy);
int j1 = j0+1;
int j2 = j0+2;
float dy = ( (j+0.5) - (j0+0.5)*hy ) / hy;
ASSERT(j0>=-1&&j0<ny+1&&dy>=0&&dy<=1);
for (int i=-1; i<Nx+1; i++){
int i0 = floor((i-0.5*hx)/hx);
int i1 = i0+1;
int i2 = i0+2;
float dx = ( (i+0.5) - (i0+0.5)*hx ) / hx;
ASSERT(i0>=-1&&i0<nx+1&&dx>=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<float> &Input, Array<float> &Mean,
const Array<float> &Distance, Array<float> &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<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
imin = max(0,i-d);
jmin = max(0,j-d);
kmin = max(0,k-d);
imax = min(Nx-1,i+d);
jmax = min(Ny-1,j+d);
kmax = min(Nz-1,k+d);
// Populate the list with values in the window
sum = 0; weight=0;
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
sum += Input(ii,jj,kk);
weight++;
}
}
}
Mean(i,j,k) = sum / weight;
}
}
}
// Compute the non-local means
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
if (fabs(Distance(i,j,k)) < THRESHOLD_DIST){
// compute the expensive non-local means
sum = 0; weight=0;
imin = max(0,i-d);
jmin = max(0,j-d);
kmin = max(0,k-d);
imax = min(Nx-1,i+d);
jmax = min(Ny-1,j+d);
kmax = min(Nz-1,k+d);
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
for (ii=imin; ii<imax; ii++){
float tmp = Mean(i,j,k) - Mean(ii,jj,kk);
sum += exp(-tmp*tmp*h)*Input(ii,jj,kk);
weight += exp(-tmp*tmp*h);
}
}
}
returnCount++;
//Output(i,j,k) = Mean(i,j,k);
Output(i,j,k) = sum / weight;
}
else{
// Just return the mean
Output(i,j,k) = Mean(i,j,k);
}
}
}
}
// Return the number of sites where NLM was applied
PROFILE_STOP("NLM3D");
return returnCount;
}
// 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);
}
// Smooth the data using the distance
void smooth( const Array<float>& VOL, const Array<float>& Dist, float sigma, Array<float>& MultiScaleSmooth, fillHalo<float>& fillFloat )
{
for (size_t i=0; i<VOL.length(); i++) {
// use exponential weight based on the distance
float dst = Dist(i);
float tmp = exp(-(dst*dst)/(sigma*sigma));
float value = dst>0 ? -1:1;
MultiScaleSmooth(i) = tmp*VOL(i) + (1-tmp)*value;
}
fillFloat.fill(MultiScaleSmooth);
}
// Segment the data
void segment( const Array<float>& data, Array<char>& ID, float tol )
{
ASSERT(data.size()==ID.size());
for (size_t i=0; i<data.length(); i++) {
if ( data(i) > tol )
ID(i) = 0;
else
ID(i) = 1;
}
}
// Remove disconnected phases
void removeDisconnected( Array<char>& 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<ID.length(); i++) {
SignDist(i) = (2*ID(i)-1);
Phase(i) = 1;
}
ComputeGlobalBlobIDs( ID.size(0)-2, ID.size(1)-2, ID.size(2)-2,
Dm.rank_info, Phase, SignDist, 0, 0, GlobalBlobID, Dm.Comm );
for (size_t i=0; i<ID.length(); i++) {
if ( GlobalBlobID(i) > 0 )
ID(i) = 0;
ID(i) = GlobalBlobID(i);
}
}
// Solve a level (without any coarse level information)
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 );
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<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, 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; i<ID.length(); i++) {
char id = Dist(i)>0 ? 1:0;
if ( id != ID(i) )
Dist(i) = 0;
}
fillFloat.fill( Dist );
std::function<float(int,const float*)> 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<imfilter::BC> BC(3,imfilter::BC::replicate);
std::vector<std::function<float(int,const float*)>> filter_set(3,filter_1D);
Dist = imfilter::imfilter_separable<float>( 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; i<ID.length(); i++) {
char id = Dist(i)>0 ? 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<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)
{
@ -648,7 +125,7 @@ int main(int argc, char **argv)
// Read the subvolume of interest on each processor
PROFILE_START("ReadVolume");
int fid = netcdf::open(filename);
int fid = netcdf::open(filename,netcdf::READ);
std::string varname("VOLUME");
netcdf::VariableType type = netcdf::getVarType( fid, varname );
std::vector<size_t> dim = netcdf::getVarDim( fid, varname );
@ -675,50 +152,18 @@ 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<imfilter::BC> BC = { imfilter::BC::replicate, imfilter::BC::replicate, imfilter::BC::replicate };
std::function<float(const Array<float>&)> filter_3D = []( const Array<float>& 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<float(const Array<float>&)> filter_1D = []( const Array<float>& 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<float>( LOCVOL[0], {1,1,1}, filter_3D, BC );
std::vector<std::function<float(const Array<float>&)>> filter_set(3,filter_1D);
LOCVOL[0] = imfilter::imfilter_separable<float>( 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<Array<float>> H(3);
H[0] = imfilter::create_filter<float>( { Nh[0] }, "gaussian", &sigma[0] );
H[1] = imfilter::create_filter<float>( { Nh[1] }, "gaussian", &sigma[1] );
H[2] = imfilter::create_filter<float>( { 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;k<Nz[0]+2;k++) {
for (int j=0;j<Ny[0]+2;j++) {
<<<<<<< HEAD
for (int i=0;i<Nx[0]+2;i++) {
=======
for (int i=0;i<Nx[0]+1;i++) {
>>>>>>> 1e8ca14d857206c8aea62af646076c8f900fc84f
int iproc = Dm[0]->iproc;
int jproc = Dm[0]->jproc;
@ -824,7 +269,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<int> 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");