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

Update to uCT

Browse Source
This commit is contained in:
James E McClure 2018-11-17 19:53:09 -05:00
parent 8c8504b496
commit e0b00873e9
3 changed files with 63 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
analysis/uCT.cpp
View File

@ -148,7 +148,7 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
// Compute the median filter on the sparse array // Compute the median filter on the sparse array
Med3D( VOL, Mean ); Med3D( VOL, Mean );
fillFloat.fill( Mean ); fillFloat.fill( Mean );
segment( Mean, ID, 0.01 ); segment( Mean, ID, threshold );
// Compute the distance using the segmented volume // Compute the distance using the segmented volume
CalcDist( Dist, ID, Dm ); CalcDist( Dist, ID, Dm );
fillFloat.fill(Dist); fillFloat.fill(Dist);

17
sample_scripts/config_build_eos
View File

@ -7,11 +7,13 @@ module unload PrgEnv-intel
module load PrgEnv-gnu/6.0.4 module load PrgEnv-gnu/6.0.4
module unload gcc cmake module unload gcc cmake
module load gcc/6.3.0 module load gcc/6.3.0
module load cray-hdf5-parallel/1.10.0.1 module load cray-hdf5-parallel/1.10.2.0
module load cray-netcdf-hdf5parallel module load cray-netcdf-hdf5parallel
module load mercurial git module load mercurial git
module load cmake3/3.6.1 module load cmake3/3.6.1
export LD_LIBRARY_PATH=/usr/lib64:/lustre/atlas/proj-shared/geo106/eos/netcdf/lib:/lustre/atlas/proj-shared/geo106/eos/zlib/lib:/lustre/atlas/proj-shared/geo106/eos/hdf5/lib:$LD_LIBRARY_PATH
export MPICH_RDMA_ENABLED_CUDA=0 export MPICH_RDMA_ENABLED_CUDA=0
echo $GNU_VERSION echo $GNU_VERSION
@ -27,7 +29,8 @@ cmake \
-D CMAKE_C_COMPILER:PATH=cc \ -D CMAKE_C_COMPILER:PATH=cc \
-D CMAKE_CXX_COMPILER:PATH=CC \ -D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_CXX_COMPILER:PATH=CC \ -D CMAKE_CXX_COMPILER:PATH=CC \
-D CXX_STD=11 \ -D CMAKE_CXX_FLAGS="-fPIC" \
-D CMAKE_CXX_STANDARD=14 \
-D USE_TIMER=false \ -D USE_TIMER=false \
-D TIMER_DIRECTORY=${HOME}/timerutility/build/opt \ -D TIMER_DIRECTORY=${HOME}/timerutility/build/opt \
-D MPI_COMPILER:BOOL=TRUE \ -D MPI_COMPILER:BOOL=TRUE \
@ -35,10 +38,12 @@ cmake \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \ -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D USE_CUDA=0 \ -D USE_CUDA=0 \
-D CUDA_FLAGS="-arch sm_35" \ -D CUDA_FLAGS="-arch sm_35" \
-DBUILD_SHARED_LIBS=OFF \ -D USE_SILO=1 \
-D USE_NETCDF=1 \ -D SILO_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/silo \
-D NETCDF_DIRECTORY=$NETCDF_DIR \ -D HDF5_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/hdf5 \
-D HDF5_DIRECTORY=$HDF5_DIR \ -D HDF5_LIB=/lustre/atlas/proj-shared/geo106/eos/hdf5/lib/libhdf5.a \
-D USE_NETCDF=1 \
-D NETCDF_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/netcdf \
-D CMAKE_SKIP_RPATH=true \ -D CMAKE_SKIP_RPATH=true \
~/LBPM-WIA ~/LBPM-WIA

82
tests/lbpm_uCT_pp.cpp
View File

