diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3d327737..34ec9f13 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -40,6 +40,7 @@ ADD_LBPM_TEST( TestWriter ) IF ( USE_NETCDF ) ADD_LBPM_TEST_PARALLEL( TestNetcdf 8 ) ADD_LBPM_EXECUTABLE( lbpm_uCT_pp ) + ADD_LBPM_EXECUTABLE( lbpm_uCT_maskfilter ) ENDIF() # Sample test that will run with 1, 2, and 4 processors, failing with 4 or more procs diff --git a/tests/lbpm_uCT_maskfilter.cpp b/tests/lbpm_uCT_maskfilter.cpp new file mode 100644 index 00000000..90306bf9 --- /dev/null +++ b/tests/lbpm_uCT_maskfilter.cpp @@ -0,0 +1,421 @@ +// Sequential blob analysis +// Reads parallel simulation data and performs connectivity analysis +// and averaging on a blob-by-blob basis +// James E. McClure 2014 + +#include +#include +#include +#include +#include +#include +#include + +#include "common/Array.h" +#include "common/Domain.h" +#include "common/Communication.h" +#include "common/MPI_Helpers.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" + + +int main(int argc, char **argv) +{ + + // Initialize MPI + int rank, nprocs; + MPI_Init(&argc,&argv); + MPI_Comm comm = MPI_COMM_WORLD; + MPI_Comm_rank(comm,&rank); + MPI_Comm_size(comm,&nprocs); + Utilities::setErrorHandlers(); + PROFILE_START("Main"); + + //std::vector filenames; + if ( argc<2 ) { + if ( rank == 0 ) + printf("At least one filename must be specified\n"); + return 1; + } + std::string filename = std::string(argv[1]); + if ( rank == 0 ) + printf("Input data file: %s\n",filename.c_str()); + + //....................................................................... + // Reading the domain information file + //....................................................................... + int nprocx, nprocy, nprocz, nx, ny, nz, nspheres; + double Lx, Ly, Lz; + read_domain( rank, nprocs, comm, nprocx, nprocy, nprocz, nx, ny, nz, nspheres, Lx, Ly, Lz ); + int BC=0; + + + // Check that the number of processors >= the number of ranks + if ( rank==0 ) { + printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); + 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); + } + if ( nprocs < nprocx*nprocy*nprocz ){ + ERROR("Insufficient number of processors"); + } + + // Determine the maximum number of levels for the desired coarsen ratio + int ratio[3] = {4,4,4}; + std::vector Nx(1,nx), Ny(1,ny), Nz(1,nz); + while ( Nx.back()%ratio[0]==0 && Nx.back()>8 && + Ny.back()%ratio[1]==0 && Ny.back()>8 && + Nz.back()%ratio[2]==0 && Nz.back()>8 ) + { + Nx.push_back( Nx.back()/ratio[0] ); + Ny.push_back( Ny.back()/ratio[1] ); + Nz.push_back( Nz.back()/ratio[2] ); + } + int N_levels = Nx.size(); + + // Initialize the domain + std::vector> Dm(N_levels); + for (int i=0; iid[n] = 1; + Dm[i]->CommInit(comm); + } + + // array containing a distance mask + Array MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2); + + // Create the level data + std::vector> ID(N_levels); + std::vector> LOCVOL(N_levels); + std::vector> Dist(N_levels); + std::vector> MultiScaleSmooth(N_levels); + std::vector> Mean(N_levels); + std::vector> NonLocalMean(N_levels); + std::vector>> fillDouble(N_levels); + std::vector>> fillFloat(N_levels); + std::vector>> fillChar(N_levels); + for (int i=0; i(Nx[i]+2,Ny[i]+2,Nz[i]+2); + LOCVOL[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + Dist[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + MultiScaleSmooth[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + Mean[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + NonLocalMean[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + ID[i].fill(0); + LOCVOL[i].fill(0); + Dist[i].fill(0); + MultiScaleSmooth[i].fill(0); + Mean[i].fill(0); + NonLocalMean[i].fill(0); + fillDouble[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],1,1,1,0,1) ); + fillFloat[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],1,1,1,0,1) ); + fillChar[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],1,1,1,0,1) ); + } + + + // Read the subvolume of interest on each processor + PROFILE_START("ReadDistance"); + int distid = netcdf::open("Distance.nc",netcdf::READ); + std::string distvarname("Distance"); + netcdf::VariableType type = netcdf::getVarType( distid, "Distance" ); + std::vector dim = netcdf::getVarDim( distid, "Distance" ); + if ( rank == 0 ) { + printf("Reading %s (%s)\n","Distance.nc",netcdf::VariableTypeName(type).c_str()); + printf(" dims = %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2])); + } + { + RankInfoStruct info( rank, nprocx, nprocy, nprocz ); + size_t x = info.ix*nx; + size_t y = info.jy*ny; + size_t z = info.kz*nz; + + // Read the local data + Array MASKVALUES = netcdf::getVar( distid, distvarname, {x,y,z}, {nx,ny,nz}, {1,1,1}); + // Copy the data and fill the halos + MASK.fill(0); + fillFloat[0]->copy( MASKVALUES, MASK ); + fillFloat[0]->fill( MASK ); + + } + netcdf::close( distid ); + MPI_Barrier(comm); + PROFILE_STOP("ReadDistance"); + if (rank==0) printf("Finished reading distance =\n"); + + + // Read the subvolume of interest on each processor + PROFILE_START("ReadVolume"); + int fid = netcdf::open(filename,netcdf::READ); + std::string varname("VOLUME"); + type = netcdf::getVarType( fid, varname ); + dim = netcdf::getVarDim( fid, varname ); + if ( rank == 0 ) { + printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str()); + printf(" dims = %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2])); + } + { + RankInfoStruct info( rank, nprocx, nprocy, nprocz ); + int x = info.ix*nx; + int y = info.jy*ny; + int z = info.kz*nz; + // Read the local data + Array VOLUME = netcdf::getVar( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} ); + // Copy the data and fill the halos + LOCVOL[0].fill(0); + fillFloat[0]->copy( VOLUME, LOCVOL[0] ); + fillFloat[0]->fill( LOCVOL[0] ); + } + netcdf::close( fid ); + MPI_Barrier(comm); + PROFILE_STOP("ReadVolume"); + if (rank==0) printf("Read complete\n"); + + + // Filter the original 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; + int kproc = Dm[0]->kproc; + + int x=iproc*Nx[0]+i-1; + int y=jproc*Ny[0]+j-1; + int z=kproc*Nz[0]+k-1; + + int cx = 0.5*nprocx*Nx[0]; + int cy = 0.5*nprocy*Ny[0]; + int cz = 0.5*nprocz*Nz[0]; + + // distance from the center line + MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) ); + + } + } + } + */ + // Compute the means for the high/low regions + // (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 mean_plus=0; + float mean_minus=0; + int count_plus=0; + int count_minus=0; + for (int k=1;k 0.f ){ + auto tmp = LOCVOL[0](i,j,k); + if ( tmp > THRESHOLD ) { + mean_plus += tmp; + count_plus++; + } else if ( tmp < -THRESHOLD ) { + mean_minus += tmp; + count_minus++; + } + } + } + } + } + mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus ); + mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus ); + if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus); + + + // Scale the source data to +-1.0 + for (size_t i=0; i= 0 ) { + LOCVOL[0](i) /= mean_plus; + LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f ); + } else { + LOCVOL[0](i) /= -mean_minus; + LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f ); + } + } + + + // Fill the source data for the coarse meshes + PROFILE_START("CoarsenMesh"); + for (int i=1; i filter(ratio[0],ratio[1],ratio[2]); + filter.fill(1.0f/filter.length()); + Array tmp(Nx[i-1],Ny[i-1],Nz[i-1]); + fillFloat[i-1]->copy( LOCVOL[i-1], tmp ); + Array coarse = tmp.coarsen( filter ); + fillFloat[i]->copy( coarse, LOCVOL[i] ); + fillFloat[i]->fill( LOCVOL[i] ); + } + PROFILE_STOP("CoarsenMesh"); + + + // Initialize the coarse level + PROFILE_START("Solve coarse mesh"); + if (rank==0) + printf("Initialize coarse mesh\n"); + solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(), + NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx ); + PROFILE_STOP("Solve coarse mesh"); + + // Refine the solution + PROFILE_START("Refine distance"); + if (rank==0) + printf("Refine mesh\n"); + for (int i=int(Nx.size())-2; i>=0; i--) { + if (rank==0) + printf(" Refining to level %i\n",int(i)); + refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i], + NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i ); + } + PROFILE_STOP("Refine distance"); + + // Perform a final filter + PROFILE_START("Filtering final domains"); + if (rank==0) + printf("Filtering final domains\n"); + Array filter_Mean, filter_Dist1, filter_Dist2; + filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 ); + PROFILE_STOP("Filtering final domains"); + + //removeDisconnected( ID[0], *Dm[0] ); + + + // Write the 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"); + std::vector meshData(N_levels); + for (size_t i=0; i( new IO::DomainMesh(Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],Lx,Ly,Lz) ); + // Source data + std::shared_ptr OrigData( new IO::Variable() ); + OrigData->name = "Source Data"; + OrigData->type = IO::VolumeVariable; + OrigData->dim = 1; + OrigData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(OrigData); + fillDouble[i]->copy( LOCVOL[i], OrigData->data ); + // Non-Local Mean + std::shared_ptr NonLocMean( new IO::Variable() ); + NonLocMean->name = "Non-Local Mean"; + NonLocMean->type = IO::VolumeVariable; + NonLocMean->dim = 1; + NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(NonLocMean); + fillDouble[i]->copy( NonLocalMean[i], NonLocMean->data ); + // Segmented Data + std::shared_ptr SegData( new IO::Variable() ); + SegData->name = "Segmented Data"; + SegData->type = IO::VolumeVariable; + SegData->dim = 1; + SegData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(SegData); + fillDouble[i]->copy( ID[i], SegData->data ); + // Signed Distance + std::shared_ptr DistData( new IO::Variable() ); + DistData->name = "Signed Distance"; + DistData->type = IO::VolumeVariable; + DistData->dim = 1; + DistData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(DistData); + fillDouble[i]->copy( Dist[i], DistData->data ); + // Smoothed Data + std::shared_ptr SmoothData( new IO::Variable() ); + SmoothData->name = "Smoothed Data"; + SmoothData->type = IO::VolumeVariable; + SmoothData->dim = 1; + SmoothData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(SmoothData); + fillDouble[i]->copy( MultiScaleSmooth[i], SmoothData->data ); + } + #if 0 + std::shared_ptr filter_Mean_var( new IO::Variable() ); + filter_Mean_var->name = "Mean"; + filter_Mean_var->type = IO::VolumeVariable; + filter_Mean_var->dim = 1; + filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]); + meshData[0].vars.push_back(filter_Mean_var); + fillDouble[0]->copy( filter_Mean, filter_Mean_var->data ); + std::shared_ptr filter_Dist1_var( new IO::Variable() ); + filter_Dist1_var->name = "Dist1"; + filter_Dist1_var->type = IO::VolumeVariable; + filter_Dist1_var->dim = 1; + filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]); + meshData[0].vars.push_back(filter_Dist1_var); + fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data ); + std::shared_ptr filter_Dist2_var( new IO::Variable() ); + filter_Dist2_var->name = "Dist2"; + filter_Dist2_var->type = IO::VolumeVariable; + filter_Dist2_var->dim = 1; + filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]); + meshData[0].vars.push_back(filter_Dist2_var); + fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data ); + #endif + + // Write visulization data + IO::writeData( 0, meshData, 2, comm ); + if (rank==0) printf("Finished. \n"); + + /* for (k=0;kiproc; int jproc = Dm[0]->jproc; int kproc = Dm[0]->kproc;