commit 7d9b3e514ca7d746bbfaf15ec40cbdff606a9743 Author: Davide Punzo Date: Wed Oct 3 11:21:21 2018 +0200 ENH: Fix RAM performance reading/writing nrrd diff --git a/Libs/vtkTeem/vtkTeemNRRDReader.cxx b/Libs/vtkTeem/vtkTeemNRRDReader.cxx index 6ed57223a..2e4299087 100644 --- a/Libs/vtkTeem/vtkTeemNRRDReader.cxx +++ b/Libs/vtkTeem/vtkTeemNRRDReader.cxx @@ -32,22 +32,23 @@ #include "vtkBitArray.h" #include "vtkCharArray.h" #include "vtkDoubleArray.h" +#include "vtkErrorCode.h" #include "vtkFloatArray.h" #include "vtkImageData.h" -#include -#include +#include "vtkInformation.h" +#include "vtkInformationVector.h" #include "vtkIntArray.h" #include "vtkLongArray.h" #include "vtkMath.h" #include "vtkNew.h" #include "vtkObjectFactory.h" #include "vtkShortArray.h" -#include +#include "vtkStreamingDemandDrivenPipeline.h" #include "vtkUnsignedCharArray.h" #include "vtkUnsignedShortArray.h" #include "vtkUnsignedIntArray.h" #include "vtkUnsignedLongArray.h" -#include +#include "vtksys/SystemTools.hxx" // Teem includes #include "teem/ten.h" @@ -66,6 +67,12 @@ vtkTeemNRRDReader::vtkTeemNRRDReader() this->PointDataType = -1; this->DataType = -1; this->NumberOfComponents = -1; +#ifdef NRRD_CHUNK_IO_AVAILABLE + this->ReadImageListAsMultipleImages = false; + this->IsCompressed = -1; + this->NumberOfImages = 1; + this->CurrentImageIndex = 0; +#endif } //---------------------------------------------------------------------------- @@ -227,11 +234,57 @@ int vtkTeemNRRDReader::CanReadFile(const char* filename) supported = false; } + this->IsCompressed = nio->encoding->isCompression; + nrrdNuke(nrrdTemp); nrrdIoStateNix(nio); return supported; } +//---------------------------------------------------------------------------- +bool vtkTeemNRRDReader::ReadImageListAsMultipleImagesOn() +{ + if (this->IsCompressed < 0) + { + vtkErrorMacro(<< "vtkTeemNRRDReader::ReadImageListAsMultipleImagesOn :" + " this method can be used only after invoking vtkTeemNRRDReader::CanReadFile."); + return false; + } + else if (this->IsCompressed) + { + vtkDebugMacro(<< "vtkTeemNRRDReader::ReadImageListAsMultipleImagesOn :" + " reading an image list as multiple image feature" + " is not available for compressed data." + " The data will be read a single multi-component image."); + return false; + } + + this->ReadImageListAsMultipleImages = true; + return true; +} + +//---------------------------------------------------------------------------- +bool vtkTeemNRRDReader::SetReadImageListAsMultipleImages(bool activate) +{ + if (this->IsCompressed < 0) + { + vtkErrorMacro(<< "vtkTeemNRRDReader::SetReadImageListAsMultipleImages :" + " this method can be used only after invoking vtkTeemNRRDReader::CanReadFile."); + return false; + } + else if (activate && this->IsCompressed) + { + vtkDebugMacro(<< "vtkTeemNRRDReader::SetReadImageListAsMultipleImages :" + " reading an image list as multiple image feature" + " is not available for compressed data." + " The data will be read a single multi-component image."); + return false; + } + + this->ReadImageListAsMultipleImages = activate; + return true; +} + //---------------------------------------------------------------------------- bool vtkTeemNRRDReader::GetPointType(Nrrd* nrrdTemp, int& pointDataType, int &numOfComponents) { @@ -325,6 +378,30 @@ bool vtkTeemNRRDReader::GetPointType(Nrrd* nrrdTemp, int& pointDataType, int &nu return false; } +//---------------------------------------------------------------------------- +bool vtkTeemNRRDReader::IsNrrdList(Nrrd *nrrdTemp) +{ + // 4D volume? + if (nrrdTemp->dim != 4) + { + return false; + } + + // 4D volume is a list? + unsigned int rangeAxisKind = nrrdTemp->axis[0].kind; + if (rangeAxisKind == nrrdKindList) + { + return true; + } + rangeAxisKind = nrrdTemp->axis[3].kind; + if (rangeAxisKind == nrrdKindList) + { + return true; + } + + return false; +} + //---------------------------------------------------------------------------- void vtkTeemNRRDReader::ExecuteInformation() { @@ -361,7 +438,8 @@ void vtkTeemNRRDReader::ExecuteInformation() if (nrrdLoad(this->nrrd, this->GetFileName(), nio) != 0) { char *err = biffGetDone(NRRD); - vtkErrorMacro("Error reading " << this->GetFileName() << ": " << err); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation: " + " reading " << this->GetFileName() << ": " << err); free(err); // err points to malloc'd data!! err = NULL; nio = nrrdIoStateNix(nio); @@ -377,12 +455,15 @@ void vtkTeemNRRDReader::ExecuteInformation() if (nrrdTypeBlock == this->nrrd->type) { - vtkErrorMacro("ReadImageInformation: Cannot currently handle nrrdTypeBlock"); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation :" + " Cannot currently handle nrrdTypeBlock"); nio = nrrdIoStateNix(nio); this->ReadStatus = 1; return; } + this->IsCompressed = nio->encoding->isCompression; + if (nio->endian == airEndianLittle) { this->SetDataByteOrderToLittleEndian(); @@ -406,7 +487,8 @@ void vtkTeemNRRDReader::ExecuteInformation() unsigned int domainAxisNum = nrrdDomainAxesGet(this->nrrd, domainAxisIdx); if (this->nrrd->spaceDim && this->nrrd->spaceDim != domainAxisNum) { - vtkErrorMacro("ReadImageInformation: this->nrrd's # independent axes (" + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation:" + " this->nrrd's # independent axes (" << domainAxisNum << ") doesn't match dimension of space" " in which orientation is defined (" << this->nrrd->spaceDim << "); not currently handled"); @@ -419,13 +501,36 @@ void vtkTeemNRRDReader::ExecuteInformation() int numOfComponents = -1; if (!vtkTeemNRRDReader::GetPointType(this->nrrd, pointDataType, numOfComponents)) { - vtkErrorMacro("ReadImageInformation: only 3 spatial dimension and 1 optional range axis is supported"); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation :" + " only 3 spatial dimension and 1 optional" + " range axis is supported"); nio = nrrdIoStateNix(nio); this->ReadStatus = 1; return; } this->SetPointDataType(pointDataType); + +#ifdef NRRD_CHUNK_IO_AVAILABLE + // Check if nrrd file has multiple images and + // ReadImageListAsMultipleImages is on. + // If yes set each component of a list (nrrd) as an individual input connection. + // In this way the reader can provide N 3D volumes instead of a huge 4D volume. + // This provides much better RAM performance. + if (this->GetReadImageListAsMultipleImages() && + vtkTeemNRRDReader::IsNrrdList(this->nrrd) && + !this->IsCompressed) + { + this->SetNumberOfComponents(1); + this->SetNumberOfImages(numOfComponents); + } + else + { + this->SetNumberOfComponents(numOfComponents); + this->SetNumberOfImages(1); + } +#else this->SetNumberOfComponents(numOfComponents); +#endif // Set type information this->SetDataType(this->NrrdToVTKScalarType(this->nrrd->type)); @@ -437,18 +542,18 @@ void vtkTeemNRRDReader::ExecuteInformation() double origin[3] = { 0 }; vtkNew ijkToRasMatrix; if (domainAxisNum > 3) - { - vtkErrorMacro("ReadImageInformation: only up to 3 domain axes are supported"); + { + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation :" + " only up to 3 domain axes are supported"); nio = nrrdIoStateNix(nio); this->ReadStatus = 1; return; - } + } for (unsigned int axii = 0; axii < domainAxisNum; axii++) { unsigned int naxi = domainAxisIdx[axii]; dataExtent[2 * axii] = 0; dataExtent[2 * axii + 1] = static_cast(this->nrrd->axis[naxi].size) - 1; - double spaceDir[NRRD_SPACE_DIM_MAX]; double axisSpacing = 1.0; int spacingStatus = nrrdSpacingCalculate(this->nrrd, naxi, &axisSpacing, spaceDir); @@ -499,10 +604,12 @@ void vtkTeemNRRDReader::ExecuteInformation() break; default: case nrrdSpacingStatusUnknown: - vtkErrorMacro("ReadImageInformation: Error interpreting nrrd spacing (nrrdSpacingStatusUnknown)"); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation : " + "Error interpreting nrrd spacing (nrrdSpacingStatusUnknown)"); break; case nrrdSpacingStatusScalarWithSpace: - vtkErrorMacro("ReadImageInformation: Error interpreting nrrd spacing (nrrdSpacingStatusScalarWithSpace)"); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation : " + "Error interpreting nrrd spacing (nrrdSpacingStatusScalarWithSpace)"); break; } } @@ -557,7 +664,8 @@ void vtkTeemNRRDReader::ExecuteInformation() default: case nrrdOriginStatusUnknown: case nrrdOriginStatusDirection: - vtkErrorMacro("ReadImageInformation: Error interpreting nrrd origin status"); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteInformation :" + " Error interpreting nrrd origin status"); nio = nrrdIoStateNix(nio); this->ReadStatus = 1; break; @@ -683,12 +791,10 @@ vtkImageData *vtkTeemNRRDReader::AllocateOutputData(vtkDataObject *out, vtkInfor // Until I can eliminate the method, I will reexecute the ExecuteInformation // before the execute. this->ExecuteInformation(); - res->SetExtent(this->GetUpdateExtent()); this->AllocatePointData(res, outInfo); return res; - } //---------------------------------------------------------------------------- @@ -697,7 +803,8 @@ void vtkTeemNRRDReader::AllocatePointData(vtkImageData *out, vtkInformation* out // if the scalar type has not been set then we have a problem if (this->DataType == VTK_VOID) { - vtkErrorMacro("Attempt to allocate scalars before scalar type was set!."); + vtkErrorMacro(<< "vtkTeemNRRDReader::AllocatePointData : " + "Attempt to allocate scalars before scalar type was set!."); return; } @@ -718,7 +825,8 @@ void vtkTeemNRRDReader::AllocatePointData(vtkImageData *out, vtkInformation* out pd = out->GetPointData()->GetTensors(); break; default: - vtkErrorMacro("Unknown PointData Type."); + vtkErrorMacro(<< "vtkTeemNRRDReader::AllocatePointData: " + "Unknown PointData Type."); return; } @@ -775,7 +883,8 @@ void vtkTeemNRRDReader::AllocatePointData(vtkImageData *out, vtkInformation* out pd = vtkSmartPointer::New(); break; default: - vtkErrorMacro("Could not allocate data type."); + vtkErrorMacro(<< "vtkTeemNRRDReader::AllocatePointData: " + "Could not allocate data type."); return; } vtkDataObject::SetPointDataActiveScalarInfo(outInfo, this->DataType, this->GetNumberOfComponents()); @@ -802,9 +911,10 @@ void vtkTeemNRRDReader::AllocatePointData(vtkImageData *out, vtkInformation* out out->GetPointData()->SetTensors(pd); break; default: - vtkErrorMacro("Unknown PointData Type."); + vtkErrorMacro(<< "vtkTeemNRRDReader::AllocatePointData :" + "Unknown PointData Type."); return; - } + } } //---------------------------------------------------------------------------- @@ -897,34 +1007,38 @@ int vtkTeemNRRDReader::tenSpaceDirectionReduce(Nrrd *nout, const Nrrd *nin, doub // are assumed to be the same as the file extent/order. void vtkTeemNRRDReader::ExecuteDataWithInformation(vtkDataObject *output, vtkInformation* outInfo) { - if (this->GetOutputInformation(0)) + // Update Extent + if (outInfo) { - this->GetOutputInformation(0)->Set( - vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), - this->GetOutputInformation(0)->Get( - vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()), 6); + outInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), + outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()), 6); } - vtkImageData *imageData = this->AllocateOutputData(output, outInfo); - - if (this->GetFileName() == NULL) +#ifdef NRRD_CHUNK_IO_AVAILABLE + // Check if ReadImageListAsMultipleImages is on + int imageIndex = 0; + int numberOfImages = 1; + if (this->GetReadImageListAsMultipleImages() && + !this->IsCompressed) { - vtkErrorMacro(<< "Either a FileName or FilePrefix must be specified."); - return; + imageIndex = this->GetCurrentImageIndex(); + numberOfImages = this->GetNumberOfImages(); } +#endif - // Read in the this->nrrd. Yes, this means that the header is being read - // twice: once by ExecuteInformation, and once here - if ( nrrdLoad(this->nrrd, this->GetFileName(), NULL) != 0 ) + // Allocate the imageData + vtkImageData *imageData = this->AllocateOutputData(output, outInfo); + if (imageData == NULL) { - char *err = biffGetDone(NRRD); // would be nice to free(err) - vtkErrorMacro("Read: Error reading " << this->GetFileName() << ":\n" << err); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation : " + "imageData not found!"); return; } - if (this->nrrd->data == NULL) + if (this->GetFileName() == NULL) { - vtkErrorMacro(<< "data is null."); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation : " + "Either a FileName or FilePrefix must be specified."); return; } @@ -952,11 +1066,186 @@ void vtkTeemNRRDReader::ExecuteDataWithInformation(vtkDataObject *output, vtkInf } this->ComputeDataIncrements(); + if (ptr == NULL) + { + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation : " + "ptr to imageData not found!"); + return; + } + + // Read in the this->nrrd. +#ifdef NRRD_CHUNK_IO_AVAILABLE + // In the case of ReadImageListAsMultipleImages is On, + // it is assumed that the nrrd file is formatted such as: + // "kinds: domain domain domain list" + // or + // "kinds: list domain domain domain" + long int numberOfImageElements = nrrdElementNumber(this->nrrd) / numberOfImages; + long int numberOfElementsPerChunk = 1E+6; + // Check if the file is compressed. In this case, accessing the data at different location is too slow. + // Therefore, the streaming is disabled and the full nrrd is read at once. + if (this->IsCompressed) + { + numberOfElementsPerChunk = numberOfImageElements; + } + long int numberOfChunks = (long int) floor(numberOfImageElements / numberOfElementsPerChunk); + long int numberOfElementsPerFinalChunk = numberOfImageElements - (numberOfChunks * numberOfElementsPerChunk); + long int elementSize = nrrdElementSize(this->nrrd); + long int chunksize = numberOfElementsPerChunk * elementSize; + long int temporalSize = numberOfImages * elementSize; + + int frameAxis = -1; + if ((this->GetReadImageListAsMultipleImages() + && numberOfImages > 1) || + this->IsCompressed) + { + unsigned int rangeAxisIdx[NRRD_DIM_MAX]; + nrrdRangeAxesGet(this->nrrd, rangeAxisIdx); + frameAxis = rangeAxisIdx[0]; + } + + // Read the chunks. + for (long int chunkIndex = 0; chunkIndex <= numberOfChunks; ++chunkIndex) + { if (chunkIndex == numberOfChunks && numberOfElementsPerFinalChunk < 1) + { + continue; + } + NrrdIoState *nio = nrrdIoStateNew(); + if (nio == NULL) + { + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation: " + "could not allocate NRRD Teem I/O struct"); + return; + } + + if (frameAxis == 0) + { + nio->chunkStartElement = chunkIndex * numberOfElementsPerChunk * numberOfImages + imageIndex; + if (chunkIndex < numberOfChunks) + { + nio->chunkElementCount = numberOfElementsPerChunk * numberOfImages - imageIndex; + } + else + { + nio->chunkElementCount = numberOfElementsPerFinalChunk * numberOfImages - imageIndex; + } + } + else if (frameAxis == 3) + { + nio->chunkStartElement = chunkIndex * numberOfElementsPerChunk + imageIndex * numberOfImageElements; + if (chunkIndex < numberOfChunks) + { + nio->chunkElementCount = numberOfElementsPerChunk; + } + else + { + nio->chunkElementCount = numberOfElementsPerFinalChunk; + } + } + else if (frameAxis == -1) + { + nio->chunkStartElement = chunkIndex * numberOfElementsPerChunk; + if (chunkIndex < numberOfChunks) + { + nio->chunkElementCount = numberOfElementsPerChunk; + } + else + { + nio->chunkElementCount = numberOfElementsPerFinalChunk; + } + } + + if (this->IsCompressed) + { + nio->chunkStartElement = 0; + nio->chunkElementCount = 0; + } + + if (nrrdLoad(this->nrrd, this->GetFileName(), nio) != 0) + { + char *err = biffGetDone(NRRD); // would be nice to free(err) + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation :" + " Error reading " << this->GetFileName() << ":\n" << err); + return; + } + + if (this->nrrd->data == NULL) + { + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation : " + "data is null."); + return; + } + + char *data_c = (char*) this->nrrd->data; + char *ptr_c = (char*) ptr; + long int subsetNumberOfElements; + long int chunkIndexSize = chunkIndex * chunksize; + long int sizeToCopy = chunksize; + if (chunkIndex < numberOfChunks) + { + subsetNumberOfElements = numberOfElementsPerChunk; + } + else + { + subsetNumberOfElements = numberOfElementsPerFinalChunk; + sizeToCopy = numberOfElementsPerFinalChunk * elementSize; + } + + if (frameAxis == 0) + { + for (long int ii = 0; ii < subsetNumberOfElements; ii++) + { + memcpy(ptr_c + (chunkIndexSize + (ii * elementSize)), + data_c + (ii * temporalSize), + elementSize); + } + } + else if (this->IsCompressed && frameAxis == 3) + { + long int oneElmentFullComponentsSize = this->NumberOfComponents * elementSize; + long int componentNumberOfElements = numberOfImageElements / this->NumberOfComponents; + long int oneFullComponentSize = elementSize * componentNumberOfElements; + for (long int jj = 0; jj < this->NumberOfComponents; jj++) + { + long int jjElementSize = jj * elementSize; + long int jjComponentSize = jj * oneFullComponentSize; + for (long int ii = 0; ii < componentNumberOfElements ; ii++) + { + memcpy(ptr_c + (ii * oneElmentFullComponentsSize + jjElementSize), + data_c + (jjComponentSize + ii * elementSize), + elementSize); + } + } + } + else + { + memcpy(ptr_c + chunkIndexSize, data_c, sizeToCopy); + } + + this->nrrd->data = airFree(this->nrrd->data); + delete nio; + } + +#else + if ( nrrdLoad(this->nrrd, this->GetFileName(), NULL) != 0 ) + { + char *err = biffGetDone(NRRD); // would be nice to free(err) + vtkErrorMacro("Read: Error reading " << this->GetFileName() << ":\n" << err); + return; + } + + if (this->nrrd->data == NULL) + { + vtkErrorMacro(<< "data is null."); + return; + } + unsigned int rangeAxisIdx[NRRD_DIM_MAX] = { 0 }; unsigned int rangeAxisNum = nrrdRangeAxesGet(this->nrrd, rangeAxisIdx); if (rangeAxisNum > 1) { - vtkErrorMacro("Read: handling more than one non-scalar axis not currently handled"); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation: " + "handling more than one non-scalar axis not currently handled"); return; } if (rangeAxisNum == 1 && rangeAxisIdx[0] != 0) @@ -973,16 +1262,16 @@ void vtkTeemNRRDReader::ExecuteDataWithInformation(vtkDataObject *output, vtkInf } // The memory size of the input and output of nrrdAxesPermute is // the same; the existing this->nrrd->data is re-used. - if (nrrdCopy(ntmp, this->nrrd) - || nrrdAxesPermute(this->nrrd, ntmp, axmap)) + if (nrrdCopy(ntmp, this->nrrd) || nrrdAxesPermute(this->nrrd, ntmp, axmap)) { char *err = biffGetDone(NRRD); // would be nice to free(err) - vtkErrorMacro("Read: Error permuting independent axis in " << this->GetFileName() << ":\n" << err); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation :" + " Error permuting independent axis in " << this->GetFileName() << ":\n" << err); return; } - nrrdNuke(ntmp); } +#endif // "the famous y-flip": we always flip along the second domain axis //Nrrd *nflip = nrrdNew(); @@ -1015,7 +1304,8 @@ void vtkTeemNRRDReader::ExecuteDataWithInformation(vtkDataObject *output, vtkInf || nrrdPad_nva(this->nrrd, ntmp, minIdx, maxIdx, nrrdBoundaryPad, 1.0)) { char *err = biffGetDone(NRRD); // would be nice to free(err) - vtkErrorMacro("Read: Error padding on conf mask in " << this->GetFileName() << ":\n" << err); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation: " + "Error padding on conf mask in " << this->GetFileName() << ":\n" << err); return; } } @@ -1077,7 +1367,8 @@ void vtkTeemNRRDReader::ExecuteDataWithInformation(vtkDataObject *output, vtkInf if (errorCode) { char *err = biffGetDone(key); // would be nice to free(err) - vtkErrorMacro("Read: Error copying, crapping or cropping:\n" << err); + vtkErrorMacro(<< "vtkTeemNRRDReader::ExecuteDataWithInformation : " + "Error copying, crapping or cropping:\n" << err); return; } nrrdNuke(ntmp); @@ -1089,6 +1380,7 @@ void vtkTeemNRRDReader::ExecuteDataWithInformation(vtkDataObject *output, vtkInf // be called here if it existed. } +#ifndef NRRD_CHUNK_IO_AVAILABLE if (ptr) { memcpy(ptr, this->nrrd->data, nrrdElementSize(this->nrrd)*nrrdElementNumber(this->nrrd)); @@ -1096,6 +1388,7 @@ void vtkTeemNRRDReader::ExecuteDataWithInformation(vtkDataObject *output, vtkInf // release the memory while keeping the struct nrrdEmpty(this->nrrd); +#endif } //---------------------------------------------------------------------------- diff --git a/Libs/vtkTeem/vtkTeemNRRDReader.h b/Libs/vtkTeem/vtkTeemNRRDReader.h index b87092007..2c9544f26 100644 --- a/Libs/vtkTeem/vtkTeemNRRDReader.h +++ b/Libs/vtkTeem/vtkTeemNRRDReader.h @@ -82,6 +82,7 @@ public: /// supports spaces in key names. const std::map GetHeaderKeysMap(); + /// /// Get a value given a key in the header const char* GetHeaderValue(const char *key); @@ -128,10 +129,54 @@ public: vtkGetMacro(DataType,int); /// - //Number of components + /// Number of components vtkSetMacro(NumberOfComponents,int); vtkGetMacro(NumberOfComponents,int); +#ifdef NRRD_CHUNK_IO_AVAILABLE + /// + /// If NumberOfComponents > 1, + /// the nrrd file will be read as a list of images. + /// The method Update will return only one component at time (as a vtkImageData). + /// The component that will be read is given by CurrentImageIndex + /// + /// These methods have to be used after CanReadFile. + /// For compressed file this option is not available + /// (the data will be read a single multi-component image) + /// and in this case the value will be forced to Off. + /// + /// It is assumed that the nrrd file is formatted such as: + /// "kinds: domain domain domain list" + /// or + /// "kinds: list domain domain domain" + /// + /// \sa NumberOfImages, CurrentImageIndex + bool ReadImageListAsMultipleImagesOn(); + + /// + /// If NumberOfComponents > 1, + /// the nrrd file will be read as a single multi-component + void ReadImageListAsMultipleImagesOff() + { + ReadImageListAsMultipleImages = false; + } + + /// + /// Set/Get ReadImageListAsMultipleImages (true/false). + /// \sa SetReadImageListAsMultipleImages(), GetReadImageListAsMultipleImages() + bool SetReadImageListAsMultipleImages(bool); + vtkGetMacro(ReadImageListAsMultipleImages,bool); + + /// + /// Number of images + vtkSetMacro(NumberOfImages,int); + vtkGetMacro(NumberOfImages,int); + + /// + /// List image index to be read/written + vtkSetMacro(CurrentImageIndex,int); + vtkGetMacro(CurrentImageIndex,int); +#endif /// /// Use image origin from the file @@ -231,6 +276,7 @@ public: break; } } + virtual vtkImageData * AllocateOutputData(vtkDataObject *out, vtkInformation* outInfo) VTK_OVERRIDE; virtual void AllocateOutputData(vtkImageData *out, vtkInformation* outInfo, int *uExtent) VTK_OVERRIDE { Superclass::AllocateOutputData(out, outInfo, uExtent); } @@ -241,6 +287,7 @@ protected: ~vtkTeemNRRDReader(); static bool GetPointType(Nrrd* nrrdTemp, int& pointDataType, int &numOfComponents); + static bool IsNrrdList(Nrrd* nrrdTemp); vtkSmartPointer RasToIjkMatrix; vtkSmartPointer MeasurementFrameMatrix; @@ -256,6 +303,12 @@ protected: int DataType; int NumberOfComponents; bool UseNativeOrigin; +#ifdef NRRD_CHUNK_IO_AVAILABLE + int NumberOfImages; + int CurrentImageIndex; + bool ReadImageListAsMultipleImages; + int IsCompressed; +#endif std::map HeaderKeyValue; std::string HeaderKeys; // buffer for returning key list diff --git a/Libs/vtkTeem/vtkTeemNRRDWriter.cxx b/Libs/vtkTeem/vtkTeemNRRDWriter.cxx index fe11d51c6..36ee4cf6a 100644 --- a/Libs/vtkTeem/vtkTeemNRRDWriter.cxx +++ b/Libs/vtkTeem/vtkTeemNRRDWriter.cxx @@ -1,8 +1,5 @@ #include - #include "vtkTeemNRRDWriter.h" - - #include "vtkImageData.h" #include "vtkPointData.h" #include "vtkObjectFactory.h" @@ -14,7 +11,6 @@ #include "itkNumberToString.h" - class AttributeMapType: public std::map {}; class AxisInfoMapType : public std::map {}; @@ -31,13 +27,20 @@ vtkTeemNRRDWriter::vtkTeemNRRDWriter() this->UseCompression = 1; // use default CompressionLevel this->CompressionLevel = -1; - this->DiffusionWeightedData = 0; + this->DiffusionWeigthedData = 0; this->FileType = VTK_BINARY; this->WriteErrorOff(); this->Attributes = new AttributeMapType; this->AxisLabels = new AxisInfoMapType; this->AxisUnits = new AxisInfoMapType; this->VectorAxisKind = nrrdKindUnknown; +#ifdef NRRD_CHUNK_IO_AVAILABLE + this->WriteMultipleImagesAsImageList = false; + this->NumberOfImages = 1; + this->CurrentImageIndex = 0; +#endif + this->nrrd = NULL; + this->nio = NULL; } //---------------------------------------------------------------------------- @@ -81,9 +84,15 @@ int vtkTeemNRRDWriter::FillInputPortInformation( // Writes all the data from the input. void vtkTeemNRRDWriter::vtkImageDataInfoToNrrdInfo(vtkImageData *in, int &kind, size_t &numComp, int &vtkType, void **buffer) { + if (in == NULL) + { + vtkErrorMacro(<< "vtkTeemNRRDWriter::vtkImageDataInfoToNrrdInfo: " + "no input vtkImageData found."); + return; + } vtkDataArray *array; - this->DiffusionWeightedData = 0; + this->DiffusionWeigthedData = 0; if ((array = static_cast (in->GetPointData()->GetScalars()))) { numComp = array->GetNumberOfComponents(); @@ -110,7 +119,7 @@ void vtkTeemNRRDWriter::vtkImageDataInfoToNrrdInfo(vtkImageData *in, int &kind, if (numGrad == numBValues && numGrad == numComp && numGrad>6) { kind = nrrdKindList; - this->DiffusionWeightedData = 1; + this->DiffusionWeigthedData = 1; } else { @@ -145,6 +154,48 @@ void vtkTeemNRRDWriter::vtkImageDataInfoToNrrdInfo(vtkImageData *in, int &kind, } } +//---------------------------------------------------------------------------- +void vtkTeemNRRDWriter::updateNRRDDatapointer(vtkImageData *in) +{ + if (in == NULL) + { + vtkErrorMacro(<< "vtkTeemNRRDWriter::updateNRRDDatapointer: " + "no input vtkImageData found."); + this->nrrd->data = NULL; + return; + } + + void *buffer; + vtkDataArray *array; + + if ((array = static_cast (in->GetPointData()->GetScalars()))) + { + buffer = array->GetVoidPointer(0); + } + else if ((array = static_cast (in->GetPointData()->GetVectors()))) + { + buffer = array->GetVoidPointer(0); + } + else if ((array = static_cast (in->GetPointData()->GetNormals()))) + { + buffer = array->GetVoidPointer(0); + } + else if ((array = static_cast (in->GetPointData()->GetTensors()))) + { + buffer = array->GetVoidPointer(0); + } + else + { + vtkErrorMacro(<< "vtkTeemNRRDWriter::updateNRRDDatapointer: " + "the type of the input vtkImageData is not compatible."); + this->nrrd->data = NULL; + return; + } + + this->nrrd->data = buffer; +} + +//---------------------------------------------------------------------------- int vtkTeemNRRDWriter::VTKToNrrdPixelType( const int vtkPixelType ) { switch( vtkPixelType ) @@ -186,6 +237,7 @@ int vtkTeemNRRDWriter::VTKToNrrdPixelType( const int vtkPixelType ) } } +//---------------------------------------------------------------------------- void* vtkTeemNRRDWriter::MakeNRRD() { Nrrd *nrrd = nrrdNew(); @@ -197,16 +249,16 @@ void* vtkTeemNRRDWriter::MakeNRRD() void *buffer; int vtkType; - // Fill in image information. - - //vtkImageData *input = this->GetInput(); - // Find Pixel type from data and select a buffer. this->vtkImageDataInfoToNrrdInfo(this->GetInput(),kind[0],size[0],vtkType, &buffer); spaceDim = 3; // VTK is always 3D volumes. +#ifdef NRRD_CHUNK_IO_AVAILABLE + if (size[0] > 1 || this->GetWriteMultipleImagesAsImageList()) +#else if (size[0] > 1) - { +#endif + { // the range axis has no space direction for (unsigned int saxi=0; saxi < spaceDim; saxi++) { @@ -221,17 +273,37 @@ void* vtkTeemNRRDWriter::MakeNRRD() nrrdDim = baseDim + spaceDim; unsigned int axi; - for (axi=0; axi < spaceDim; axi++) + for (axi = 0; axi < spaceDim; axi++) { +#ifndef NRRD_CHUNK_IO_AVAILABLE size[axi+baseDim] = this->GetInput()->GetDimensions()[axi]; kind[axi+baseDim] = nrrdKindDomain; - origin[axi] = this->IJKToRASMatrix->GetElement((int) axi,3); +#else + size[axi] = this->GetInput()->GetDimensions()[axi]; + kind[axi] = nrrdKindDomain; +#endif + origin[axi] = this->IJKToRASMatrix->GetElement((int) axi, 3); //double spacing = this->GetInput()->GetSpacing()[axi]; for (unsigned int saxi=0; saxi < spaceDim; saxi++) { +#ifndef NRRD_CHUNK_IO_AVAILABLE spaceDir[axi+baseDim][saxi] = this->IJKToRASMatrix->GetElement(saxi,axi); +#else + spaceDir[axi][saxi] = this->IJKToRASMatrix->GetElement(saxi,axi); +#endif } } +#ifdef NRRD_CHUNK_IO_AVAILABLE + if (this->GetWriteMultipleImagesAsImageList()) + { + size[3] = this->GetNumberOfImages(); + kind[3] = nrrdKindList; + for (unsigned int saxi=0; saxi < spaceDim; saxi++) + { + spaceDir[3][saxi] = AIR_NAN; + } + } +#endif if (nrrdWrap_nva(nrrd, const_cast (buffer), this->VTKToNrrdPixelType( vtkType ), @@ -240,8 +312,8 @@ void* vtkTeemNRRDWriter::MakeNRRD() || nrrdSpaceOriginSet(nrrd, origin)) { char *err = biffGetDone(NRRD); // would be nice to free(err) - vtkErrorMacro("Write: Error wrapping nrrd for " - << this->GetFileName() << ":\n" << err); + vtkErrorMacro(<< "vtkTeemNRRDWriter::MakeNRRD : Error wrapping nrrd for " + << this->GetFileName() << ":\n" << err); // Free the nrrd struct but don't touch nrrd->data nrrd = nrrdNix(nrrd); this->WriteErrorOn(); @@ -312,7 +384,7 @@ void* vtkTeemNRRDWriter::MakeNRRD() itk::NumberToString DoubleConvert; // 2. Take care about diffusion data - if (this->DiffusionWeightedData) + if (this->DiffusionWeigthedData) { unsigned int numGrad = this->DiffusionGradients->GetNumberOfTuples(); unsigned int numBValues = this->BValues->GetNumberOfTuples(); @@ -356,7 +428,19 @@ void* vtkTeemNRRDWriter::MakeNRRD() } } return nrrd; - } +} + +//---------------------------------------------------------------------------- +Nrrd *vtkTeemNRRDWriter::GetNRRDTeem() +{ + return this->nrrd; +} + +//---------------------------------------------------------------------------- +NrrdIoState *vtkTeemNRRDWriter::GetNRRDIoTeem() +{ + return this->nio; +} //---------------------------------------------------------------------------- // Writes all the data from the input. @@ -365,58 +449,94 @@ void vtkTeemNRRDWriter::WriteData() this->WriteErrorOff(); if (this->GetFileName() == NULL) { - vtkErrorMacro("FileName has not been set. Cannot save file"); + vtkErrorMacro(<< "vtkTeemNRRDWriter::WriteData : " + "FileName has not been set. Cannot save file"); this->WriteErrorOn(); return; } - Nrrd* nrrd = (Nrrd*)this->MakeNRRD(); - if (nrrd == NULL) + // Check if GetWriteMultipleImagesAsImageList is on +#ifdef NRRD_CHUNK_IO_AVAILABLE + int imageIndex = 0; + int numberOfImages = 1; + if (this->GetWriteMultipleImagesAsImageList()) { - return; + imageIndex = this->GetCurrentImageIndex(); + numberOfImages = this->GetNumberOfImages(); + } +#endif + + if (this->nrrd == NULL || this->nio == NULL) + { + this->InitializeNRRDTeem(); } - NrrdIoState *nio = nrrdIoStateNew(); + this->updateNRRDDatapointer(this->GetInput()); // set encoding for data: compressed (raw), (uncompressed) raw, or ascii - if ( this->GetUseCompression() && nrrdEncodingGzip->available() ) + if (this->GetUseCompression() && nrrdEncodingGzip->available()) { // this is necessarily gzip-compressed *raw* data - nio->encoding = nrrdEncodingGzip; - nio->zlibLevel = this->CompressionLevel; + this->nio->encoding = nrrdEncodingGzip; + this->nio->zlibLevel = this->CompressionLevel; } else { int fileType = this->GetFileType(); - switch ( fileType ) + switch (fileType) { default: case VTK_BINARY: - nio->encoding = nrrdEncodingRaw; + this->nio->encoding = nrrdEncodingRaw; break; case VTK_ASCII: - nio->encoding = nrrdEncodingAscii; + this->nio->encoding = nrrdEncodingAscii; break; } } // set endianness as unknown of output - nio->endian = airEndianUnknown; + this->nio->endian = airEndianUnknown; + +#ifdef NRRD_CHUNK_IO_AVAILABLE + if (this->GetWriteMultipleImagesAsImageList()) + { + long int numberOfImageElements = nrrdElementNumber(this->nrrd) / numberOfImages; + this->nio->chunkElementCount = numberOfImageElements; + this->nio->chunkStartElement = numberOfImageElements * imageIndex; + if (imageIndex == 0) + { + this->nio->keepNrrdDataFileOpen = 0; + } + else + { + this->nio->keepNrrdDataFileOpen = 1; + } + } +#endif // Write the nrrd to file. - if (nrrdSave(this->GetFileName(), nrrd, nio)) + if (nrrdSave(this->GetFileName(), this->nrrd, this->nio)) { char *err = biffGetDone(NRRD); // would be nice to free(err) - vtkErrorMacro("Write: Error writing " - << this->GetFileName() << ":\n" << err); + vtkErrorMacro(<< "vtkTeemNRRDWriter::WriteData : Error writing " + << this->GetFileName() << ":\n" << err); this->WriteErrorOn(); } - // Free the nrrd struct but don't touch nrrd->data - nrrd = nrrdNix(nrrd); - nio = nrrdIoStateNix(nio); + +#ifdef NRRD_CHUNK_IO_AVAILABLE + if (!this->GetWriteMultipleImagesAsImageList()) + { + // Free the nrrd struct but don't touch nrrd->data + this->nrrd = nrrdNix(this->nrrd); + this->nio = nrrdIoStateNix(this->nio); + } +#endif + return; } +//---------------------------------------------------------------------------- void vtkTeemNRRDWriter::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os,indent); @@ -427,6 +547,7 @@ void vtkTeemNRRDWriter::PrintSelf(ostream& os, vtkIndent indent) this->MeasurementFrameMatrix->PrintSelf(os,indent); } +//---------------------------------------------------------------------------- void vtkTeemNRRDWriter::SetAttribute(const std::string& name, const std::string& value) { if (!this->Attributes) @@ -437,6 +558,18 @@ void vtkTeemNRRDWriter::SetAttribute(const std::string& name, const std::string& (*this->Attributes)[name] = value; } +//---------------------------------------------------------------------------- +std::string vtkTeemNRRDWriter::GetAttribute(const std::string& name) +{ + if (!this->Attributes) + { + return ""; + } + + return (*this->Attributes)[name]; +} + +//---------------------------------------------------------------------------- void vtkTeemNRRDWriter::SetAxisLabel(unsigned int axis, const char* label) { if (!this->AxisLabels) @@ -453,6 +586,7 @@ void vtkTeemNRRDWriter::SetAxisLabel(unsigned int axis, const char* label) } } +//---------------------------------------------------------------------------- void vtkTeemNRRDWriter::SetAxisUnit(unsigned int axis, const char* unit) { if (!this->AxisUnits) @@ -469,7 +603,26 @@ void vtkTeemNRRDWriter::SetAxisUnit(unsigned int axis, const char* unit) } } +//---------------------------------------------------------------------------- void vtkTeemNRRDWriter::SetVectorAxisKind(int kind) { this->VectorAxisKind = kind; } + +//---------------------------------------------------------------------------- +void vtkTeemNRRDWriter::InitializeNRRDTeem() +{ + this->nrrd = (Nrrd*)this->MakeNRRD(); + if (this->nrrd == NULL) + { + vtkErrorMacro(<< "vtkTeemNRRDWriter::InitializeNRRDTeem : " + "failed to instantiate nrrd teem pointer.") + } + + this->nio = nrrdIoStateNew(); + if (this->nio == NULL) + { + vtkErrorMacro(<< "vtkTeemNRRDWriter::InitializeNRRDTeem : " + "failed to instantiate nrrdIo teem pointer.") + } +} diff --git a/Libs/vtkTeem/vtkTeemNRRDWriter.h b/Libs/vtkTeem/vtkTeemNRRDWriter.h index b732c11ba..538fdf0c7 100644 --- a/Libs/vtkTeem/vtkTeemNRRDWriter.h +++ b/Libs/vtkTeem/vtkTeemNRRDWriter.h @@ -63,6 +63,46 @@ public: void SetFileTypeToASCII() {this->SetFileType(VTK_ASCII);}; void SetFileTypeToBinary() {this->SetFileType(VTK_BINARY);}; +#ifdef NRRD_CHUNK_IO_AVAILABLE + /// + /// the nrrd file will be written as a list of images. + /// The method vtkWriter::Write will write only one component at time on file. + /// The volume that will be written on file is given by CurrentImageIndex + /// + /// It is assumed that the nrrd file will be formatted such as: + /// "kinds: domain domain domain list" + /// + /// \sa NumberOfImages, CurrentImageIndex + void WriteMultipleImagesAsImageListsOn() + { + WriteMultipleImagesAsImageList = true; + } + + /// + /// the volume to be saved is not a multi-component image list + void WriteMultipleImagesAsImageListOff() + { + WriteMultipleImagesAsImageList = false; + } + + /// + /// Set/Get ReadImageListAsMultipleImages (true/false). + /// \sa SetReadImageListAsMultipleImages(), GetReadImageListAsMultipleImages() + vtkSetMacro(WriteMultipleImagesAsImageList,bool); + vtkGetMacro(WriteMultipleImagesAsImageList,bool); + + /// + /// Number of images + vtkSetMacro(NumberOfImages,int); + vtkGetMacro(NumberOfImages,int); + + /// + /// Sequences frame to be read/written. + /// The first image index must be zero + vtkSetMacro(CurrentImageIndex,int); + vtkGetMacro(CurrentImageIndex,int); +#endif + vtkBooleanMacro(WriteError, int); vtkSetMacro(WriteError, int); vtkGetMacro(WriteError, int); @@ -71,6 +111,10 @@ public: /// file on write void SetAttribute(const std::string& name, const std::string& value); + /// Method to get an attribute that will be passed into the NRRD + /// file on write + std::string GetAttribute(const std::string& name); + /// Method to set label for each axis void SetAxisLabel(unsigned int axis, const char* label); @@ -85,6 +129,12 @@ public: /// Utility function to return image as a Nrrd* void* MakeNRRD(); + /// Get the teem nrrd pointer + Nrrd* GetNRRDTeem(); + + /// Get the teem nrrdIo state pointer + NrrdIoState* GetNRRDIoTeem(); + protected: vtkTeemNRRDWriter(); ~vtkTeemNRRDWriter(); @@ -95,12 +145,18 @@ protected: /// Write method. It is called by vtkWriter::Write(); void WriteData() VTK_OVERRIDE; + /// Initialize the teem nrrd pointer + void InitializeNRRDTeem(); + /// /// Flag to set to on when a write error occurred int WriteError; char *FileName; + Nrrd* nrrd; + NrrdIoState *nio; + vtkDoubleArray* BValues; vtkDoubleArray* DiffusionGradients; @@ -109,6 +165,11 @@ protected: int UseCompression; int CompressionLevel; +#ifdef NRRD_CHUNK_IO_AVAILABLE + int NumberOfImages; + int CurrentImageIndex; + bool WriteMultipleImagesAsImageList; +#endif int FileType; AttributeMapType *Attributes; @@ -120,8 +181,9 @@ private: vtkTeemNRRDWriter(const vtkTeemNRRDWriter&); /// Not implemented. void operator=(const vtkTeemNRRDWriter&); /// Not implemented. void vtkImageDataInfoToNrrdInfo(vtkImageData *in, int &nrrdKind, size_t &numComp, int &vtkType, void **buffer); + void updateNRRDDatapointer(vtkImageData *in); int VTKToNrrdPixelType( const int vtkPixelType ); - int DiffusionWeightedData; + int DiffusionWeigthedData; }; #endif