For another project, please consult its home page to locate the current issue tracker.
View Issue Details [ Jump to Notes ] | [ Print ] | ||||||||
ID | Project | Category | View Status | Date Submitted | Last Update | ||||
0010178 | ParaView | Feature | public | 2010-01-22 13:24 | 2010-03-09 22:08 | ||||
Reporter | Ken Moreland | ||||||||
Assigned To | Ken Moreland | ||||||||
Priority | immediate | Severity | minor | Reproducibility | have not tried | ||||
Status | closed | Resolution | fixed | ||||||
Platform | OS | OS Version | |||||||
Product Version | |||||||||
Target Version | 3.8 | Fixed in Version | 3.8 | ||||||
Summary | 0010178: Efficient raw image reads on parallel file systems | ||||||||
Description | When reading raw image data from a file, the current image reader does uncoordinated reads of the data. When many processes are each trying to read a 3D subblock from a file, each process requests many small strips of data from the file. When a parallel file system gets saturated with lots of these small requests, the devices get swamped trying to satisfy all the independent requests. In contrast, reading from a parallel file system can be very efficient when reading large blocks of data. If the reads requests from the independent process are consolidated, they can be satisfied with one or a few large reads. MPI version 2 has a set of IO commands to allow you to do just this. I propose creating a subset of the image reader that takes advantage of the MPI IO calls to make parallel reads efficient. | ||||||||
Tags | No tags attached. | ||||||||
Project | |||||||||
Topic Name | |||||||||
Type | |||||||||
Attached Files | MPIImageReader.patch [^] (65,498 bytes) 2010-01-22 13:27 [Show Content] [Hide Content]diff --git a/Applications/ParaView/ParaViewReaders.xml b/Applications/ParaView/ParaViewReaders.xml index 7bbb906..424f242 100644 --- a/Applications/ParaView/ParaViewReaders.xml +++ b/Applications/ParaView/ParaViewReaders.xml @@ -40,6 +40,7 @@ <Proxy group="sources" name="popreader"/> <Proxy group="sources" name="AVSucdSeriesReader" /> <Proxy group="sources" name="MetaImageReader" /> + <Proxy group="sources" name="NrrdReader" /> <Proxy group="sources" name="Facet Proxy" /> <Proxy group="sources" name="PNGReader" /> <Proxy group="sources" name="TIFFReader" /> diff --git a/Servers/ServerManager/Resources/readers.xml b/Servers/ServerManager/Resources/readers.xml index 43e28a1..e572c45 100644 --- a/Servers/ServerManager/Resources/readers.xml +++ b/Servers/ServerManager/Resources/readers.xml @@ -3588,7 +3588,7 @@ </SourceProxy> <SourceProxy name="ImageReader" - class="vtkImageReader" + class="vtkMPIImageReader" label="Image Reader"> <Documentation short_help="Read raw regular rectilinear grid data from a file." @@ -3741,6 +3741,41 @@ <!-- End ImageReader --> </SourceProxy> + <SourceProxy name="NrrdReader" + class="vtkPNrrdReader" + label="Nrrd Reader"> + <Documentation short_help="Read raw image files with Nrrd meta data." + long_help="Read raw image files with Nrrd meta data."> + The Nrrd reader reads raw image data much like the Raw Image Reader + except that it will also read metadata information in the Nrrd format. + This means that the reader will automatically set information like file + dimensions. + + There are several limitations on what type of nrrd files we can + read. This reader only supports nrrd files in raw format. Other + encodings like ascii and hex will result in errors. When reading in + detached headers, this only supports reading one file that is + detached. + </Documentation> + + <StringVectorProperty name="FileName" + command="SetFileName" + animateable="0" + number_of_elements="1"> + <FileListDomain name="files" /> + <Documentation> + The name of the file to read (or the meta data file that will + point to the actual file). + </Documentation> + </StringVectorProperty> + + <Hints> + <ReaderFactory extensions="nrrd nhdr" + file_description="Nrrd Raw Image Files" /> + </Hints> + + </SourceProxy> <!-- NrrdReader --> + <SourceProxy name="XdmfReader2" class="vtkXdmfReader2" label="XDMF Reader2"> diff --git a/VTK/IO/vtkImageReader.cxx b/VTK/IO/vtkImageReader.cxx index 5375efd..05f848d 100644 --- a/VTK/IO/vtkImageReader.cxx +++ b/VTK/IO/vtkImageReader.cxx @@ -46,9 +46,8 @@ vtkImageReader::vtkImageReader() { this->DataVOI[idx*2] = this->DataVOI[idx*2 + 1] = 0; } - - // Left over from short reader - this->DataMask = 0xffff; + + this->DataMask = static_cast<vtkTypeUInt64>(~0UL); this->Transform = NULL; this->ScalarArrayName = NULL; @@ -211,7 +210,7 @@ void vtkImageReaderUpdate2(vtkImageReader *self, vtkImageData *data, int comp, pixelSkip; long filePos, correction = 0; unsigned long count = 0; - unsigned short DataMask; + vtkTypeUInt64 DataMask; unsigned long target; // Get the requested extents. @@ -327,7 +326,7 @@ void vtkImageReaderUpdate2(vtkImageReader *self, vtkImageData *data, for (idx0 = dataExtent[0]; idx0 <= dataExtent[1]; ++idx0) { // Copy pixel into the output. - if (DataMask == 0xffff) + if (DataMask == static_cast<vtkTypeUInt64>(~0UL)) { for (comp = 0; comp < pixelSkip; comp++) { @@ -336,10 +335,9 @@ void vtkImageReaderUpdate2(vtkImageReader *self, vtkImageData *data, } else { - // left over from short reader (what about other types. for (comp = 0; comp < pixelSkip; comp++) { - outPtr0[comp] = (OT)((short)(inPtr[comp]) & DataMask); + outPtr0[comp] = (OT)((vtkTypeUInt64)(inPtr[comp]) & DataMask); } } // move to next pixel diff --git a/VTK/IO/vtkImageReader.h b/VTK/IO/vtkImageReader.h index 68e5b28..95b1b85 100644 --- a/VTK/IO/vtkImageReader.h +++ b/VTK/IO/vtkImageReader.h @@ -45,17 +45,13 @@ public: vtkGetVector6Macro(DataVOI,int); // Description: - // Set/Get the Data mask. - vtkGetMacro(DataMask,unsigned short); - void SetDataMask(int val) - { - if (val == this->DataMask) - { - return; - } - this->DataMask = static_cast<unsigned short>(val); - this->Modified(); - } + // Set/Get the Data mask. The data mask is a simply integer whose bits are + // treated as a mask to the bits read from disk. That is, the data mask is + // bitwise-and'ed to the numbers read from disk. This ivar is stored as 64 + // bits, the largest mask you will need. The mask will be truncated to the + // data size required to be read (using the least significant bits). + vtkGetMacro(DataMask, vtkTypeUInt64); + vtkSetMacro(DataMask, vtkTypeUInt64); // Description: // Set/Get transformation matrix to transform the data from slice space @@ -82,7 +78,7 @@ protected: vtkImageReader(); ~vtkImageReader(); - unsigned short DataMask; // Mask each pixel with ... + vtkTypeUInt64 DataMask; vtkTransform *Transform; diff --git a/VTK/Parallel/CMakeLists.txt b/VTK/Parallel/CMakeLists.txt index ae3d156..6fdbb22 100644 --- a/VTK/Parallel/CMakeLists.txt +++ b/VTK/Parallel/CMakeLists.txt @@ -45,6 +45,7 @@ vtkExtractPiece.cxx vtkExtractPolyDataPiece.cxx vtkExtractUnstructuredGridPiece.cxx vtkExtractUserDefinedPiece.cxx +vtkMPIImageReader.cxx vtkMultiProcessController.cxx vtkMultiProcessStream.cxx vtkParallelFactory.cxx @@ -58,6 +59,7 @@ vtkPieceRequestFilter.cxx vtkPieceScalars.cxx vtkPKdTree.cxx vtkPLinearExtrusionFilter.cxx +vtkPNrrdReader.cxx vtkPOPReader.cxx ${_VTK_OPENFOAM_SOURCES} vtkPOutlineCornerFilter.cxx diff --git a/VTK/Parallel/Testing/Cxx/CMakeLists.txt b/VTK/Parallel/Testing/Cxx/CMakeLists.txt index 89a2240..2d65af5 100644 --- a/VTK/Parallel/Testing/Cxx/CMakeLists.txt +++ b/VTK/Parallel/Testing/Cxx/CMakeLists.txt @@ -349,6 +349,16 @@ IF(VTK_USE_DISPLAY AND VTK_USE_RENDERING) ${VTK_MPI_POSTFLAGS}) ENDIF(VTK_MPI_MAX_NUMPROCS GREATER 1) + IF(VTK_MPI_MAX_NUMPROCS GREATER 6) + ADD_TEST(ParallelIso-7proc-image + ${VTK_MPIRUN_EXE} ${VTK_MPI_PRENUMPROC_FLAGS} ${VTK_MPI_NUMPROC_FLAG} 7 ${VTK_MPI_PREFLAGS} + ${CXX_TEST_PATH}/ParallelIsoTest + -D ${VTK_DATA_ROOT} + -T ${VTK_BINARY_DIR}/Testing/Temporary + -V Baseline/Parallel/ParallelIso.7proc.png + ${VTK_MPI_POSTFLAGS}) + ENDIF(VTK_MPI_MAX_NUMPROCS GREATER 6) + ENDIF (VTK_MPIRUN_EXE) # # If we do not have the data, still run the tests that we can diff --git a/VTK/Parallel/Testing/Cxx/ExerciseMultiProcessController.cxx b/VTK/Parallel/Testing/Cxx/ExerciseMultiProcessController.cxx index bef35e7..8ea5db6 100644 --- a/VTK/Parallel/Testing/Cxx/ExerciseMultiProcessController.cxx +++ b/VTK/Parallel/Testing/Cxx/ExerciseMultiProcessController.cxx @@ -1052,6 +1052,22 @@ int ExerciseMultiProcessController(vtkMultiProcessController *controller) args.retval = 1; } + int color = (group1->GetLocalProcessId() >= 0) ? 1 : 2; + vtkMultiProcessController *subcontroller + = controller->PartitionController(color, 0); + subcontroller->SetSingleMethod(Run, &args); + subcontroller->SingleMethodExecute(); + subcontroller->Delete(); + + try + { + CheckSuccess(controller, !args.retval); + } + catch (ExerciseMultiProcessControllerError) + { + args.retval = 1; + } + return args.retval; } diff --git a/VTK/Parallel/Testing/Cxx/ParallelIso.cxx b/VTK/Parallel/Testing/Cxx/ParallelIso.cxx index 4ae5b4b..ce0478e 100644 --- a/VTK/Parallel/Testing/Cxx/ParallelIso.cxx +++ b/VTK/Parallel/Testing/Cxx/ParallelIso.cxx @@ -13,7 +13,7 @@ =========================================================================*/ // This example demonstrates the use of data parallelism in VTK. The -// pipeline ( vtkImageReader -> vtkContourFilter -> vtkElevationFilter ) +// pipeline ( vtkMPIImageReader -> vtkContourFilter -> vtkElevationFilter ) // is created in parallel and each process is assigned 1 piece to process. // All satellite processes send the result to the first process which // collects and renders them. @@ -27,10 +27,10 @@ #include "vtkContourFilter.h" #include "vtkDataSet.h" #include "vtkElevationFilter.h" -#include "vtkImageReader.h" #include "vtkMath.h" #include "vtkMPIController.h" #include "vtkParallelFactory.h" +#include "vtkPNrrdReader.h" #include "vtkPolyData.h" #include "vtkPolyDataMapper.h" #include "vtkTestUtilities.h" @@ -87,7 +87,7 @@ void SetIsoValueRMI(void *localArg, void* vtkNotUsed(remoteArg), // This will be called by all processes void MyMain( vtkMultiProcessController *controller, void *arg ) { - vtkImageReader *reader; + vtkPNrrdReader *reader; vtkContourFilter *iso; vtkElevationFilter *elev; int myid, numProcs; @@ -102,12 +102,13 @@ void MyMain( vtkMultiProcessController *controller, void *arg ) // Create the reader, the data file name might have // to be changed depending on where the data files are. char* fname = vtkTestUtilities::ExpandDataFileName(args->argc, args->argv, - "Data/headsq/quarter"); - reader = vtkImageReader::New(); - reader->SetDataByteOrderToLittleEndian(); - reader->SetDataExtent(0, 63, 0, 63, 1, 93); - reader->SetFilePrefix(fname); - reader->SetDataSpacing(3.2, 3.2, 1.5); + "Data/headsq/quarter.nhdr"); + reader = vtkPNrrdReader::New(); + reader->SetFileName(fname); + // reader->SetDataByteOrderToLittleEndian(); + // reader->SetDataExtent(0, 63, 0, 63, 1, 93); + // reader->SetFilePrefix(fname); + // reader->SetDataSpacing(3.2, 3.2, 1.5); delete[] fname; // Iso-surface. @@ -129,6 +130,9 @@ void MyMain( vtkMultiProcessController *controller, void *arg ) exec->SetUpdateNumberOfPieces(exec->GetOutputInformation(0), numProcs); exec->SetUpdatePiece(exec->GetOutputInformation(0), myid); + // Make sure all processes update at the same time. + elev->Update(); + if (myid != 0) { // If I am not the root process diff --git a/VTK/Parallel/vtkMPI.h b/VTK/Parallel/vtkMPI.h index 1e95eb5..f59a09b 100644 --- a/VTK/Parallel/vtkMPI.h +++ b/VTK/Parallel/vtkMPI.h @@ -52,6 +52,13 @@ public: MPI_Comm* Handle; }; +class VTK_PARALLEL_EXPORT vtkMPIOpaqueFileHandle +{ +public: + vtkMPIOpaqueFileHandle() : Handle(MPI_FILE_NULL) { } + MPI_File Handle; +}; + //----------------------------------------------------------------------------- class vtkMPICommunicatorOpaqueRequest { diff --git a/VTK/Parallel/vtkMPICommunicator.cxx b/VTK/Parallel/vtkMPICommunicator.cxx index 8a321b0..9fe3943 100644 --- a/VTK/Parallel/vtkMPICommunicator.cxx +++ b/VTK/Parallel/vtkMPICommunicator.cxx @@ -552,6 +552,53 @@ int vtkMPICommunicator::Initialize(vtkProcessGroup *group) return 1; } +//----------------------------------------------------------------------------- +int vtkMPICommunicator::SplitInitialize(vtkCommunicator *oldcomm, + int color, int key) +{ + if (this->Initialized) return 0; + + vtkMPICommunicator *mpiComm = vtkMPICommunicator::SafeDownCast(oldcomm); + if (!mpiComm) + { + vtkErrorMacro("Split communicator must be an MPI communicator."); + return 0; + } + + // If mpiComm has been initialized, it is guaranteed (unless the MPI calls + // return an error somewhere) to have valid Communicator. + if (!mpiComm->Initialized) + { + vtkWarningMacro("The communicator passed has not been initialized!"); + return 0; + } + + this->KeepHandleOff(); + + this->MPIComm->Handle = new MPI_Comm; + int err; + if ( (err = MPI_Comm_split(*(mpiComm->MPIComm->Handle), color, key, + this->MPIComm->Handle)) + != MPI_SUCCESS ) + { + delete this->MPIComm->Handle; + this->MPIComm->Handle = 0; + + char *msg = vtkMPIController::ErrorString(err); + vtkErrorMacro("MPI error occured: " << msg); + delete[] msg; + + return 0; + } + + this->InitializeNumberOfProcesses(); + this->Initialized = 1; + + this->Modified(); + + return 1; +} + //---------------------------------------------------------------------------- // Start the copying process void vtkMPICommunicator::InitializeCopy(vtkMPICommunicator* source) diff --git a/VTK/Parallel/vtkMPICommunicator.h b/VTK/Parallel/vtkMPICommunicator.h index d92f96d..4b30bfc 100644 --- a/VTK/Parallel/vtkMPICommunicator.h +++ b/VTK/Parallel/vtkMPICommunicator.h @@ -86,6 +86,11 @@ public: VTK_LEGACY(int Initialize(vtkMPICommunicator* mpiComm, vtkMPIGroup* group)); // Description: + // Used to initialize the communicator (i.e. create the underlying MPI_Comm) + // using MPI_Comm_split on the given communicator. + int SplitInitialize(vtkCommunicator *oldcomm, int color, int key); + + // Description: // Performs the actual communication. You will usually use the convenience // Send functions defined in the superclass. virtual int SendVoidArray(const void *data, vtkIdType length, int type, diff --git a/VTK/Parallel/vtkMPIController.cxx b/VTK/Parallel/vtkMPIController.cxx index 312243f..a8bc757 100644 --- a/VTK/Parallel/vtkMPIController.cxx +++ b/VTK/Parallel/vtkMPIController.cxx @@ -324,3 +324,19 @@ vtkMPIController *vtkMPIController::CreateSubController(vtkProcessGroup *group) controller->SetCommunicator(subcomm); return controller; } + +//----------------------------------------------------------------------------- +vtkMPIController *vtkMPIController::PartitionController(int localColor, + int localKey) +{ + VTK_CREATE(vtkMPICommunicator, subcomm); + + if (!subcomm->SplitInitialize(this->Communicator, localColor, localKey)) + { + return NULL; + } + + vtkMPIController *controller = vtkMPIController::New(); + controller->SetCommunicator(subcomm); + return controller; +} diff --git a/VTK/Parallel/vtkMPIController.h b/VTK/Parallel/vtkMPIController.h index ca3b635..02668b0 100644 --- a/VTK/Parallel/vtkMPIController.h +++ b/VTK/Parallel/vtkMPIController.h @@ -114,6 +114,8 @@ public: virtual vtkMPIController *CreateSubController(vtkProcessGroup *group); + virtual vtkMPIController *PartitionController(int localColor, int localKey); + //BTX // Description: diff --git a/VTK/Parallel/vtkMPIImageReader.cxx b/VTK/Parallel/vtkMPIImageReader.cxx new file mode 100644 index 0000000..0025c48 --- /dev/null +++ b/VTK/Parallel/vtkMPIImageReader.cxx @@ -0,0 +1,522 @@ +// -*- c++ -*- +/*========================================================================= + + Program: Visualization Toolkit + Module: $RCSfile$ + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/*---------------------------------------------------------------------------- + Copyright (c) Sandia Corporation + See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details. +----------------------------------------------------------------------------*/ + +#include "vtkMPIImageReader.h" + +#include "vtkMultiProcessController.h" +#include "vtkObjectFactory.h" +#include "vtkToolkits.h" + +#ifdef VTK_USE_MPI + +// Include the MPI headers and then determine if MPIIO is available. +#include "vtkMPI.h" + +#ifdef MPI_VERSION +# if (MPI_VERSION >= 2) +# define VTK_USE_MPI_IO 1 +# endif +#endif +#if !defined(VTK_USE_MPI_IO) && defined(ROMIO_VERSION) +# define VTK_USE_MPI_IO 1 +#endif +#if !defined(VTK_USE_MPI_IO) && defined(MPI_SEEK_SET) +# define VTK_USE_MPI_IO 1 +#endif + +#endif // VTK_USE_MPI + +// If VTK_USE_MPI_IO is set, that means we will read the data ourself using +// MPIIO. Otherwise, just delegate everything to the superclass. + +#ifdef VTK_USE_MPI_IO + +// We only need these includes if we are actually loading the data. +#include "vtkByteSwap.h" +#include "vtkDataArray.h" +#include "vtkImageData.h" +#include "vtkMPIController.h" +#include "vtkPointData.h" +#include "vtkTransform.h" +#include "vtkType.h" + +// This macro can be wrapped around MPI function calls to easily report errors. +// Reporting errors is more important with file I/O because, unlike network I/O, +// they usually don't terminate the program. +#define MPICall(funcall) \ + { \ + int __my_result = funcall; \ + if (__my_result != MPI_SUCCESS) \ + { \ + char errormsg[MPI_MAX_ERROR_STRING]; \ + int dummy; \ + MPI_Error_string(__my_result, errormsg, &dummy); \ + vtkErrorMacro(<< "Received error when calling" << endl \ + << #funcall << endl << endl \ + << errormsg); \ + } \ + } + +#endif // VTK_USE_MPI_IO + +//============================================================================= +vtkCxxRevisionMacro(vtkMPIImageReader, "$Revision$"); +vtkStandardNewMacro(vtkMPIImageReader); + +vtkCxxSetObjectMacro(vtkMPIImageReader, Controller, vtkMultiProcessController); +vtkCxxSetObjectMacro(vtkMPIImageReader, GroupedController, + vtkMultiProcessController); + +//----------------------------------------------------------------------------- +#ifdef VTK_USE_MPI_IO +template<class T> +inline void vtkMPIImageReaderMaskBits(T *data, vtkIdType length, + vtkTypeUInt64 _mask) +{ + T mask = (T)_mask; + + // If the mask is the identity, just return. + if ((_mask == (vtkTypeUInt64)~0UL) || (mask == (T)~0) || (_mask == 0)) return; + + for (vtkIdType i = 0; i < length; i++) + { + data[i] &= mask; + } +} + +// Override float and double because masking bits for them makes no sense. +VTK_TEMPLATE_SPECIALIZE +void vtkMPIImageReaderMaskBits(float *, vtkIdType, vtkTypeUInt64) +{ + return; +} +VTK_TEMPLATE_SPECIALIZE +void vtkMPIImageReaderMaskBits(double *, vtkIdType, vtkTypeUInt64) +{ + return; +} +#endif //VTK_USE_MPI_IO + +//----------------------------------------------------------------------------- +#ifdef VTK_USE_MPI_IO +namespace { + template<class T> + inline T MY_ABS(T x) { return (x < 0) ? -x : x; } + + template<class T> + inline T MY_MIN(T x, T y) { return (x < y) ? x : y; } +}; +#endif //VTK_USE_MPI_IO + +//============================================================================= +vtkMPIImageReader::vtkMPIImageReader() +{ + this->Controller = NULL; + this->SetController(vtkMultiProcessController::GetGlobalController()); + + this->GroupedController = NULL; +} + +vtkMPIImageReader::~vtkMPIImageReader() +{ + this->SetController(NULL); + this->SetGroupedController(NULL); +} + +void vtkMPIImageReader::PrintSelf(ostream &os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); + + os << indent << "Controller: " << this->Controller << endl; +} + +//----------------------------------------------------------------------------- +int vtkMPIImageReader::GetDataScalarTypeSize() +{ + switch (this->GetDataScalarType()) + { + vtkTemplateMacro(return sizeof(VTK_TT)); + default: + vtkErrorMacro("Unknown data type."); + return 0; + } +} + +//----------------------------------------------------------------------------- +#ifdef VTK_USE_MPI_IO +void vtkMPIImageReader::PartitionController(const int extent[6]) +{ + // Number of points in the z direction of the whole data. + int numZ = this->DataExtent[5] - this->DataExtent[4] + 1; + + if ((this->GetFileDimensionality() == 3) || (numZ == 1)) + { + // Everyone reads from the same single file. No need to partion controller. + this->SetGroupedController(this->Controller); + return; + } + + // The following algorithm will have overflow problems if there are more + // than 2^15 files. I doubt anyone will ever be crazy enough to set up a + // large 3D image with that many slice files, but just in case... + if (numZ >= 32768) + { + vtkErrorMacro("I do not support more than 32768 files."); + return; + } + + // Hash the Z extent. This is guaranteed to be unique for any pair of + // extents (within the constraint given above). + int extentHash = ( extent[4]+this->DataExtent[4] + + (extent[5]+this->DataExtent[4])*numZ ); + + vtkMultiProcessController *subController + = this->Controller->PartitionController(extentHash, 0); + this->SetGroupedController(subController); + subController->Delete(); +} +#else // VTK_USE_MPI_IO +void vtkMPIImageReader::PartitionController(const int *) +{ + vtkErrorMacro(<< "vtkMPIImageReader::PartitionController() called when MPIIO " + << "not available."); +} +#endif // VTK_USE_MPI_IO + +//----------------------------------------------------------------------------- +// Technically we should be returning a 64 bit number, but I doubt any header +// will be bigger than the value stored in an unsigned int. Thus, we just +// follow the convention of the superclass. +#ifdef VTK_USE_MPI_IO +unsigned long vtkMPIImageReader::GetHeaderSize(vtkMPIOpaqueFileHandle &file) +{ + if (this->ManualHeaderSize) + { + return this->HeaderSize; + } + else + { + this->ComputeDataIncrements(); + + MPI_Offset size; + MPICall(MPI_File_get_size(file.Handle, &size)); + return static_cast<unsigned long> + (size - this->DataIncrements[this->GetFileDimensionality()]); + } +} +#else // VTK_USE_MPI_IO +unsigned long vtkMPIImageReader::GetHeaderSize(vtkMPIOpaqueFileHandle &) +{ + vtkErrorMacro(<< "vtkMPIImageReader::GetHeaderSize() called when MPIIO " + << "not available."); + return 0; +} +#endif // VTK_USE_MPI_IO + +//----------------------------------------------------------------------------- +#ifdef VTK_USE_MPI_IO +void vtkMPIImageReader::SetupFileView(vtkMPIOpaqueFileHandle &file, + const int extent[6]) +{ + int arrayOfSizes[3]; + int arrayOfSubSizes[3]; + int arrayOfStarts[3]; + + for (int i = 0; i < this->GetFileDimensionality(); i++) + { + arrayOfSizes[i] = this->DataExtent[i*2+1] - this->DataExtent[i*2] + 1; + arrayOfSubSizes[i] = extent[i*2+1] - extent[i*2] + 1; + arrayOfStarts[i] = extent[i*2]; + } + // Adjust for base size of data type and tuple size. + int baseSize = this->GetDataScalarTypeSize() * this->NumberOfScalarComponents; + arrayOfSizes[0] *= baseSize; + arrayOfSubSizes[0] *= baseSize; + arrayOfStarts[0] *= baseSize; + + // Create a view in MPIIO. + MPI_Datatype view; + MPICall(MPI_Type_create_subarray(this->GetFileDimensionality(), + arrayOfSizes, arrayOfSubSizes, arrayOfStarts, + MPI_ORDER_FORTRAN, MPI_BYTE, &view)); + MPICall(MPI_Type_commit(&view)); + MPICall(MPI_File_set_view(file.Handle, this->GetHeaderSize(file), MPI_BYTE, + view, const_cast<char *>("native"), MPI_INFO_NULL)); + MPICall(MPI_Type_free(&view)); +} +#else // VTK_USE_MPI_IO +void vtkMPIImageReader::SetupFileView(vtkMPIOpaqueFileHandle &, const int[6]) +{ + vtkErrorMacro(<< "vtkMPIImageReader::SetupFileView() called when MPIIO " + << "not available."); +} +#endif // VTK_USE_MPI_IO + +//----------------------------------------------------------------------------- +#ifdef VTK_USE_MPI_IO +void vtkMPIImageReader::ReadSlice(int slice, const int extent[6], void *buffer) +{ + this->ComputeInternalFileName(slice); + + vtkMPICommunicator *mpiComm = vtkMPICommunicator::SafeDownCast( + this->GroupedController->GetCommunicator()); + + // Open the file for this slice. + vtkMPIOpaqueFileHandle file; + int result; + result = MPI_File_open(*mpiComm->GetMPIComm()->GetHandle(), + this->InternalFileName, MPI_MODE_RDONLY, + MPI_INFO_NULL, &file.Handle); + if (!(result == MPI_SUCCESS)) + { + vtkErrorMacro("Could not open file: " << this->InternalFileName); + return; + } + + // Set up the file view based on the extents. + this->SetupFileView(file, extent); + + // Figure out how many bytes to read. + vtkIdType length = this->GetDataScalarTypeSize(); + length *= this->NumberOfScalarComponents; + length *= extent[1]-extent[0]+1; + length *= extent[3]-extent[2]+1; + if (this->GetFileDimensionality() == 3) length *= extent[5]-extent[4]+1; + + if (length > VTK_INT_MAX) + { + vtkErrorMacro(<< "Cannot read more than " << VTK_INT_MAX + << " bytes at a time."); + } + + // Do the read. This is a coordinated parallel operation for efficiency. + MPICall(MPI_File_read_all(file.Handle, buffer, static_cast<int>(length), + MPI_BYTE, MPI_STATUS_IGNORE)); + + MPICall(MPI_File_close(&file.Handle)); +} +#else // VTK_USE_MPI_IO +void vtkMPIImageReader::ReadSlice(int, const int [6], void *) +{ + vtkErrorMacro(<< "vtkMPIImageReader::ReadSlice() called with MPIIO " + << "not available."); +} +#endif // VTK_USE_MPI_IO + +//----------------------------------------------------------------------------- +#ifdef VTK_USE_MPI_IO +// This method could be made a lot more efficient. +void vtkMPIImageReader::TransformData(vtkImageData *data) +{ + if (!this->Transform) return; + + vtkDataArray *fileData = data->GetPointData()->GetScalars(); + vtkDataArray *dataData = fileData->NewInstance(); + dataData->SetName(fileData->GetName()); + dataData->SetNumberOfComponents(fileData->GetNumberOfComponents()); + dataData->SetNumberOfTuples(fileData->GetNumberOfTuples()); + + int dataExtent[6]; + data->GetExtent(dataExtent); + + int fileExtent[6]; + this->ComputeInverseTransformedExtent(dataExtent, fileExtent); + + vtkIdType dataMinExtent[3]; + vtkIdType fileMinExtent[3]; + vtkIdType dataExtentSize[3]; + vtkIdType fileExtentSize[3]; + for (int i = 0; i < 3; i++) + { + dataMinExtent[i] = MY_MIN(dataExtent[2*i], dataExtent[2*i+1]); + fileMinExtent[i] = MY_MIN(fileExtent[2*i], fileExtent[2*i+1]); + dataExtentSize[i] = MY_ABS(dataExtent[2*i+1] - dataExtent[2*i]) + 1; + fileExtentSize[i] = MY_ABS(fileExtent[2*i+1] - fileExtent[2*i]) + 1; + } + + for (vtkIdType file_k = 0; file_k < fileExtentSize[2]; file_k++) + { + for (vtkIdType file_j = 0; file_j < fileExtentSize[1]; file_j++) + { + for (vtkIdType file_i = 0; file_i < fileExtentSize[0]; file_i++) + { + double fileXYZ[3]; + fileXYZ[0] = file_i + fileMinExtent[0]; + fileXYZ[1] = file_j + fileMinExtent[1]; + fileXYZ[2] = file_k + fileMinExtent[2]; + double dataXYZ[3]; + this->Transform->TransformPoint(fileXYZ, dataXYZ); + vtkIdType data_i = dataXYZ[0] - dataMinExtent[0]; + vtkIdType data_j = dataXYZ[1] - dataMinExtent[1]; + vtkIdType data_k = dataXYZ[2] - dataMinExtent[2]; + + vtkIdType fileTuple + = ((file_k*fileExtentSize[1] + file_j)*fileExtentSize[0]) + file_i; + vtkIdType dataTuple + = ((data_k*dataExtentSize[1] + data_j)*dataExtentSize[0]) + data_i; + + dataData->SetTuple(dataTuple, fileTuple, fileData); + } + } + } + + data->GetPointData()->SetScalars(dataData); + dataData->Delete(); +} +#else // VTK_USE_MPI_IO +void vtkMPIImageReader::TransformData(vtkImageData *) +{ + vtkErrorMacro(<< "vtkMPIImageReader::TransformData() called with MPIIO " + << "not available."); +} +#endif // VTK_USE_MPI_IO + +//----------------------------------------------------------------------------- +void vtkMPIImageReader::ExecuteData(vtkDataObject *output) +{ +#ifdef VTK_USE_MPI_IO + vtkMPIController *MPIController + = vtkMPIController::SafeDownCast(this->Controller); + if (!MPIController) + { + this->Superclass::ExecuteData(output); + return; + } + + vtkImageData *data = this->AllocateOutputData(output); + + if (!this->FileName && !this->FilePattern && !this->FileNames) + { + vtkErrorMacro("Either a valid FileName, FileNames, or FilePattern" + " must be specified."); + return; + } + + // VTK stores images in traditional "right handed" coordinates. That is, the + // origin is in the lower left corner. Many images, especially those with RGB + // colors, have the origin in the upper right corner. In this case, we have + // to flip the y axis. + vtkTransform *saveTransform = this->Transform; + if (!this->FileLowerLeft) + { + vtkTransform *newTransform = vtkTransform::New(); + if (this->Transform) + { + newTransform->Concatenate(this->Transform); + } + else + { + newTransform->Identity(); + } + newTransform->Scale(1.0, -1.0, 1.0); + this->Transform = newTransform; + } + + // Get information on data partion requested. + int inExtent[6]; + vtkIdType inIncrements[3]; + data->GetExtent(inExtent); + data->GetIncrements(inIncrements); + vtkDataArray *outputDataArray = data->GetPointData()->GetScalars(); + vtkIdType numValues = ( outputDataArray->GetNumberOfComponents() + + outputDataArray->GetNumberOfTuples() ); + + outputDataArray->SetName(this->ScalarArrayName); + + vtkDebugMacro("Reading extent: " + << inExtent[0] << ", " << inExtent[1] << ", " + << inExtent[2] << ", " << inExtent[3] << ", " + << inExtent[4] << ", " << inExtent[5]); + + // Respect the Transform. + int outExtent[6]; + vtkIdType outIncrements[3]; + this->ComputeInverseTransformedExtent(inExtent, outExtent); + + // The superclass' ComputeInverseTransformedIncrements does not give us + // increments we can use. It just reorders the inIncrements, (offsets in the + // target data structure). This does not give us valid offsets for the file. + // Instead, we just recompute them. + //this->ComputeInverseTransformedIncrements(inIncrements, outIncrements); + outIncrements[0] = inIncrements[0]; + outIncrements[1] = outIncrements[0]*(MY_ABS(outExtent[1]-outExtent[0])+1); + outIncrements[2] = outIncrements[1]*(MY_ABS(outExtent[3]-outExtent[2])+1); + + this->ComputeDataIncrements(); + + // Get information on data type. + int typeSize = this->GetDataScalarTypeSize(); + + // Group processes based on which files they read. + this->PartitionController(outExtent); + + // Get the pointer to the data buffer. Don't worry. We support all the + // data types. I am just casting it to a char (byte) so that I can do + // byte arithmetic on the data. + char *dataBuffer = reinterpret_cast<char *>(data->GetScalarPointer()); + + if (this->GetFileDimensionality() == 3) + { + // Everything is in one big file. Read it all in one shot. + this->ReadSlice(0, outExtent, dataBuffer); + } + else // this->GetFileDimensionality() == 2 + { + // Read everything slice-by-slice. + char *ptr = dataBuffer; + for (int slice = outExtent[4]; slice <= outExtent[5]; slice++) + { + this->UpdateProgress( (0.9*(slice-outExtent[4])) + / (outExtent[5]-outExtent[4]+1)); + this->ReadSlice(slice, outExtent, ptr); + ptr += typeSize*outIncrements[2]; + } + } + + this->UpdateProgress(0.9); + + // Swap bytes as necessary. + if (this->GetSwapBytes() && typeSize > 1) + { + vtkByteSwap::SwapVoidRange(dataBuffer, numValues, typeSize); + } + + // Mask bits as necessary. + switch (this->GetDataScalarType()) + { + vtkTemplateMacro(vtkMPIImageReaderMaskBits((VTK_TT *)dataBuffer, numValues, + this->DataMask)); + } + + // Perform permutation transformation of data if necessary. + this->TransformData(data); + + if (!this->FileLowerLeft) + { + this->Transform->Delete(); + this->Transform = saveTransform; + } + + // Done with this for now. + this->SetGroupedController(NULL); +#else // VTK_USE_MPI_IO + this->Superclass::ExecuteData(output); +#endif // VTK_USE_MPI_IO +} diff --git a/VTK/Parallel/vtkMPIImageReader.h b/VTK/Parallel/vtkMPIImageReader.h new file mode 100644 index 0000000..dc30faf --- /dev/null +++ b/VTK/Parallel/vtkMPIImageReader.h @@ -0,0 +1,120 @@ +// -*- c++ -*- +/*========================================================================= + + Program: Visualization Toolkit + Module: $RCSfile$ + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/*---------------------------------------------------------------------------- + Copyright (c) Sandia Corporation + See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details. +----------------------------------------------------------------------------*/ + +// .NAME vtkMPIImageReader Superclass of parallel binary image file readers. +// +// .SECTION Description +// +// vtkMPIImageReader provides the mechanism to read a brick of bytes (or shorts, +// or ints, or floats, or doubles, ...) from a file or series of files. You can +// use it to read raw image data from files. You may also be able to subclass +// this to read simple file formats. +// +// What distinguishes this class from vtkImageReader and vtkImageReader2 is that +// it performs synchronized parallel I/O using the MPIIO layer. This can make a +// huge difference in file read times, especially when reading in parallel from +// a parallel file system. +// +// Dispite the name of this class, vtkMPIImageReader will work even if MPI is +// not available. If MPI is not available or MPIIO is not available or the +// given Controller is not a vtkMPIController (or NULL), then this class will +// silently work exactly like its superclass. The point is that you can safely +// use this class in applications that may or may not be compiled with MPI (or +// may or may not actually be run with MPI). +// +// .SECTION See Also +// vtkMultiProcessController, vtkImageReader, vtkImageReader2 +// + +#ifndef __vtkMPIImageReader_h +#define __vtkMPIImageReader_h + +#include "vtkImageReader.h" + +class vtkMPIOpaqueFileHandle; +class vtkMultiProcessController; + +class VTK_PARALLEL_EXPORT vtkMPIImageReader : public vtkImageReader +{ +public: + vtkTypeRevisionMacro(vtkMPIImageReader, vtkImageReader); + static vtkMPIImageReader *New(); + virtual void PrintSelf(ostream &os, vtkIndent indent); + + // Description: + // Get/set the multi process controller to use for coordinated reads. By + // default, set to the global controller. + vtkGetObjectMacro(Controller, vtkMultiProcessController); + virtual void SetController(vtkMultiProcessController *); + +protected: + vtkMPIImageReader(); + ~vtkMPIImageReader(); + + vtkMultiProcessController *Controller; + + // Description: + // Returns the size, in bytes of the scalar data type (GetDataScalarType). + int GetDataScalarTypeSize(); + + // Description: + // Break up the controller based on the files each process reads. Each group + // comprises the processes that read the same files in the same order. + // this->GroupedController is set to the group for the current process. + virtual void PartitionController(const int extent[6]); + + // Description: + // Get the header size of the given open file. This should be used in liu of + // the GetHeaderSize methods of the superclass. + virtual unsigned long GetHeaderSize(vtkMPIOpaqueFileHandle &file); + + // Description: + // Set up a "view" on the open file that will allow you to read the 2D or 3D + // subarray from the file in one read. Once you call this method, the file + // will look as if it contains only the data the local process needs to read + // in. + virtual void SetupFileView(vtkMPIOpaqueFileHandle &file, const int extent[6]); + + // Description: + // Given a slice of the data, open the appropriate file, read the data into + // given buffer, and close the file. For three dimensional data, always + // use "slice" 0. Make sure the GroupedController is properly created before + // calling this using the PartitionController method. + virtual void ReadSlice(int slice, const int extent[6], void *buffer); + + // Description: + // Transform the data from the order read from a file to the order to place + // in the output data (as defined by the transform). + virtual void TransformData(vtkImageData *data); + + // Description: + // A group of processes that are reading the same file (as determined by + // PartitionController. + void SetGroupedController(vtkMultiProcessController *); + vtkMultiProcessController *GroupedController; + + virtual void ExecuteData(vtkDataObject *data); + +private: + vtkMPIImageReader(const vtkMPIImageReader &); // Not implemented + void operator=(const vtkMPIImageReader &); // Not implemented +}; + +#endif //__vtkMPIImageReader_h diff --git a/VTK/Parallel/vtkMultiProcessController.cxx b/VTK/Parallel/vtkMultiProcessController.cxx index 17f850e..f839ad8 100644 --- a/VTK/Parallel/vtkMultiProcessController.cxx +++ b/VTK/Parallel/vtkMultiProcessController.cxx @@ -29,6 +29,12 @@ #include "vtkMPIController.h" #endif +#include "vtkSmartPointer.h" +#define VTK_CREATE(type, name) \ + vtkSmartPointer<type> name = vtkSmartPointer<type>::New() + +#include <vtkstd/list> +#include <vtkstd/vector> #include <vtksys/hash_map.hxx> //---------------------------------------------------------------------------- @@ -318,6 +324,60 @@ vtkMultiProcessController *vtkMultiProcessController::CreateSubController( return subcontroller; } +//----------------------------------------------------------------------------- +vtkMultiProcessController *vtkMultiProcessController::PartitionController( + int localColor, + int localKey) +{ + vtkMultiProcessController *subController = NULL; + + int numProc = this->GetNumberOfProcesses(); + + vtkstd::vector<int> allColors(numProc); + this->AllGather(&localColor, &allColors[0], 1); + + vtkstd::vector<int> allKeys(numProc); + this->AllGather(&localKey, &allKeys[0], 1); + + vtkstd::vector<bool> inPartition; + inPartition.assign(numProc, false); + + for (int i = 0; i < numProc; i++) + { + if (inPartition[i]) continue; + int targetColor = allColors[i]; + vtkstd::list<int> partitionIds; // Make sorted list, then put in group. + for (int j = i; j < numProc; j++) + { + if (allColors[j] != targetColor) continue; + inPartition[j] = true; + vtkstd::list<int>::iterator iter = partitionIds.begin(); + while ((iter != partitionIds.end()) && (allKeys[*iter] <= allKeys[j])) + { + iter++; + } + partitionIds.insert(iter, j); + } + // Copy list into process group. + VTK_CREATE(vtkProcessGroup, group); + group->Initialize(this); + group->RemoveAllProcessIds(); + for (vtkstd::list<int>::iterator iter = partitionIds.begin(); + iter != partitionIds.end(); iter++) + { + group->AddProcessId(*iter); + } + // Use group to create controller. + vtkMultiProcessController *sc = this->CreateSubController(group); + if (sc) + { + subController = sc; + } + } + + return subController; +} + //---------------------------------------------------------------------------- int vtkMultiProcessController::RemoveFirstRMI(int tag) { diff --git a/VTK/Parallel/vtkMultiProcessController.h b/VTK/Parallel/vtkMultiProcessController.h index 8dff30b..c4ec81a 100644 --- a/VTK/Parallel/vtkMultiProcessController.h +++ b/VTK/Parallel/vtkMultiProcessController.h @@ -161,6 +161,18 @@ public: // returned on all process not in the group. virtual vtkMultiProcessController *CreateSubController( vtkProcessGroup *group); + + // Description: + // Partitions this controller based on a coloring. That is, each process + // passes in a color. All processes with the same color are grouped into the + // same partition. The processes are ordered by their self-assigned key. + // Lower keys have lower process ids. Ties are broken by the current process + // ids. (For example, if all the keys are 0, then the resulting processes + // will be ordered in the same way.) This method returns a new controller to + // each process that represents the local partition. This is basically the + // same operation as MPI_Comm_split. + virtual vtkMultiProcessController *PartitionController(int localColor, + int localKey); //------------------ RMIs -------------------- //BTX diff --git a/VTK/Parallel/vtkPNrrdReader.cxx b/VTK/Parallel/vtkPNrrdReader.cxx new file mode 100644 index 0000000..3949b82 --- /dev/null +++ b/VTK/Parallel/vtkPNrrdReader.cxx @@ -0,0 +1,625 @@ +// -*- c++ -*- +/*========================================================================= + + Program: Visualization Toolkit + Module: $RCSfile$ + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/*---------------------------------------------------------------------------- + Copyright (c) Sandia Corporation + See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details. +----------------------------------------------------------------------------*/ + +#include "vtkPNrrdReader.h" + +#include "vtkCharArray.h" +#include "vtkDummyController.h" +#include "vtkImageData.h" +#include "vtkObjectFactory.h" +#include "vtkStringArray.h" + +#include <vtkstd/string> +#include <vtkstd/vector> +#include <vtksys/SystemTools.hxx> +#include <istream> +#include <sstream> + +#include <math.h> +#include <string.h> + +#include "vtkSmartPointer.h" +#define VTK_CREATE(type, name) \ + vtkSmartPointer<type> name = vtkSmartPointer<type>::New() + +//============================================================================= +static vtkstd::string trim(vtkstd::string s) +{ + size_t start = 0; + while ((start < s.length()) && (s[start] == ' ')) start++; + size_t end = s.length(); + while ((end > start) && (s[end-1] == ' ')) end--; + return s.substr(start, end-start); +} + +//----------------------------------------------------------------------------- +static vtkstd::vector<vtkstd::string> split(vtkstd::string s) +{ + vtkstd::vector<vtkstd::string> result; + size_t startValue = 0; + while (true) + { + while ((startValue != s.npos) && (s[startValue] == ' ')) startValue++; + if (startValue == s.npos) return result; + size_t endValue = s.find(' ', startValue); + if (endValue != s.npos) + { + result.push_back(s.substr(startValue, endValue-startValue)); + } + else + { + result.push_back(s.substr(startValue)); + } + startValue = endValue; + } +} + +//----------------------------------------------------------------------------- +static void GetVector(vtkstd::string s, vtkstd::vector<int> &dest) +{ + vtkstd::vector<vtkstd::string> strlist = split(s); + for (size_t i = 0; i < dest.size(); i++) + { + if (i < strlist.size()) + { + dest[i] = atoi(strlist[i].c_str()); + } + else + { + dest[i] = 0; + } + } +} + +//----------------------------------------------------------------------------- +static void GetVector(vtkstd::string s, vtkstd::vector<double> &dest) +{ + vtkstd::vector<vtkstd::string> strlist = split(s); + for (size_t i = 0; i < dest.size(); i++) + { + if (i < strlist.size()) + { + dest[i] = atof(strlist[i].c_str()); + } + else + { + dest[i] = 0.0; + } + } +} + +//----------------------------------------------------------------------------- +static vtkstd::vector<double> ParseVector(vtkstd::string s) +{ + vtkstd::vector<double> result; + + s = trim(s); + if ((s[0] != '(') || (s[s.length()-1] != ')')) return result; + s = s.substr(1, s.length()-2); + while (true) + { + size_t i = s.find(','); + vtkstd::string value = s.substr(0, i); + result.push_back(atof(value.c_str())); + if (i == s.npos) break; + s = s.substr(i+1); + } + + return result; +} + +//----------------------------------------------------------------------------- +static int NrrdType2VTKType(vtkstd::string nrrdType) +{ + nrrdType = trim(nrrdType); + if ( (nrrdType == "signed char") + || (nrrdType == "int8") || (nrrdType == "int8_t") ) + { + return VTK_CHAR; + } + else if ( (nrrdType == "uchar") || (nrrdType == "unsigned char") + || (nrrdType == "uint8") || (nrrdType == "uint8_t") ) + { + return VTK_UNSIGNED_CHAR; + } + else if ( (nrrdType == "short") || (nrrdType == "short int") + || (nrrdType == "signed short") || (nrrdType == "signed short int") + || (nrrdType == "int16") || (nrrdType == "int16_t") ) + { + return VTK_SHORT; + } + else if ( (nrrdType == "ushort") || (nrrdType == "unsigned short") + || (nrrdType == "unsigned short int") + || (nrrdType == "uint16") || (nrrdType == "uint16_t") ) + { + return VTK_UNSIGNED_SHORT; + } + else if ( (nrrdType == "int") || (nrrdType == "signed int") + || (nrrdType == "int32") || (nrrdType == "int32_t") ) + { + return VTK_INT; + } + else if ( (nrrdType == "uint") || (nrrdType == "unsigned int") + || (nrrdType == "uint32") || (nrrdType == "uint32_t") ) + { + return VTK_UNSIGNED_INT; + } + else if ( (nrrdType == "longlong") || (nrrdType == "long long") + || (nrrdType == "long long int") || (nrrdType == "signed long long") + || (nrrdType == "signed long long int") || (nrrdType == "int64") + || (nrrdType == "int64_t") ) + { + return VTK_TYPE_INT64; + } + else if ( (nrrdType == "ulonglong") || (nrrdType == "unsigned long long") + || (nrrdType == "unsigned long long int") + || (nrrdType == "uint64") || (nrrdType == "uint64_t") ) + { + return VTK_TYPE_UINT64; + } + else if (nrrdType == "float") + { + return VTK_FLOAT; + } + else if (nrrdType == "double") + { + return VTK_DOUBLE; + } + else if (nrrdType == "block") + { + vtkGenericWarningMacro(<< "Reading blocks not supported."); + return VTK_VOID; + } + else + { + vtkGenericWarningMacro(<< "Unknown type: " << nrrdType); + return VTK_VOID; + } +} + +//============================================================================= +vtkCxxRevisionMacro(vtkPNrrdReader, "$Revision: 14830 $"); +vtkStandardNewMacro(vtkPNrrdReader); + +//----------------------------------------------------------------------------- +vtkPNrrdReader::vtkPNrrdReader() +{ + this->DataFiles = vtkStringArray::New(); +} + +vtkPNrrdReader::~vtkPNrrdReader() +{ + this->DataFiles->Delete(); + this->DataFiles = NULL; +} + +void vtkPNrrdReader::PrintSelf(ostream &os, vtkIndent indent) +{ + this->Superclass::PrintSelf(os, indent); +} + +//----------------------------------------------------------------------------- +int vtkPNrrdReader::CanReadFile(const char *filename) +{ + vtkstd::ifstream file(filename, + vtkstd::ios_base::in | vtkstd::ios_base::binary); + vtkstd::string firstLine; + vtkstd::getline(file, firstLine); + if (firstLine.substr(0, 4) == "NRRD") + { + return 2; + } + else + { + return 0; + } +} + +//----------------------------------------------------------------------------- +int vtkPNrrdReader::RequestInformation(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) +{ + if (!this->ReadHeader()) return 0; + + return this->Superclass::RequestInformation(request, + inputVector, outputVector); +} + +//----------------------------------------------------------------------------- +int vtkPNrrdReader::ReadHeader() +{ + if (!this->FileName) + { + vtkErrorMacro(<< "No filename set."); + return 0; + } + + VTK_CREATE(vtkCharArray, headerBuffer); + + // Having a dummy controller means less cases later. + if (!this->Controller) this->Controller = vtkDummyController::New(); + + // Read the header on process 0 and broadcast to everyone else. + if (this->Controller->GetLocalProcessId() == 0) + { + vtkstd::ifstream file(this->FileName, + vtkstd::ios_base::in | vtkstd::ios_base::binary); + // Read in 4 MB. Assuming that the header will be smaller than that. + headerBuffer->SetNumberOfTuples(0x400000); + file.read(headerBuffer->GetPointer(0), 0x400000-1); + vtkIdType buffersize = file.gcount(); + headerBuffer->SetValue(buffersize, '\0'); + headerBuffer->SetNumberOfTuples(buffersize+1); + + // Find a blank line (which signals the end of the header). Be careful + // because line endings can be "\n", "\r\n", or both. It is possible that + // the entire file is the header (happens with "detatched headers"). + char *bufferStart = headerBuffer->GetPointer(0); + char *s = bufferStart; + while ((s = strchr(s+1, '\n')) != NULL) + { + // If you find a double line ending, terminate the string and shorten the + // buffer. + if (s[1] == '\n') + { + s[2] = '\0'; + headerBuffer->SetNumberOfTuples( + static_cast<vtkIdType>(s + 3 - bufferStart)); + break; + } + if ((s[1] == '\r') && (s[2] == '\n')) + { + s[3] = '\0'; + headerBuffer->SetNumberOfTuples( + static_cast<vtkIdType>(s + 4 - bufferStart)); + break; + } + } + } + + this->Controller->Broadcast(headerBuffer, 0); + + return this->ReadHeader(headerBuffer); +} + +//----------------------------------------------------------------------------- +int vtkPNrrdReader::ReadHeader(vtkCharArray *headerBuffer) +{ + this->HeaderSize = headerBuffer->GetNumberOfTuples(); + + vtkstd::string headerStringBuffer = headerBuffer->GetPointer(0); + vtkstd::stringstream header(headerStringBuffer); + + vtkstd::string line; + getline(header, line); + if (line.substr(0, 4) != "NRRD") + { + vtkErrorMacro(<< this->FileName << " is not a nrrd file."); + return 0; + } + + this->DataFiles->Initialize(); + int byteskip = 0; + int numDimensions = 0; + int subDimension = -1; + vtkstd::vector<int> dimSizes; + vtkstd::vector<double> dimSpacing; + this->FileLowerLeft = 1; + while (1) + { + getline(header, line); + if (line.length() < 1) break; + + if (line[0] == '#') + { + // Comment. Ignore. + continue; + } + + size_t delm = line.find(": "); + if (delm != line.npos) + { + // A field/description pair. + vtkstd::string field = line.substr(0, delm); + vtkstd::string description = line.substr(delm+2); + if (field == "dimension") + { + numDimensions = atoi(description.c_str()); + } + else if (field == "sizes") + { + dimSizes.resize(numDimensions); + GetVector(description, dimSizes); + } + else if (field == "spacings") + { + dimSpacing.resize(numDimensions); + GetVector(description, dimSpacing); + } + else if (field == "type") + { + this->DataScalarType = NrrdType2VTKType(description); + if (this->DataScalarType == VTK_VOID) return 0; + // The superclass does this, but I'm not sure it's necessary. + this->GetOutput()->SetScalarType(this->DataScalarType); + } + else if (field == "encoding") + { + if (description != "raw") + { + vtkErrorMacro(<< "Unsupported encoding: " << description); + return 0; + } + } + else if ((field == "data file") || (field == "datafile")) + { + vtkstd::vector<vtkstd::string> filepatterninfo = split(description); + if (filepatterninfo[0] == "LIST") + { + // After LIST there is an optional subdimension (see next case below). + subDimension = ( (filepatterninfo.size() > 1) + ? atoi(filepatterninfo[1].c_str()) : numDimensions ); + + // In this mode files are listed one per line to the end of the file. + while (1) + { + getline(header, line); + trim(line); + if (line.length() < 1) break; + this->DataFiles->InsertNextValue(line); + } + break; + } + else if (filepatterninfo.size() >=4) + { + // description should be "<format> <min> <max> <step> [<subdim>]" + // where <format> is a string to be processed by sprintf and <min>, + // <max>, and <step> form the numbers. <subdim> defines on which + // dimension the files are split up. + vtkstd::string format = filepatterninfo[0]; + int min = atoi(filepatterninfo[1].c_str()); + int max = atoi(filepatterninfo[2].c_str()); + int step = atoi(filepatterninfo[3].c_str()); + subDimension = ( (filepatterninfo.size() > 4) + ? atoi(filepatterninfo[4].c_str()) : numDimensions ); + char *filename = new char[format.size() + 20]; + for (int i = min; i <= max; i += step) + { + sprintf(filename, format.c_str(), i); + this->DataFiles->InsertNextValue(filename); + } + delete[] filename; + } + else + { + // description is simply a filename + this->DataFiles->InsertNextValue(description); + } + } + else if (field == "space") + { + // All spaces are either 3D or 3D with time. + if (description.find("time") != description.npos) + { + vtkErrorMacro(<< "Time in NRRD array not supported (yet)."); + return 0; + } + if ( (description == "left-anterior-superior") + || (description == "LAS") + || (description == "3D-left-handed") ) + { + this->FileLowerLeft = 0; + } + numDimensions = 3; + } + else if (field == "labels") + { + vtkstd::string dataname = description.substr(description.find("\"")+1); + dataname = dataname.substr(0, dataname.find("\"")); + delete[] this->ScalarArrayName; + this->ScalarArrayName = new char[dataname.size()+1]; + strcpy(this->ScalarArrayName, dataname.c_str()); + } + else if (field == "space dimension") + { + numDimensions = atoi(description.c_str()); + } + else if (field == "space origin") + { + vtkstd::vector<double> origin = ParseVector(description); + for (size_t i = 0; (i < 3) && (i < origin.size()); i++) + { + this->DataOrigin[i] = origin[i]; + } + } + else if (field == "space directions") + { + vtkstd::vector<vtkstd::string> directions = split(description); + dimSpacing.resize(0); + for (size_t j = 0; j < directions.size(); j++) + { + if (directions[j] == "none") + { + dimSpacing.push_back(0); + continue; + } + vtkstd::vector<double> dir = ParseVector(directions[j]); + // We don't support orientation, but we do support spacing. + double mag = 0.0; + for (size_t k = 0; k < dir.size(); k++) + { + mag += dir[k]*dir[k]; + } + dimSpacing.push_back(sqrt(mag)); + } + } + else if (field == "endian") + { +#ifdef VTK_WORDS_BIGENDIAN + this->SwapBytes = static_cast<int>(description == "little"); +#else + this->SwapBytes = static_cast<int>(description == "big"); +#endif + } + else if ((field == "line skip") || (field == "lineskip")) + { + if (atoi(description.c_str()) != 0) + { + vtkErrorMacro(<< "line skip not supported"); + return 0; + } + } + else if ((field == "byte skip") || (field == "byteskip")) + { + byteskip = atoi(description.c_str()); + } + else if ( (field == "space units") + || (field == "sample units") || (field == "sampleunits") + || (field == "measurement frame") + || (field == "block size") || (field == "blocksize") + || (field == "content") + || (field == "thicknesses") + || (field == "axis mins") || (field == "axismins") + || (field == "axis maxs") || (field == "axismaxs") + || (field == "centers") || (field == "centerings") + || (field == "units") + || (field == "kinds") + || (field == "min") || (field == "max") + || (field == "old min") || (field == "oldmin") + || (field == "old max") || (field == "oldmax") + || (field == "number") ) + { + // Ignored. + } + else + { + vtkWarningMacro(<< "Unknown field: '" << field << "'"); + } + continue; + } + + delm = line.find(":="); + if (delm != line.npos) + { + // A key/value pair. + continue; + } + } + + // NRRD does not distinguish between vector entries and dimensions. For + // example, RGB tuples are represented by adding a dimension of size 3. + // VTK really needs to know the difference. Here we are going to guess. + // If the fastest changing dimension is 9 or less we consider that a tuple. + // We will also consider any 4th dimension as a tuple. + if ( (dimSizes.size() > 3) + || ((dimSizes.size() > 0) && (dimSizes[0] <= 9)) + || ((dimSpacing.size() > 0) && (dimSpacing[0] == 0)) ) + { + this->NumberOfScalarComponents = dimSizes[0]; + dimSizes.erase(dimSizes.begin()); + if (!dimSpacing.empty()) dimSpacing.erase(dimSpacing.begin()); + subDimension--; + } + else + { + this->NumberOfScalarComponents = 1; + } + + // Record the dimensions. + this->FileDimensionality = dimSizes.size(); + for (unsigned int i = 0; i < 3; i++) + { + this->DataExtent[i*2+0] = 0; + this->DataExtent[i*2+1] = (i < dimSizes.size()) ? dimSizes[i]-1 : 0; + this->DataSpacing[i] = (i < dimSpacing.size()) ? dimSpacing[i] : 1; + } + + if (this->DataFiles->GetNumberOfValues() > 0) + { + if (this->DataFiles->GetNumberOfValues() > 1) + { + this->FileDimensionality--; + if (this->FileDimensionality != 2) + { + vtkErrorMacro(<< "Data split into multiple files is only supported" + << " when each file is 2D (+ an optional vector" + << " dimension)."); + } + if (subDimension != 3) + { + vtkErrorMacro(<< "Data split into multiple files is only supported" + << " when each file is 2D (+ an optional vector" + << " dimension). This means the subdim must be on" + << " that third (or fourth in the case of a vector)" + << " dimension."); + } + } + vtkstd::string parentDir + = vtksys::SystemTools::GetParentDirectory(this->FileName); + for (vtkIdType i = 0; i < this->DataFiles->GetNumberOfValues(); i++) + { + vtkstd::string relativePath = this->DataFiles->GetValue(i); + vtkstd::string fullPath + = vtksys::SystemTools::CollapseFullPath(relativePath.c_str(), + parentDir.c_str()); + this->DataFiles->SetValue(i, fullPath); + vtkWarningMacro(<< "Loading data file: " <<fullPath.c_str()); + } + } + + vtkstd::stringstream message; + this->Print(message); + vtkWarningMacro(<< message.str().c_str()); + + return 1; +} + +//============================================================================= +int vtkPNrrdReader::RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) +{ + // Get rid of superclasses FileNames. We don't support that functionality, + // but we exploit that of the superclass. + if (this->FileNames != NULL) + { + this->FileNames->Delete(); + this->FileNames = NULL; + } + + char *saveFileName = this->FileName; + + if (this->DataFiles->GetNumberOfValues() == 1) + { + this->FileName = const_cast<char *>(this->DataFiles->GetValue(0).c_str()); + } + else if (this->DataFiles->GetNumberOfValues() > 1) + { + this->FileNames = this->DataFiles; + } + + this->Superclass::RequestData(request, inputVector, outputVector); + + this->FileName = saveFileName; + this->FileNames = NULL; + + return 1; +} diff --git a/VTK/Parallel/vtkPNrrdReader.h b/VTK/Parallel/vtkPNrrdReader.h new file mode 100644 index 0000000..7b03270 --- /dev/null +++ b/VTK/Parallel/vtkPNrrdReader.h @@ -0,0 +1,75 @@ +// -*- c++ -*- +/*========================================================================= + + Program: Visualization Toolkit + Module: $RCSfile$ + + Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + See Copyright.txt or http://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ +/*---------------------------------------------------------------------------- + Copyright (c) Sandia Corporation + See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details. +----------------------------------------------------------------------------*/ + +// .NAME vtkPNrrdReader - Read nrrd files efficiently from parallel file systems (and reasonably well elsewhere). +// +// .SECTION Description +// +// vtkPNrrdReader is a subclass of vtkMPIImageReader that will read Nrrd format +// header information of the image before reading the data. This means that the +// reader will automatically set information like file dimensions. +// +// .SECTION Bugs +// +// There are several limitations on what type of nrrd files we can read. This +// reader only supports nrrd files in raw format. Other encodings like ascii +// and hex will result in errors. When reading in detached headers, this only +// supports reading one file that is detached. +// + +#ifndef __vtkPNrrdReader_h +#define __vtkPNrrdReader_h + +#include "vtkMPIImageReader.h" + +class vtkCharArray; + +class VTK_PARALLEL_EXPORT vtkPNrrdReader : public vtkMPIImageReader +{ +public: + vtkTypeRevisionMacro(vtkPNrrdReader, vtkMPIImageReader); + static vtkPNrrdReader *New(); + virtual void PrintSelf(ostream &os, vtkIndent indent); + + virtual int CanReadFile(const char *filename); + +protected: + vtkPNrrdReader(); + ~vtkPNrrdReader(); + + virtual int RequestInformation(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector); + + virtual int RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector); + + virtual int ReadHeader(); + virtual int ReadHeader(vtkCharArray *headerBuffer); + + vtkStringArray *DataFiles; + +private: + vtkPNrrdReader(const vtkPNrrdReader &); // Not implemented. + void operator=(const vtkPNrrdReader &); // Not implemented. +}; + +#endif //__vtkPNrrdReader_h | ||||||||
Relationships | |
Relationships |
Notes | |
(0019273) Ken Moreland (manager) 2010-01-22 13:34 |
The issue is fixed. Because this fix touches some of the base IO and parallel communication files, someone from Kitware should look at these changes and sign off before I check it in. I have assigned the bug to Utkarsh so that he knows he needs to look at the change. He should assign the bug back to me once he reviews the code. The change is attached as a patch to this report. I have also posted the change to a branch of my github repository. The repository can be cloned from the url <git://github.com/kmorel/ParaView.git> [^] and is in the <MPIImageReader> branch. The branch is also summarized at this github.com web site. http://github.com/kmorel/ParaView/tree/MPIImageReader [^] |
(0019729) Ken Moreland (manager) 2010-03-05 12:30 |
/cvsroot/VTKData/VTKData/Baseline/Parallel/ParallelIso.7proc.png,v <-- ParallelIso.7proc.png initial revision: 1.1 /cvsroot/VTKData/VTKData/Data/headsq/quarter.nhdr,v <-- quarter.nhdr initial revision: 1.1 /cvsroot/ParaView3/ParaView3/Applications/ParaView/ParaViewReaders.xml,v <-- Applications/ParaView/ParaViewReaders.xml new revision: 1.6; previous revision: 1.5 /cvsroot/ParaView3/ParaView3/Servers/ServerManager/Resources/readers.xml,v <-- Servers/ServerManager/Resources/readers.xml new revision: 1.197; previous revision: 1.196 /cvsroot/ParaView3/ParaView3/VTK/IO/vtkImageReader.cxx,v <-- VTK/IO/vtkImageReader.cxx new revision: 1.123; previous revision: 1.122 /cvsroot/ParaView3/ParaView3/VTK/IO/vtkImageReader.h,v <-- VTK/IO/vtkImageReader.h new revision: 1.74; previous revision: 1.73 /cvsroot/ParaView3/ParaView3/VTK/Parallel/CMakeLists.txt,v <-- VTK/Parallel/CMakeLists.txt new revision: 1.204; previous revision: 1.203 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPI.h,v <-- VTK/Parallel/vtkMPI.h new revision: 1.7; previous revision: 1.6 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPICommunicator.cxx,v <-- VTK/Parallel/vtkMPICommunicator.cxx new revision: 1.56; previous revision: 1.55 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPICommunicator.h,v <-- VTK/Parallel/vtkMPICommunicator.h new revision: 1.43; previous revision: 1.42 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPIController.cxx,v <-- VTK/Parallel/vtkMPIController.cxx new revision: 1.29; previous revision: 1.28 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPIController.h,v <-- VTK/Parallel/vtkMPIController.h new revision: 1.26; previous revision: 1.25 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPIImageReader.cxx,v <-- VTK/Parallel/vtkMPIImageReader.cxx initial revision: 1.1 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPIImageReader.h,v <-- VTK/Parallel/vtkMPIImageReader.h initial revision: 1.1 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMultiProcessController.cxx,v <-- VTK/Parallel/vtkMultiProcessController.cxx new revision: 1.39; previous revision: 1.38 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMultiProcessController.h,v <-- VTK/Parallel/vtkMultiProcessController.h new revision: 1.45; previous revision: 1.44 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkPNrrdReader.cxx,v <-- VTK/Parallel/vtkPNrrdReader.cxx initial revision: 1.1 /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkPNrrdReader.h,v <-- VTK/Parallel/vtkPNrrdReader.h initial revision: 1.1 /cvsroot/ParaView3/ParaView3/VTK/Parallel/Testing/Cxx/CMakeLists.txt,v <-- VTK/Parallel/Testing/Cxx/CMakeLists.txt new revision: 1.93; previous revision: 1.92 /cvsroot/ParaView3/ParaView3/VTK/Parallel/Testing/Cxx/ExerciseMultiProcessController.cxx,v <-- VTK/Parallel/Testing/Cxx/ExerciseMultiProcessController.cxx new revision: 1.11; previous revision: 1.10 /cvsroot/ParaView3/ParaView3/VTK/Parallel/Testing/Cxx/ParallelIso.cxx,v <-- VTK/Parallel/Testing/Cxx/ParallelIso.cxx new revision: 1.28; previous revision: 1.27 cvs update: warning: `vtkMPIImageReader.h' was lost =================================================================== Checking out vtkMPIImageReader.h RCS: /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPIImageReader.h,v VERS: 1.1 *************** cvs update: warning: `vtkPNrrdReader.h' was lost =================================================================== Checking out vtkPNrrdReader.h RCS: /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkPNrrdReader.h,v VERS: 1.1 *************** cvs update: warning: `vtkPNrrdReader.cxx' was lost =================================================================== Checking out vtkPNrrdReader.cxx RCS: /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkPNrrdReader.cxx,v VERS: 1.1 *************** cvs update: warning: `vtkMPIImageReader.cxx' was lost =================================================================== Checking out vtkMPIImageReader.cxx RCS: /cvsroot/ParaView3/ParaView3/VTK/Parallel/vtkMPIImageReader.cxx,v VERS: 1.1 *************** |
(0019792) Alan Scott (manager) 2010-03-09 22:08 |
I tried this back about a week ago on the small supernova, and it worked fine. Remote servers, Trunk. |
Notes |
Issue History | |||
Date Modified | Username | Field | Change |
2010-01-22 13:24 | Ken Moreland | New Issue | |
2010-01-22 13:24 | Ken Moreland | Status | backlog => tabled |
2010-01-22 13:24 | Ken Moreland | Assigned To | => Ken Moreland |
2010-01-22 13:27 | Ken Moreland | File Added: MPIImageReader.patch | |
2010-01-22 13:34 | Ken Moreland | Note Added: 0019273 | |
2010-01-22 13:34 | Ken Moreland | Assigned To | Ken Moreland => Utkarsh Ayachit |
2010-02-04 13:04 | Utkarsh Ayachit | Priority | normal => immediate |
2010-02-10 09:34 | Utkarsh Ayachit | Assigned To | Utkarsh Ayachit => Berk Geveci |
2010-03-01 09:11 | Utkarsh Ayachit | Assigned To | Berk Geveci => Ken Moreland |
2010-03-05 12:30 | Ken Moreland | Note Added: 0019729 | |
2010-03-05 12:30 | Ken Moreland | Status | tabled => @80@ |
2010-03-05 12:30 | Ken Moreland | Fixed in Version | => 3.8 |
2010-03-05 12:30 | Ken Moreland | Resolution | open => fixed |
2010-03-09 22:08 | Alan Scott | Note Added: 0019792 | |
2010-03-09 22:08 | Alan Scott | Status | @80@ => closed |
2011-06-16 13:10 | Zack Galbreath | Category | Feature Request => Feature |
Issue History |
Copyright © 2000 - 2018 MantisBT Team |