Commit b78d13d1 authored by Dan Lipsa's avatar Dan Lipsa

Add an option to create a StructuredGrid for a 3D SegY dataset.

parent 08c1967e
......@@ -68,9 +68,6 @@ int TestSegY3DReader(int argc, char* argv[])
ren->AddActor(actor);
ren->ResetCamera();
//ren->GetActiveCamera()->Elevation(90);
//ren->GetActiveCamera()->Zoom(5);
// interact with data
renWin->Render();
......
......@@ -43,6 +43,7 @@ vtkSegYReader::vtkSegYReader()
this->DataExtent[1] = this->DataExtent[3] = this->DataExtent[5] = 0;
this->XYCoordMode = VTK_SEGY_SOURCE;
this->StructuredGrid = 0;
this->XCoordByte = 73;
this->YCoordByte = 77;
......@@ -125,15 +126,15 @@ int vtkSegYReader::RequestData(vtkInformation* vtkNotUsed(request),
}
}
this->Reader->LoadTraces();
if (this->Is3D)
if (this->Is3D && ! this->StructuredGrid)
{
this->Reader->ExportData3D(vtkImageData::SafeDownCast(output),
this->DataExtent, this->DataOrigin, this->DataSpacing);
this->Reader->ExportData(vtkImageData::SafeDownCast(output),
this->DataExtent, this->DataOrigin, this->DataSpacing);
}
else
{
vtkStructuredGrid* grid = vtkStructuredGrid::SafeDownCast(output);
this->Reader->ExportData2D(grid);
this->Reader->ExportData(grid, this->DataExtent);
grid->Squeeze();
}
this->Reader->In.close();
......@@ -158,7 +159,7 @@ int vtkSegYReader::RequestInformation(vtkInformation * vtkNotUsed(request),
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
this->DataExtent, 6);
std::cout << std::endl;
if (this->Is3D)
if (this->Is3D && ! this->StructuredGrid)
{
std::cout << "DataOrigin: ";
std::copy(this->DataOrigin, this->DataOrigin + 3, std::ostream_iterator<double>(std::cout, " "));
......@@ -175,7 +176,8 @@ int vtkSegYReader::RequestInformation(vtkInformation * vtkNotUsed(request),
int vtkSegYReader::FillOutputPortInformation(
int vtkNotUsed(port), vtkInformation* info)
{
const char* outputTypeName = this->Is3D ? "vtkImageData" : "vtkStructuredGrid";
const char* outputTypeName =
(this->Is3D && ! this->StructuredGrid) ? "vtkImageData" : "vtkStructuredGrid";
info->Set(vtkDataObject::DATA_TYPE_NAME(), outputTypeName);
std::cout << "FillOutputPortInformation" << std::endl;
return 1;
......@@ -210,7 +212,7 @@ int vtkSegYReader::RequestDataObject(
if (!output || !output->IsA(outputTypeName))
{
vtkDataSet* newOutput = nullptr;
if (this->Is3D)
if (this->Is3D && ! this->StructuredGrid)
{
newOutput = vtkImageData::New();
}
......
......@@ -25,9 +25,13 @@ class vtkImageData;
class vtkSegYReaderInternal;
/**
* Reader for SegY data. We create a vtkStructuredGrid for a 2.5D data
* and an vtkImageData for a 3D data. The structure of the data is created
* in the following order: crossline number, inline number, depth
* @class vtkSegYReader
* @brief Reads SegY data files.
*
* vtkSegYReader read the SegY data format. We create a
* vtkStructuredGrid for a 2.5D data and an vtkImageData for a 3D
* data. The structure of the data is created in the following order:
* crossline number, inline number, depth
*/
class VTKIOSEGY_EXPORT vtkSegYReader : public vtkDataSetAlgorithm
{
......@@ -60,7 +64,6 @@ public:
*/
vtkSetClampMacro(XYCoordMode, int, VTK_SEGY_SOURCE, VTK_SEGY_CUSTOM);
vtkGetMacro(XYCoordMode, int);
vtkBooleanMacro(XYCoordMode, int);
void SetXYCoordModeToSource();
void SetXYCoordModeToCDP();
void SetXYCoordModeToCustom();
......@@ -100,6 +103,19 @@ public:
vtkGetMacro(VerticalCRS, int);
//@}
//@{
/**
* Specify if we create a vtkStructuredGrid even if the data is
* 3D. Note this consumes more memory but it shows the precise
* location of the data. It may be useful if we want to show several
* datasets in the same window. The default value is false
* which means that we create a vtkImageData for a SegY 3D dataset.
*/
vtkSetMacro(StructuredGrid, int);
vtkGetMacro(StructuredGrid, int);
vtkBooleanMacro(StructuredGrid, int);
//@}
protected:
int RequestData(vtkInformation* request,
vtkInformationVector** inputVector,
......@@ -122,6 +138,7 @@ protected:
int DataExtent[6];
int XYCoordMode;
int StructuredGrid;
// Custom XY coordinate byte positions
int XCoordByte;
......
......@@ -24,7 +24,6 @@
#include "vtkSegYBinaryHeaderBytesPositions.h"
#include "vtkSegYIOUtils.h"
#include "vtkSegYTraceReader.h"
#include "vtkSmartPointer.h"
#include "vtkStructuredGrid.h"
#include <algorithm>
......@@ -219,8 +218,8 @@ bool vtkSegYReaderInternal::Is3DComputeParameters(
//-----------------------------------------------------------------------------
void vtkSegYReaderInternal::ExportData3D(vtkImageData* imageData,
int* extent, double* origin, double* spacing)
void vtkSegYReaderInternal::ExportData(vtkImageData* imageData,
int* extent, double* origin, double* spacing)
{
imageData->SetExtent(extent);
imageData->SetOrigin(origin);
......@@ -232,57 +231,59 @@ void vtkSegYReaderInternal::ExportData3D(vtkImageData* imageData,
scalars->SetNumberOfTuples(dims[0] * dims[1] * dims[2]);
scalars->SetName("trace");
imageData->GetPointData()->SetScalars(scalars);
for (int k = 0; k < this->SampleCountPerTrace; ++k)
int id = 0;
for (int k = 0; k < dims[2]; ++k)
{
for (int j = 0; j < dims[1]; ++j)
{
for (int i = 0; i < dims[0]; ++i)
{
vtkSegYTrace* trace = this->Traces[j * dims[0] + i];
scalars->SetValue(k * dims[1] * dims[0] + j * dims[0] + i,
trace->Data[k]);
scalars->SetValue(id++, trace->Data[k]);
}
}
}
}
//-----------------------------------------------------------------------------
void vtkSegYReaderInternal::ExportData2D(vtkStructuredGrid* grid)
void vtkSegYReaderInternal::ExportData(vtkStructuredGrid* grid, int* extent)
{
if (!grid)
{
return;
}
int crosslineCount = this->Traces.size();
grid->SetDimensions(crosslineCount, 1, this->SampleCountPerTrace);
vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
grid->SetExtent(extent);
int* dims = grid->GetDimensions();
vtkNew<vtkPoints> points;
vtkSmartPointer<vtkFloatArray> outScalars =
vtkSmartPointer<vtkFloatArray>::New();
outScalars->SetName("trace");
outScalars->SetNumberOfComponents(1);
outScalars->Allocate(crosslineCount * this->SampleCountPerTrace);
vtkNew<vtkFloatArray> scalars;
scalars->SetName("trace");
scalars->SetNumberOfComponents(1);
scalars->Allocate(dims[0] * dims[1] * dims[2]);
int sign = this->VerticalCRS == 0 ? -1 : 1;
int j = 0;
for (int k = 0; k < this->SampleCountPerTrace; k++)
int id = 0;
for (int k = 0; k < dims[2]; ++k)
{
for (unsigned int i = 0; i < crosslineCount; i++)
for (int j = 0; j < dims[1]; ++j)
{
auto trace = this->Traces[i];
double coordinateMultiplier = decodeMultiplier(trace->CoordinateMultiplier);
float x = trace->XCoordinate * coordinateMultiplier;
float y = trace->YCoordinate * coordinateMultiplier;
for (int i = 0; i < dims[0]; ++i)
{
auto trace = this->Traces[j * dims[0] + i];
double coordinateMultiplier = decodeMultiplier(trace->CoordinateMultiplier);
float x = trace->XCoordinate * coordinateMultiplier;
float y = trace->YCoordinate * coordinateMultiplier;
// The samples are uniformly placed at sample interval depths
// Dividing by 1000.0 to convert from microseconds to milliseconds.
float z = sign * k * (trace->SampleInterval / 1000.0);
points->InsertNextPoint(x, y, z);
// The samples are uniformly placed at sample interval depths
// Dividing by 1000.0 to convert from microseconds to milliseconds.
float z = sign * k * (trace->SampleInterval / 1000.0);
points->InsertNextPoint(x, y, z);
outScalars->InsertValue(j++, trace->Data[k]);
scalars->InsertValue(id++, trace->Data[k]);
}
}
}
grid->SetPoints(points);
grid->GetPointData()->SetScalars(outScalars);
grid->GetPointData()->SetScalars(scalars);
}
......@@ -40,8 +40,8 @@ public:
bool Is3DComputeParameters(int* extent, double* origin, double* spacing);
void LoadTraces();
void ExportData3D(vtkImageData*, int* extent, double* origin, double* spacing);
void ExportData2D(vtkStructuredGrid*);
void ExportData(vtkImageData*, int* extent, double* origin, double* spacing);
void ExportData(vtkStructuredGrid*, int* extent);
void SetXYCoordBytePositions(int x, int y);
void SetVerticalCRS(int);
......
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