Commit 90dfa378 authored by David Gobbi's avatar David Gobbi

Allow NIFTI reader and writer to handle planar RGB files.

The NIFTI format uses packed RGB, but the older Analyze format, in its
earliest incarnations, used planar RGB.  For the sake of compatibility
with these old Analyze files, or with NIFTI files that were created from
them without a planar-to-packed conversion, a PlanarRGB option has been
added to the vtkNIFTIImageReader and the vtkNIFTIImageWriter.
parent ecff3add
......@@ -39,9 +39,10 @@ then re-reading it to ensure that the contents are identical.
#include <string>
static const char *testfiles[6][2] = {
static const char *testfiles[7][2] = {
{ "Data/minimal.nii.gz", "out_minimal.nii.gz" },
{ "Data/minimal.img.gz", "out_minimal.hdr" },
{ "Data/planar_rgb.nii.gz", "out_planar_rgb.nii" },
{ "Data/nifti_rgb.nii.gz", "out_nifti_rgb.nii" },
{ "Data/filtered_func_data.nii.gz", "out_filtered_func_data.nii.gz" },
{ "Data/minimal.hdr.gz", "out_minimal_2.nii" },
......@@ -123,7 +124,8 @@ static void TestDisplay(vtkRenderWindow *renwin, const char *infile)
};
static double TestReadWriteRead(
const char *infile, const char *infile2, const char *outfile)
const char *infile, const char *infile2, const char *outfile,
bool planarRGB)
{
// read a NIFTI file
vtkNew<vtkNIFTIImageReader> reader;
......@@ -139,6 +141,10 @@ static double TestReadWriteRead(
reader->SetFileNames(filenames.GetPointer());
}
reader->TimeAsVectorOn();
if (planarRGB)
{
reader->PlanarRGBOn();
}
reader->Update();
vtkNew<vtkNIFTIImageWriter> writer;
......@@ -163,6 +169,10 @@ static double TestReadWriteRead(
{
writer->SetSFormMatrix(reader->GetQFormMatrix());
}
if (planarRGB)
{
writer->PlanarRGBOn();
}
writer->Write();
// to exercise PrintSelf
......@@ -174,6 +184,10 @@ static double TestReadWriteRead(
vtkNew<vtkNIFTIImageReader> reader2;
reader2->SetFileName(outfile);
reader2->TimeAsVectorOn();
if (planarRGB)
{
reader2->PlanarRGBOn();
}
reader2->Update();
// the images should be identical
......@@ -323,10 +337,11 @@ int TestNIFTIReaderWriter(int argc, char *argv[])
char *infile2 = 0;
char *infile =
vtkTestUtilities::ExpandDataFileName(argc, argv, testfiles[i][0]);
if (i == 4)
bool planarRGB = (i == 2);
if (i == 5)
{
infile2 =
vtkTestUtilities::ExpandDataFileName(argc, argv, testfiles[5][0]);
vtkTestUtilities::ExpandDataFileName(argc, argv, testfiles[6][0]);
}
if (!infile)
{
......@@ -356,7 +371,8 @@ int TestNIFTIReaderWriter(int argc, char *argv[])
return 1;
}
double err = TestReadWriteRead(infile, infile2, outpath.c_str());
double err = TestReadWriteRead(infile, infile2, outpath.c_str(),
planarRGB);
if (err != 0.0)
{
cerr << "Input " << infile << " differs from output " << outpath.c_str()
......
......@@ -63,6 +63,7 @@ vtkNIFTIImageReader::vtkNIFTIImageReader()
this->QFormMatrix = 0;
this->SFormMatrix = 0;
this->NIFTIHeader = 0;
this->PlanarRGB = false;
}
//----------------------------------------------------------------------------
......@@ -222,6 +223,7 @@ void vtkNIFTIImageReader::PrintSelf(ostream& os, vtkIndent indent)
}
os << indent << "NIFTIHeader:" << (this->NIFTIHeader ? "\n" : " (none)\n");
os << indent << "PlanarRGB: " << (this->PlanarRGB ? "On\n" : "Off\n");
}
//----------------------------------------------------------------------------
......@@ -1149,6 +1151,11 @@ int vtkNIFTIImageReader::RequestData(
return 0;
}
// check if planar RGB is applicable (Analyze only)
bool planarRGB = (this->PlanarRGB &&
(this->NIFTIHeader->GetDataType() == NIFTI_TYPE_RGB24 ||
this->NIFTIHeader->GetDataType() == NIFTI_TYPE_RGBA32));
int swapBytes = this->GetSwapBytes();
int scalarSize = data->GetScalarSize();
int numComponents = data->GetNumberOfScalarComponents();
......@@ -1165,6 +1172,7 @@ int vtkNIFTIImageReader::RequestData(
z_off_t fileVoxelIncr = scalarSize*numComponents/vectorDim;
z_off_t fileRowIncr = fileVoxelIncr*this->Dim[1];
z_off_t filePlaneIncr = fileRowIncr*this->Dim[2];
z_off_t fileSliceIncr = fileRowIncr*this->Dim[2];
z_off_t fileTimeIncr = fileSliceIncr*this->Dim[3];
z_off_t fileVectorIncr = fileTimeIncr*this->Dim[4];
......@@ -1173,9 +1181,19 @@ int vtkNIFTIImageReader::RequestData(
fileVectorIncr = fileTimeIncr;
}
// planar RGB requires different increments
int planarSize = 1; // if > 1, indicates planar RGB
if (planarRGB)
{
planarSize = numComponents/vectorDim;
fileVoxelIncr = scalarSize;
fileRowIncr = fileVoxelIncr*this->Dim[1];
filePlaneIncr = fileRowIncr*this->Dim[2];
}
// add a buffer for planar-vector to packed-vector conversion
unsigned char *rowBuffer = 0;
if (vectorDim > 1)
if (vectorDim > 1 || planarRGB)
{
rowBuffer = new unsigned char[outSizeX*fileVoxelIncr];
}
......@@ -1192,11 +1210,23 @@ int vtkNIFTIImageReader::RequestData(
dataPtr += sliceOffset*(outSizeZ - 1);
}
// special increment to handle planar RGB
vtkIdType planarOffset = 0;
vtkIdType planarEndOffset = 0;
if (planarRGB)
{
planarOffset = scalarSize*numComponents;
planarOffset *= outSizeX;
planarOffset *= outSizeY;
planarOffset -= scalarSize;
planarEndOffset = planarOffset - scalarSize*(planarSize - 1);
}
// report progress every 2% of the way to completion
this->InvokeEvent(vtkCommand::StartEvent);
this->UpdateProgress(0.0);
vtkIdType target =
static_cast<vtkIdType>(0.02*outSizeY*outSizeZ*vectorDim) + 1;
static_cast<vtkIdType>(0.02*planarSize*outSizeY*outSizeZ*vectorDim) + 1;
vtkIdType count = 0;
// seek to the start of the data
......@@ -1207,10 +1237,11 @@ int vtkNIFTIImageReader::RequestData(
// read the data one row at a time, do planar-to-packed conversion
// of vector components if NIFTI file has a vector dimension
int rowSize = numComponents/vectorDim*outSizeX;
int rowSize = fileVoxelIncr/scalarSize*outSizeX;
int t = 0; // counter for time
int c = 0; // counter for vector components
int j = 0; // counter for rows
int p = 0; // counter for planes (planar RGB)
int k = 0; // counter for slices
unsigned char *ptr = dataPtr;
......@@ -1232,7 +1263,7 @@ int vtkNIFTIImageReader::RequestData(
}
}
if (vectorDim == 1)
if (vectorDim == 1 && !planarRGB)
{
// read directly into the output instead of into a buffer
rowBuffer = ptr;
......@@ -1254,7 +1285,7 @@ int vtkNIFTIImageReader::RequestData(
vtkByteSwap::SwapVoidRange(rowBuffer, rowSize, scalarSize);
}
if (vectorDim == 1)
if (vectorDim == 1 && !planarRGB)
{
// advance the pointer to the next row
ptr += outSizeX*numComponents*scalarSize;
......@@ -1286,36 +1317,44 @@ int vtkNIFTIImageReader::RequestData(
if (++j == outSizeY)
{
j = 0;
offset += fileSliceIncr - outSizeY*fileRowIncr;
ptr -= 2*sliceOffset; // for reverse slice order
if (++k == outSizeZ)
offset += filePlaneIncr - outSizeY*fileRowIncr;
// back up for next plane (R, G, or B) if planar mode
ptr -= planarOffset;
if (++p == planarSize)
{
k = 0;
offset += fileVectorIncr - outSizeZ*fileSliceIncr;
if (++t == timeDim)
{
t = 0;
}
if (++c == vectorDim)
{
break;
}
// back up the ptr to the beginning of the image,
// then increment to the next vector component
ptr = dataPtr + c*fileVoxelIncr;
if (this->TimeAsVector)
p = 0;
ptr += planarEndOffset; // advance to start of next slice
ptr -= 2*sliceOffset; // for reverse slice order
if (++k == outSizeZ)
{
// if timeDim is included in the vectorDim (and hence in the
// VTK scalar components) then we have to make sure that
// the vector components are packed before the time steps
ptr = dataPtr + (c + t*(vectorDim - 1))/timeDim*fileVoxelIncr;
k = 0;
offset += fileVectorIncr - outSizeZ*fileSliceIncr;
if (++t == timeDim)
{
t = 0;
}
if (++c == vectorDim)
{
break;
}
// back up the ptr to the beginning of the image,
// then increment to the next vector component
ptr = dataPtr + c*fileVoxelIncr*planarSize;
if (this->TimeAsVector)
{
// if timeDim is included in the vectorDim (and hence in the
// VTK scalar components) then we have to make sure that
// the vector components are packed before the time steps
ptr = dataPtr + (c + t*(vectorDim - 1))/timeDim*
fileVoxelIncr*planarSize;
}
}
}
}
}
if (vectorDim > 1)
if (vectorDim > 1 || planarRGB)
{
delete [] rowBuffer;
}
......
......@@ -91,6 +91,16 @@ public:
double GetRescaleSlope() { return this->RescaleSlope; }
double GetRescaleIntercept() { return this->RescaleIntercept; }
// Description:
// Read planar RGB (separate R, G, and B planes), rather than packed RGB.
// The NIFTI format should always use packed RGB. The Analyze format,
// however, was used to store both planar RGB and packed RGB depending
// on the software, without any indication in the header about which
// convention was being used. Use this if you have a planar RGB file.
vtkGetMacro(PlanarRGB, bool);
vtkSetMacro(PlanarRGB, bool);
vtkBooleanMacro(PlanarRGB, bool);
// Description:
// QFac gives the slice order in the NIFTI file versus the VTK image.
// If QFac is -1, then the VTK slice index J is related to the NIFTI
......@@ -204,6 +214,10 @@ protected:
// A copy of the header from the file that was most recently read.
vtkNIFTIImageHeader *NIFTIHeader;
// Description:
// Use planar RGB instead of the default (packed).
bool PlanarRGB;
private:
vtkNIFTIImageReader(const vtkNIFTIImageReader&); // Not implemented.
void operator=(const vtkNIFTIImageReader&); // Not implemented.
......
......@@ -71,6 +71,8 @@ vtkNIFTIImageWriter::vtkNIFTIImageWriter()
strncpy(this->Description, "VTK", 3);
strncpy(&this->Description[3], version, l);
this->Description[l + 3] = '\0';
// Planar RGB (NIFTI doesn't allow this, it's here for Analyze)
this->PlanarRGB = false;
}
//----------------------------------------------------------------------------
......@@ -159,6 +161,7 @@ void vtkNIFTIImageWriter::PrintSelf(ostream& os, vtkIndent indent)
os << "(none)\n";
}
os << indent << "NIFTIVersion: " << this->NIFTIVersion << "\n";
os << indent << "PlanarRGB: " << (this->PlanarRGB ? "On\n" : "Off\n");
}
//----------------------------------------------------------------------------
......@@ -768,6 +771,11 @@ int vtkNIFTIImageWriter::RequestData(
unsigned char *dataPtr =
static_cast<unsigned char *>(data->GetScalarPointer());
// check if planar RGB is applicable (Analyze only)
bool planarRGB = (this->PlanarRGB &&
(this->OwnHeader->GetDataType() == NIFTI_TYPE_RGB24 ||
this->OwnHeader->GetDataType() == NIFTI_TYPE_RGBA32));
int swapBytes = 0;
int scalarSize = data->GetScalarSize();
int numComponents = data->GetNumberOfScalarComponents();
......@@ -781,10 +789,16 @@ int vtkNIFTIImageWriter::RequestData(
vectorDim *= timeDim;
z_off_t fileVoxelIncr = scalarSize*numComponents/vectorDim;
int planarSize = 1;
if (planarRGB)
{
planarSize = numComponents/vectorDim;
fileVoxelIncr = scalarSize;
}
// add a buffer for planar-vector to packed-vector conversion
unsigned char *rowBuffer = 0;
if (vectorDim > 1 || swapBytes)
if (vectorDim > 1 || planarRGB || swapBytes)
{
rowBuffer = new unsigned char[outSizeX*fileVoxelIncr];
}
......@@ -801,16 +815,29 @@ int vtkNIFTIImageWriter::RequestData(
dataPtr += sliceOffset*(outSizeZ - 1);
}
// special increment to handle planar RGB
vtkIdType planarOffset = 0;
vtkIdType planarEndOffset = 0;
if (planarRGB)
{
planarOffset = scalarSize*numComponents;
planarOffset *= outSizeX;
planarOffset *= outSizeY;
planarOffset -= scalarSize;
planarEndOffset = planarOffset - scalarSize*(planarSize - 1);
}
// report progress every 2% of the way to completion
vtkIdType target =
static_cast<vtkIdType>(0.02*outSizeY*outSizeZ*vectorDim) + 1;
static_cast<vtkIdType>(0.02*planarSize*outSizeY*outSizeZ*vectorDim) + 1;
vtkIdType count = 0;
// write the data one row at a time, do planar-to-packed conversion
// of vector components if NIFTI file has a vector dimension
int rowSize = numComponents/vectorDim*outSizeX;
int rowSize = fileVoxelIncr/scalarSize*outSizeX;
int c = 0; // counter for vector components
int j = 0; // counter for rows
int p = 0; // counter for planes (planar RGB)
int k = 0; // counter for slices
int t = 0; // counter for time
......@@ -818,7 +845,7 @@ int vtkNIFTIImageWriter::RequestData(
while (!this->AbortExecute && !this->ErrorCode)
{
if (vectorDim == 1 && swapBytes == 0)
if (vectorDim == 1 && !planarRGB && !swapBytes)
{
// write directly from input, instead of using a buffer
rowBuffer = ptr;
......@@ -867,28 +894,36 @@ int vtkNIFTIImageWriter::RequestData(
if (++j == outSizeY)
{
j = 0;
ptr -= 2*sliceOffset; // for reverse slice order
if (++k == outSizeZ)
// back up for next plane (R, G, or B) if planar mode
ptr -= planarOffset;
if (++p == planarSize)
{
k = 0;
if (++t == timeDim)
{
t = 0;
}
if (++c == vectorDim)
{
break;
}
// back up the ptr to the beginning of the image,
// then increment to the next vector component
ptr = dataPtr + c*fileVoxelIncr;
if (timeDim > 1)
p = 0;
ptr += planarEndOffset; // advance to start of next slice
ptr -= 2*sliceOffset; // for reverse slice order
if (++k == outSizeZ)
{
// if timeDim is included in the vectorDim (and hence in the
// VTK scalar components) then we have to make sure that
// the vector components are packed before the time steps
ptr = dataPtr + (c + t*(vectorDim - 1))/timeDim*fileVoxelIncr;
k = 0;
if (++t == timeDim)
{
t = 0;
}
if (++c == vectorDim)
{
break;
}
// back up the ptr to the beginning of the image,
// then increment to the next vector component
ptr = dataPtr + c*fileVoxelIncr*planarSize;
if (timeDim > 1)
{
// if timeDim is included in the vectorDim (and hence in the
// VTK scalar components) then we have to make sure that
// the vector components are packed before the time steps
ptr = dataPtr + (c + t*(vectorDim - 1))/timeDim*
fileVoxelIncr*planarSize;
}
}
}
}
......@@ -896,7 +931,7 @@ int vtkNIFTIImageWriter::RequestData(
// only delete this if it was alloced (if it was not alloced, it
// would have been set directly to a row out the output image)
if (vectorDim > 1 || swapBytes)
if (vectorDim > 1 || swapBytes || planarRGB)
{
delete [] rowBuffer;
}
......
......@@ -84,6 +84,16 @@ public:
vtkSetMacro(RescaleIntercept, double);
vtkGetMacro(RescaleIntercept, double);
// Description:
// Write planar RGB (separate R, G, and B planes), rather than packed RGB.
// Use this option with extreme caution: the NIFTI standard requires RGB
// pixels to be packed. The Analyze format, however, was used to store
// both planar RGB and packed RGB depending on the software, without any
// indication in the header about which convention was being used.
vtkGetMacro(PlanarRGB, bool);
vtkSetMacro(PlanarRGB, bool);
vtkBooleanMacro(PlanarRGB, bool);
// Description:
// The QFac sets the ordering of the slices in the NIFTI file.
// If QFac is -1, then the slice ordering in the file will be reversed
......@@ -169,6 +179,10 @@ protected:
vtkNIFTIImageHeader *OwnHeader;
int NIFTIVersion;
// Description:
// Use planar RGB instead of the default (packed).
bool PlanarRGB;
private:
vtkNIFTIImageWriter(const vtkNIFTIImageWriter&); // Not implemented.
void operator=(const vtkNIFTIImageWriter&); // Not implemented.
......
d4efcf01edb656d21c4d2b7f69fa1c80
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment