Skip to content

GitLab

Commit 3522289d authored by Charles Law's avatar Charles Law
Browse files

Now supports images.

parent 776ea9f1
Hide whitespace changes
Inline Side-by-side
Parallel/vtkPDataSetReader.cxx
View file @ 3522289d
......@@ -42,6 +42,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include "vtkPDataSetReader.h"
#include "vtkDataSetReader.h"
#include "vtkStructuredPointsReader.h"
#include "vtkAppendPolyData.h"
#include "vtkAppendFilter.h"
#include "vtkPolyData.h"
......@@ -74,9 +75,11 @@ vtkPDataSetReader::vtkPDataSetReader()
{
this->FileName = NULL;
this->VTKFileFlag = 0;
this->StructuredFlag = 0;
this->NumberOfPieces = 0;
this->DataType = -1;
this->PieceFileNames = NULL;
this->PieceExtents = NULL;
}
//----------------------------------------------------------------------------
......@@ -90,10 +93,10 @@ vtkPDataSetReader::~vtkPDataSetReader()
for (int i = 0; i < this->NumberOfPieces; ++i)
{
if (this->PieceFileNames[i])
{
delete [] this->PieceFileNames[i];
this->PieceFileNames[i] = NULL;
}
{
delete [] this->PieceFileNames[i];
this->PieceFileNames[i] = NULL;
}
}
}
}
......@@ -109,7 +112,7 @@ void vtkPDataSetReader::SetNumberOfPieces(int num)
return;
}
// Delete the previous file names.
// Delete the previous file names/extents.
for (i = 0; i < this->NumberOfPieces; ++i)
{
if (this->PieceFileNames[i])
......@@ -117,12 +120,22 @@ void vtkPDataSetReader::SetNumberOfPieces(int num)
delete [] this->PieceFileNames[i];
this->PieceFileNames[i] = NULL;
}
if (this->PieceExtents && this->PieceExtents[i])
{
delete [] this->PieceExtents[i];
this->PieceExtents[i] = NULL;
}
}
if (this->PieceFileNames)
{
delete [] this->PieceFileNames;
this->PieceFileNames = NULL;
}
if (this->PieceExtents)
{
delete [] this->PieceExtents;
this->PieceExtents = NULL;
}
this->NumberOfPieces = 0;
......@@ -137,6 +150,16 @@ void vtkPDataSetReader::SetNumberOfPieces(int num)
{
this->PieceFileNames[i] = new char[512];
}
if (this->StructuredFlag)
{
this->PieceExtents = new int*[num];
for (i = 0; i < num; ++i)
{
this->PieceExtents[i] = new int[6];
}
}
this->NumberOfPieces = num;
}
......@@ -213,7 +236,7 @@ vtkDataSet *vtkPDataSetReader::CheckOutput()
newOutput = vtkImageData::New();
break;
case VTK_STRUCTURED_POINTS:
newOutput = vtkStructuredPoints::New();
newOutput = vtkImageData::New();
break;
default:
vtkErrorMacro("Unknown data type.");
......@@ -240,8 +263,6 @@ void vtkPDataSetReader::ExecuteInformation()
char str[1024];
vtkDataSet *output;
int i;
int dx, dy, dz;
float x, y, z;
file = this->OpenFile();
if (file == NULL)
......@@ -262,109 +283,15 @@ void vtkPDataSetReader::ExecuteInformation()
delete file;
return;
}
// Try to find the line that specifies the dataset type.
i = 0;
while (strncmp(str, "DATASET", 7) != 0 && i < 6)
{
file->getline(str, 1024);
++i;
}
if (strncmp(str, "DATASET POLYDATA", 16) == 0)
{
this->DataType = VTK_POLY_DATA;
}
else if (strncmp(str, "DATASET UNSTRUCTURED_GRID", 25) == 0)
{
this->DataType = VTK_UNSTRUCTURED_GRID;
}
else if (strncmp(str, "DATASET STRUCTURED_GRID", 23) == 0)
{
this->DataType = VTK_STRUCTURED_GRID;
vtkStructuredGrid *grid = (vtkStructuredGrid*)(this->CheckOutput());
file->getline(str, 1024, ' ');
if (strncmp(str, "DIMENSIONS", 10) != 0)
{
vtkErrorMacro("Expecting 'DIMENSIONS' insted of: " << str);
file->close();
delete file;
return;
}
*file >> dx;
*file >> dy;
*file >> dz;
grid->SetWholeExtent(0, dx-1, 0, dy-1, 0, dz-1);
}
else if (strncmp(str, "DATASET RECTILINEAR_GRID", 24) == 0)
{
this->DataType = VTK_RECTILINEAR_GRID;
vtkRectilinearGrid *grid = (vtkRectilinearGrid*)(this->CheckOutput());
file->getline(str, 1024, ' ');
if (strncmp(str, "DIMENSIONS", 10) != 0)
{
vtkErrorMacro("Expecting 'DIMENSIONS' insted of: " << str);
file->close();
delete file;
return;
}
*file >> dx;
*file >> dy;
*file >> dz;
grid->SetWholeExtent(0, dx-1, 0, dy-1, 0, dz-1);
}
else if (strncmp(str, "DATASET STRUCTURED_POINTS", 25) == 0)
{
this->DataType = VTK_STRUCTURED_POINTS;
vtkImageData *image = (vtkImageData*)(this->CheckOutput());
file->getline(str, 1024, ' ');
// hack to stop reading.
while (strlen(str) > 5)
{
if (strncmp(str, "DIMENSIONS", 10) == 0)
{
*file >> dx;
*file >> dy;
*file >> dz;
image->SetWholeExtent(0, dx-1, 0, dy-1, 0, dz-1);
}
if (strncmp(str, "SPACING", 7) == 0 || strncmp(str, "ASPECT_RATIO", 12) == 0)
{
*file >> x;
*file >> y;
*file >> z;
image->SetSpacing(x, y, z);
}
if (strncmp(str, "ORIGIN", 6) == 0)
{
*file >> x;
*file >> y;
*file >> z;
image->SetOrigin(x, y, z);
}
file->getline(str, 1024);
file->getline(str, 1024, ' ');
}
}
else
{
vtkErrorMacro("I can not figure out what type of data set this is: " << str);
delete file;
return;
}
// This is a vtk file not a PVTK file.
this->ReadVTKFileInformation(str, file);
this->VTKFileFlag = 1;
file->close();
delete file;
output = this->CheckOutput();
if (output->IsA("vtkPolyData") || output->IsA("vtkUnstructuredGrid"))
{
output->SetMaximumNumberOfPieces(1);
}
this->VTKFileFlag = 1;
return;
}
// Read PVTK meta data.
this->VTKFileFlag = 0;
file->getline(str, 1024);
if (strncmp(str, " dataType=", 15) != 0)
......@@ -377,22 +304,27 @@ void vtkPDataSetReader::ExecuteInformation()
if (strcmp(str+16, "vtkPolyData") == 0)
{
this->DataType = VTK_POLY_DATA;
this->StructuredFlag = 0;
}
else if (strcmp(str+10, "vtkUnstructuredGrid") == 0)
else if (strcmp(str+16, "vtkUnstructuredGrid") == 0)
{
this->DataType = VTK_UNSTRUCTURED_GRID;
this->StructuredFlag = 0;
}
else if (strcmp(str+10, "vtkstructuredGrid") == 0)
else if (strcmp(str+16, "vtkStructuredGrid") == 0)
{
this->DataType = VTK_STRUCTURED_GRID;
this->StructuredFlag = 1;
}
else if (strcmp(str+10, "vtkRectilinearGrid") == 0)
else if (strcmp(str+16, "vtkRectilinearGrid") == 0)
{
this->DataType = VTK_RECTILINEAR_GRID;
this->StructuredFlag = 1;
}
else if (strcmp(str+10, "vtkImageData") == 0)
else if (strcmp(str+16, "vtkImageData") == 0)
{
this->DataType = VTK_IMAGE_DATA;
this->StructuredFlag = 1;
}
else
{
......@@ -406,27 +338,237 @@ void vtkPDataSetReader::ExecuteInformation()
return;
}
// Read Image specific information
if (this->DataType == VTK_IMAGE_DATA)
{
this->ReadImageInformation(vtkImageData::SafeDownCast(output), str, file);
}
// Read whole extent for structured data.
if (this->StructuredFlag)
{
this->ReadWholeExtent(output, str, file);
}
// Read the number of pieces.
file->getline(str, 1024);
if (strncmp(str, " numberOfPieces=", 21) != 0)
{
vtkErrorMacro("Expecting numberofpieces.");
vtkErrorMacro("Expecting numberOfPieces.");
return;
}
// Take the last > off.
str[strlen(str)-2] = '\0';
this->SetNumberOfPieces(atoi(str+21));
str[strlen(str)-3] = '\0';
this->SetNumberOfPieces(atoi(str+22));
if (output->IsA("vtkPolyData") || output->IsA("vtkUnstructuredGrid"))
{
output->SetMaximumNumberOfPieces(this->NumberOfPieces);
}
// Read the filename and extents for each piece.
for (i = 0; i < this->NumberOfPieces; ++i)
{
file->getline(str, 512);
// Take all characters after the quote off.
str[strlen(str)-4] = '\0';
strcpy(this->PieceFileNames[i], str+19);
if (this->StructuredFlag)
{
int *pi = this->PieceExtents[i];
file->getline(str, 512);
// Take all characters after the quote off.
str[strlen(str)-1] = '\0';
strcpy(this->PieceFileNames[i], str+19);
// Now read the extent.
file->getline(str, 512);
// Take all characters after the quote off.
str[strlen(str)-4] = '\0';
sscanf(str+14, "%d %d %d %d %d %d", pi, pi+1, pi+2, pi+3, pi+4, pi+5);
}
else
{
file->getline(str, 512);
// Take all characters after the quote off.
str[strlen(str)-4] = '\0';
strcpy(this->PieceFileNames[i], str+19);
}
}
}
//----------------------------------------------------------------------------
// Does not read whole extent.
void vtkPDataSetReader::ReadImageInformation(vtkImageData *output,
char *str, ifstream *fp)
{
float origin[3];
float spacing[3];
if (output == NULL)
{
vtkErrorMacro("Expecting and image.");
return;
}
// scalarType="10"
// origin="-1.75 -1.25 0"
// spacing="0.1 0.1 0.1"
// Read the scalar type.
fp->getline(str, 1024);
if (strncmp(str, " scalarType=", 17) != 0)
{
vtkErrorMacro("Expecting scalarType.");
return;
}
// Take the quotes off.
str[strlen(str)-1] = '\0';
output->SetScalarType(atoi(str+18));
// Read the origin.
fp->getline(str, 1024);
if (strncmp(str, " origin=", 13) != 0)
{
vtkErrorMacro("Expecting origin.");
return;
}
// Take the quotes off.
str[strlen(str)-1] = '\0';
sscanf(str+14, "%f %f %f", origin, origin+1, origin+2);
output->SetOrigin(origin);
// Read the spacing.
fp->getline(str, 1024);
if (strncmp(str, " spacing=", 14) != 0)
{
vtkErrorMacro("Expecting spacing.");
return;
}
// Take the quotes off.
str[strlen(str)-1] = '\0';
sscanf(str+15, "%f %f %f", spacing, spacing+1, spacing+2);
output->SetSpacing(spacing);
}
//----------------------------------------------------------------------------
void vtkPDataSetReader::ReadWholeExtent(vtkDataSet *output,
char *str, ifstream *fp)
{
int ext[6];
if (output == NULL)
{
vtkErrorMacro("Expecting a data set.");
return;
}
// wholeExtent="0 25 0 25 0 10"
// Read the wholeExtent.
fp->getline(str, 1024);
if (strncmp(str, " wholeExtent=", 18) != 0)
{
vtkErrorMacro("Expecting wholeExtent.");
return;
}
// Take the quotes off.
str[strlen(str)-1] = '\0';
sscanf(str+19, "%d %d %d %d %d %d", ext, ext+1, ext+2, ext+3, ext+4, ext+5);
output->SetWholeExtent(ext);
}
//----------------------------------------------------------------------------
void vtkPDataSetReader::ReadVTKFileInformation(char *str, ifstream *file)
{
int i;
int dx, dy, dz;
float x, y, z;
vtkDataSet *output;
// Try to find the line that specifies the dataset type.
i = 0;
while (strncmp(str, "DATASET", 7) != 0 && i < 6)
{
file->getline(str, 1024);
++i;
}
if (strncmp(str, "DATASET POLYDATA", 16) == 0)
{
this->DataType = VTK_POLY_DATA;
}
else if (strncmp(str, "DATASET UNSTRUCTURED_GRID", 25) == 0)
{
this->DataType = VTK_UNSTRUCTURED_GRID;
}
else if (strncmp(str, "DATASET STRUCTURED_GRID", 23) == 0)
{
this->DataType = VTK_STRUCTURED_GRID;
vtkStructuredGrid *grid = (vtkStructuredGrid*)(this->CheckOutput());
file->getline(str, 1024, ' ');
if (strncmp(str, "DIMENSIONS", 10) != 0)
{
vtkErrorMacro("Expecting 'DIMENSIONS' insted of: " << str);
return;
}
*file >> dx;
*file >> dy;
*file >> dz;
grid->SetWholeExtent(0, dx-1, 0, dy-1, 0, dz-1);
}
else if (strncmp(str, "DATASET RECTILINEAR_GRID", 24) == 0)
{
this->DataType = VTK_RECTILINEAR_GRID;
vtkRectilinearGrid *grid = (vtkRectilinearGrid*)(this->CheckOutput());
file->getline(str, 1024, ' ');
if (strncmp(str, "DIMENSIONS", 10) != 0)
{
vtkErrorMacro("Expecting 'DIMENSIONS' insted of: " << str);
return;
}
*file >> dx;
*file >> dy;
*file >> dz;
grid->SetWholeExtent(0, dx-1, 0, dy-1, 0, dz-1);
}
else if (strncmp(str, "DATASET STRUCTURED_POINTS", 25) == 0)
{
this->DataType = VTK_IMAGE_DATA;
vtkImageData *image = (vtkImageData*)(this->CheckOutput());
file->getline(str, 1024, ' ');
// hack to stop reading.
while (strlen(str) > 5)
{
if (strncmp(str, "DIMENSIONS", 10) == 0)
{
*file >> dx;
*file >> dy;
*file >> dz;
image->SetWholeExtent(0, dx-1, 0, dy-1, 0, dz-1);
}
if (strncmp(str, "SPACING", 7) == 0 || strncmp(str, "ASPECT_RATIO", 12) == 0)
{
*file >> x;
*file >> y;
*file >> z;
image->SetSpacing(x, y, z);
}
if (strncmp(str, "ORIGIN", 6) == 0)
{
*file >> x;
*file >> y;
*file >> z;
image->SetOrigin(x, y, z);
}
file->getline(str, 1024);
file->getline(str, 1024, ' ');
}
}
else
{
vtkErrorMacro("I can not figure out what type of data set this is: " << str);
return;
}
output = this->CheckOutput();
if (output->IsA("vtkPolyData") || output->IsA("vtkUnstructuredGrid"))
{
output->SetMaximumNumberOfPieces(1);
}
}
......@@ -463,7 +605,9 @@ void vtkPDataSetReader::Execute()
reader->SetFileName(this->FileName);
reader->Update();
vtkDataSet *data = reader->GetOutput();
this->DataType = data->GetDataObjectType();
// Structured points giving me a pain.
//this->DataType = data->GetDataObjectType();
output = this->CheckOutput();
if (output == NULL)
{
......@@ -486,8 +630,11 @@ void vtkPDataSetReader::Execute()
case VTK_UNSTRUCTURED_GRID:
this->UnstructuredGridExecute();
break;
case VTK_IMAGE_DATA:
this->ImageDataExecute();
break;
default:
vtkErrorMacro("We do not handle structured data yet.");
vtkErrorMacro("We do not handle vtkStructuredGrid yet.");
return;
}
}
......@@ -618,6 +765,176 @@ void vtkPDataSetReader::UnstructuredGridExecute()
}
//----------------------------------------------------------------------------
// Structured data is trickier. Which files to load?
void vtkPDataSetReader::ImageDataExecute()
{
vtkStructuredPointsReader *reader;
vtkImageData *output;
int uExt[6];
int ext[6];
int *pieceMask;
int i, j;
// Use out internal method to get the output because GetOutput calls
// UpdateInformation.
output = vtkImageData::SafeDownCast(this->CheckOutput());
if (output == NULL)
{
vtkErrorMacro("Could not create output.");
return;
}
// Allocate the data object.
output->GetUpdateExtent(uExt);
output->SetExtent(uExt);
output->AllocateScalars();
// Get the pieces that will be read.
pieceMask = new int[this->NumberOfPieces];
for (i = 0; i < this->NumberOfPieces; ++i)
{
pieceMask[i] = 0;
}
this->CoverExtent(uExt, pieceMask);
// Now read and append
reader = vtkStructuredPointsReader::New();
for (i = 0; i < this->NumberOfPieces; ++i)
{
if (pieceMask[i])
{
reader->SetFileName(this->PieceFileNames[i]);
reader->Update();
// Sanity check: extent is correct. Ignore electric slide.
reader->GetOutput()->GetExtent(ext);
if (ext[1] - ext[0] != this->PieceExtents[i][1] - this->PieceExtents[i][0] ||
ext[3] - ext[2] != this->PieceExtents[i][3] - this->PieceExtents[i][2] ||
ext[5] - ext[4] != this->PieceExtents[i][5] - this->PieceExtents[i][4])
{
vtkErrorMacro("Unexpected extent in VTK file: " << this->PieceFileNames[i]);
}
else
{
// Reverse the electric slide.
reader->GetOutput()->SetExtent(this->PieceExtents[i]);
// Intersect extent and output extent
reader->GetOutput()->GetExtent(ext);
for (j = 0; j < 3; ++j)
{
if (ext[j*2] < uExt[j*2])
{
ext[j*2] = uExt[j*2];
}
if (ext[j*2+1] > uExt[j*2+1])
{
ext[j*2+1] = uExt[j*2+1];
}
}
output->CopyAndCastFrom(reader->GetOutput(), ext);
}
}
}
delete [] pieceMask;
reader->Delete();
}
//----------------------------------------------------------------------------
void vtkPDataSetReader::CoverExtent(int ext[6], int *pieceMask)
{
int bestArea;
int area;
int best;
int cExt[6]; // Covered
int rExt[6]; // Remainder piece
int i, j;
// Pick the piece with the largest coverage.
// Greedy search should be good enough.
best = -1;
bestArea = 0;
for (i = 0; i < this->NumberOfPieces; ++i)
{
// Compute coverage.
area = 1;
for (j = 0; j < 3; ++