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 @@ + 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 @@ + + + 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. + + + + + + The name of the file to read (or the meta data file that will + point to the actual file). + + + + + + + + + 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(~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(~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(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 +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 + inline T MY_ABS(T x) { return (x < 0) ? -x : x; } + + template + 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 + (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("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(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(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 name = vtkSmartPointer::New() + +#include +#include #include //---------------------------------------------------------------------------- @@ -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 allColors(numProc); + this->AllGather(&localColor, &allColors[0], 1); + + vtkstd::vector allKeys(numProc); + this->AllGather(&localKey, &allKeys[0], 1); + + vtkstd::vector inPartition; + inPartition.assign(numProc, false); + + for (int i = 0; i < numProc; i++) + { + if (inPartition[i]) continue; + int targetColor = allColors[i]; + vtkstd::list 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::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::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 +#include +#include +#include +#include + +#include +#include + +#include "vtkSmartPointer.h" +#define VTK_CREATE(type, name) \ + vtkSmartPointer name = vtkSmartPointer::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 split(vtkstd::string s) +{ + vtkstd::vector 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 &dest) +{ + vtkstd::vector 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 &dest) +{ + vtkstd::vector 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 ParseVector(vtkstd::string s) +{ + vtkstd::vector 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(s + 3 - bufferStart)); + break; + } + if ((s[1] == '\r') && (s[2] == '\n')) + { + s[3] = '\0'; + headerBuffer->SetNumberOfTuples( + static_cast(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 dimSizes; + vtkstd::vector 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 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 " []" + // where is a string to be processed by sprintf and , + // , and form the numbers. 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 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 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 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(description == "little"); +#else + this->SwapBytes = static_cast(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: " <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(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