@ -73,8 +73,8 @@ int main(int argc, char **argv)
int nprocz = nproc[2]; int nprocz = nproc[2];
auto InputFile = uct_db->getScalar<std::string>( "InputFile" ); auto InputFile = uct_db->getScalar<std::string>( "InputFile" );
auto target = uct_db->getScalar<int>("target"); auto target = uct_db->getScalar<float>("target");
auto background = uct_db->getScalar<int>("background"); auto background = uct_db->getScalar<float>("background");
auto rough_cutoff = uct_db->getScalar<float>( "rough_cutoff" ); auto rough_cutoff = uct_db->getScalar<float>( "rough_cutoff" );
auto lamda = uct_db->getScalar<float>( "lamda" ); auto lamda = uct_db->getScalar<float>( "lamda" );
auto nlm_sigsq = uct_db->getScalar<float>( "nlm_sigsq" ); auto nlm_sigsq = uct_db->getScalar<float>( "nlm_sigsq" );
@ -91,6 +91,10 @@ int main(int argc, char **argv)
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs); printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz); printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
printf("target value = %f \n",target);
printf("background value = %f \n",background);
printf("cylinder center = %i, %i, %i \n",center[0],center[1],center[2]);
printf("cylinder radius = %f \n",CylRad);
} }
if ( nprocs < nprocx*nprocy*nprocz ){ if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors"); ERROR("Insufficient number of processors");
@ -188,22 +192,23 @@ int main(int argc, char **argv)
PROFILE_STOP("ReadVolume"); PROFILE_STOP("ReadVolume");
if (rank==0) printf("Read complete\n"); if (rank==0) printf("Read complete\n");
// Filter the original data // Filter the original data
filter_src( *Dm[0], LOCVOL[0] ); filter_src( *Dm[0], LOCVOL[0] );
// Set up the mask to be distance to cylinder (crop outside cylinder) // Set up the mask to be distance to cylinder (crop outside cylinder)
if (rank==0) printf("Cropping with cylinder: %i, %i, %i, radius=%f \n",Dm[0]->nprocx()*Nx[0],Dm[0]->nprocy()*Ny[0],Dm[0]->nprocz()*Nz[0],CylRad);
for (int k=0;k<Nz[0]+2;k++) { for (int k=0;k<Nz[0]+2;k++) {
for (int j=0;j<Ny[0]+2;j++) { for (int j=0;j<Ny[0]+2;j++) {
for (int i=0;i<Nx[0]+2;i++) { for (int i=0;i<Nx[0]+2;i++) {
int x=Dm[0]->iproc()*Nx[0]+i-1; float x= float(Dm[0]->iproc()*Nx[0]+i-1);
int y=Dm[0]->jproc()*Ny[0]+j-1; float y= float (Dm[0]->jproc()*Ny[0]+j-1);
int z=Dm[0]->kproc()*Nz[0]+k-1; float z= float(Dm[0]->kproc()*Nz[0]+k-1);
int cx = center[0] - offset[0]; float cx = float(center[0] - offset[0]);
int cy = center[1] - offset[1]; float cy = float(center[1] - offset[1]);
int cz = center[2] - offset[2]; float cz = float(center[2] - offset[2]);
// distance from the center line // distance from the center line
MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) ); MASK(i,j,k) = sqrt((z-cz)*(z-cz) + (y-cy)*(y-cy));
//if (sqrt(((z-cz)*(z-cz) + (y-cy)*(y-cy)) ) > CylRad) LOCVOL[0](i,j,k)=background;
} }
} }
} }
@ -211,16 +216,21 @@ int main(int argc, char **argv)
// Compute the means for the high/low regions // Compute the means for the high/low regions
// (should use automated mixture model to approximate histograms) // (should use automated mixture model to approximate histograms)
//float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) ); //float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) );
double THRESHOLD=0.5*(target+background); float THRESHOLD=0.5*(target+background);
double mean_plus=0; float mean_plus=0;
double mean_minus=0; float mean_minus=0;
float min_value = LOCVOL[0](0);
float max_value = LOCVOL[0](0);
int count_plus=0; int count_plus=0;
int count_minus=0; int count_minus=0;
for (int k=1;k<Nz[0]+1;k++) { for (int k=1;k<Nz[0]+1;k++) {
for (int j=1;j<Ny[0]+1;j++) { for (int j=1;j<Ny[0]+1;j++) {
for (int i=1;i<Nx[0]+1;i++) { for (int i=1;i<Nx[0]+1;i++) {
if (MASK(i,j,k) > 0.f ){
auto tmp = LOCVOL[0](i,j,k);
//LOCVOL[0](i,j,k) = MASK(i,j,k);
if (MASK(i,j,k) < CylRad ){
auto tmp = LOCVOL[0](i,j,k);
/* if ((tmp-background)*(tmp-target) > 0){ /* if ((tmp-background)*(tmp-target) > 0){
// direction to background / target is the same // direction to background / target is the same
if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background
@ -230,37 +240,47 @@ int main(int argc, char **argv)
if ( tmp > THRESHOLD ) { if ( tmp > THRESHOLD ) {
mean_plus += tmp; mean_plus += tmp;
count_plus++; count_plus++;
} else if ( tmp < -THRESHOLD ) { }
else {
mean_minus += tmp; mean_minus += tmp;
count_minus++; count_minus++;
} }
} if (tmp < min_value) min_value = tmp;
if (tmp > max_value) max_value = tmp;
}
} }
} }
} }
count_plus=sumReduce( Dm[0]->Comm, count_plus);
count_minus=sumReduce( Dm[0]->Comm, count_minus);
if (rank==0) printf("minimum value=%f, max value=%f \n",min_value,max_value);
if (rank==0) printf("plus=%i, minus=%i \n",count_plus,count_minus);
ASSERT( count_plus > 0 && count_minus > 0 ); ASSERT( count_plus > 0 && count_minus > 0 );
mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus ); MPI_Barrier(comm);
mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus ); mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / count_plus;
mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / count_minus;
MPI_Barrier(comm);
if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus); if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus);
//if (rank==0) printf("Scale the input data (size = %i) \n",LOCVOL[0].length());
MPI_Barrier(comm);
// Scale the source data to +-1.0
for (size_t i=0; i<LOCVOL[0].length(); i++) { for (size_t i=0; i<LOCVOL[0].length(); i++) {
if (MASK(i) < 0.f){ if ( MASK(i) > CylRad ){
LOCVOL[0](i) = 1.0; LOCVOL[0](i)=background;
} }
else if ( LOCVOL[0](i) >= 0 ) { if ( LOCVOL[0](i) >= THRESHOLD ) {
LOCVOL[0](i) /= mean_plus; auto tmp = LOCVOL[0](i)/ mean_plus;
LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f ); LOCVOL[0](i) = std::min( tmp, 1.0f );
} else {
LOCVOL[0](i) /= -mean_minus;
LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f );
} }
else {
auto tmp = -LOCVOL[0](i)/mean_minus;
LOCVOL[0](i) = std::max( tmp, -1.0f );
}
//LOCVOL[0](i) = MASK(i);
} }
// Fill the source data for the coarse meshes // Fill the source data for the coarse meshes
if (rank==0) printf("Coarsen the mesh for N_levels=%i \n",N_levels);
MPI_Barrier(comm);
PROFILE_START("CoarsenMesh"); PROFILE_START("CoarsenMesh");
for (int i=1; i<N_levels; i++) { for (int i=1; i<N_levels; i++) {
Array<float> filter(ratio[0],ratio[1],ratio[2]); Array<float> filter(ratio[0],ratio[1],ratio[2]);