Commit 7e882559 authored by Dan Lipsa's avatar Dan Lipsa

Fix SegY3D reader.

parent b68e22d5
......@@ -5,5 +5,6 @@ ExternalData_Expand_Arguments(VTKData _
vtk_add_test_cxx(vtkIOSegYCxxTests tests
TestSegY2DReader.cxx
TestSegY2DReaderZoom.cxx
TestSegY3DReader.cxx
)
vtk_test_cxx_executable(vtkIOSegYCxxTests tests)
/*=========================================================================
Program: Visualization Toolkit
Module: TestSegY2DReader.cxx
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.
=========================================================================*/
// .NAME Test of vtkSegY2DReader
// .SECTION Description
//
#include "vtkSegY3DReader.h"
#include "vtkActor.h"
#include "vtkCamera.h"
#include "vtkColorTransferFunction.h"
#include "vtkDataSetMapper.h"
#include "vtkIdList.h"
#include "vtkNew.h"
#include "vtkRegressionTestImage.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderer.h"
#include "vtkStructuredGrid.h"
#include "vtkTestUtilities.h"
int TestSegY3DReader(int argc, char* argv[])
{
// Basic visualisation.
vtkNew<vtkRenderWindow> renWin;
renWin->SetSize(300, 300);
vtkNew<vtkRenderer> ren;
renWin->AddRenderer(ren);
vtkNew<vtkRenderWindowInteractor> iren;
iren->SetRenderWindow(renWin);
// Read file name.
char* fname;
fname =
vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/SegY/waha8.sgy");
vtkNew<vtkColorTransferFunction> lut;
lut->AddRGBPoint(-127, 0.23, 0.30, 0.75);
lut->AddRGBPoint(0.0, 0.86, 0.86, 0.86);
lut->AddRGBPoint(126, 0.70, 0.02, 0.15);
vtkNew<vtkSegY3DReader> reader;
vtkNew<vtkDataSetMapper> mapper;
vtkNew<vtkActor> actor;
reader->SetFileName(fname);
reader->Update();
delete[] fname;
mapper->SetInputConnection(reader->GetOutputPort());
mapper->SetLookupTable(lut);
mapper->SetColorModeToMapScalars();
actor->SetMapper(mapper);
ren->AddActor(actor);
ren->ResetCamera();
ren->GetActiveCamera()->Elevation(90);
ren->GetActiveCamera()->Zoom(5);
// interact with data
renWin->Render();
int retVal = vtkRegressionTestImage(renWin);
if (retVal == vtkRegressionTester::DO_INTERACTOR)
{
iren->Start();
}
return !retVal;
}
......@@ -9,7 +9,7 @@ vtk_module(vtkIOSegY
DEPENDS
vtkCommonDataModel
vtkCommonExecutionModel
vtkIOCore
vtkIOImage
PRIVATE_DEPENDS
vtkCommonCore
)
......@@ -13,45 +13,214 @@
=========================================================================*/
#include "vtkSegY3DReader.h"
#include "vtkImageData.h"
#include "vtkInformationVector.h"
#include "vtkInformation.h"
#include "vtkObjectFactory.h"
#include "vtkSegYReader.h"
#include "vtkSmartPointer.h"
#include "vtkSegY3DReader.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStructuredGrid.h"
#include <algorithm>
#include <iostream>
#include <iterator>
vtkStandardNewMacro(vtkSegY3DReader);
//-----------------------------------------------------------------------------
vtkSegY3DReader::vtkSegY3DReader()
{
this->image = nullptr;
this->SetNumberOfInputPorts(0);
this->Reader = new vtkSegYReader();
this->FileName = nullptr;
this->Is3D = false;
this->DataOrigin[0] = this->DataOrigin[1] = this->DataOrigin[2] = 0.0;
this->DataSpacing[0] = this->DataSpacing[1] = this->DataSpacing[2] = 1.0;
this->DataExtent[0] = this->DataExtent[2] = this->DataExtent[4] = 0;
this->DataExtent[1] = this->DataExtent[3] = this->DataExtent[5] = 0;
this->XYCoordMode = VTK_SEGY_SOURCE;
this->XCoordByte = 73;
this->YCoordByte = 77;
this->VerticalCRS = VTK_SEGY_VERTICAL_HEIGHTS;
}
//-----------------------------------------------------------------------------
vtkSegY3DReader::~vtkSegY3DReader()
{
if (this->image)
{
this->image = nullptr;
}
delete this->Reader;
}
//-----------------------------------------------------------------------------
void vtkSegY3DReader::SetXYCoordModeToSource()
{
this->SetXYCoordMode(VTK_SEGY_SOURCE);
}
//-----------------------------------------------------------------------------
void vtkSegY3DReader::SetXYCoordModeToCDP()
{
this->SetXYCoordMode(VTK_SEGY_CDP);
}
//-----------------------------------------------------------------------------
void vtkSegY3DReader::SetXYCoordModeToCustom()
{
this->SetXYCoordMode(VTK_SEGY_CUSTOM);
}
//-----------------------------------------------------------------------------
void vtkSegY3DReader::PrintSelf(ostream& os, vtkIndent indent)
{
os << indent << "FileName: " << (this->FileName ? this->FileName : "(none)")
<< "\n";
Superclass::PrintSelf(os, indent);
}
//-----------------------------------------------------------------------------
int vtkSegY3DReader::RequestData(vtkInformation* vtkNotUsed(request),
vtkInformationVector** vtkNotUsed(inputVector),
vtkInformationVector* outputVector)
{
vtkInformation* outInfo = outputVector->GetInformationObject(0);
if (!outInfo)
{
return 0;
}
vtkDataObject* output = outInfo->Get(vtkDataObject::DATA_OBJECT());
if (!output)
{
return 0;
}
if (this->Is3D)
{
this->Reader->LoadTraces();
this->Reader->ExportData3D(vtkImageData::SafeDownCast(output),
this->DataExtent, this->DataOrigin, this->DataSpacing);
}
else
{
switch (this->XYCoordMode)
{
case VTK_SEGY_SOURCE:
{
this->Reader->SetXYCoordBytePositions(72, 76);
break;
}
case VTK_SEGY_CDP:
{
this->Reader->SetXYCoordBytePositions(180, 184);
break;
}
case VTK_SEGY_CUSTOM:
{
this->Reader->SetXYCoordBytePositions(this->XCoordByte - 1,
this->YCoordByte - 1);
break;
}
default:
{
vtkErrorMacro(<< "Unknown value for XYCoordMode " << this->XYCoordMode);
return 1;
}
}
vtkStructuredGrid* grid = vtkStructuredGrid::SafeDownCast(output);
this->Reader->SetVerticalCRS(this->VerticalCRS);
this->Reader->LoadTraces();
this->Reader->ExportData2D(grid);
grid->Squeeze();
}
this->Reader->In.close();
std::cout << "RequestData" << std::endl;
return 1;
}
//-----------------------------------------------------------------------------
vtkImageData* vtkSegY3DReader::GetImage()
int vtkSegY3DReader::RequestInformation(vtkInformation * vtkNotUsed(request),
vtkInformationVector **vtkNotUsed(inputVector),
vtkInformationVector *outputVector)
{
this->reader.LoadFromFile(FileName);
if (this->Is3D)
{
vtkInformation* outInfo = outputVector->GetInformationObject(0);
if (!outInfo)
{
vtkErrorMacro("Invalid output information object");
return 0;
}
image = vtkSmartPointer<vtkImageData>::New();
if (!reader.ExportData3D(image))
cout << "Failed to export 3D image from reader" << endl;
std::cout << "DataExtent: ";
std::copy(this->DataExtent, this->DataExtent + 6, std::ostream_iterator<int>(std::cout, " "));
std::cout << "\nDataOrigin: ";
std::copy(this->DataOrigin, this->DataOrigin + 3, std::ostream_iterator<double>(std::cout, " "));
std::cout << "\nDataSpacing: ";
std::copy(this->DataSpacing, this->DataSpacing + 3, std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
this->DataExtent, 6);
outInfo->Set(vtkDataObject::ORIGIN(), this->DataOrigin, 3);
outInfo->Set(vtkDataObject::SPACING(), this->DataSpacing, 3);
}
return 1;
}
return image;
//----------------------------------------------------------------------------
int vtkSegY3DReader::FillOutputPortInformation(
int vtkNotUsed(port), vtkInformation* info)
{
const char* outputTypeName = this->Is3D ? "vtkImageData" : "vtkStructuredGrid";
info->Set(vtkDataObject::DATA_TYPE_NAME(), outputTypeName);
std::cout << "FillOutputPortInformation" << std::endl;
return 1;
}
//----------------------------------------------------------------------------
int vtkSegY3DReader::RequestDataObject(
vtkInformation*,
vtkInformationVector** vtkNotUsed(inputVector),
vtkInformationVector* outputVector)
{
vtkInformation* info = outputVector->GetInformationObject(0);
vtkDataSet *output = vtkDataSet::SafeDownCast(
info->Get(vtkDataObject::DATA_OBJECT()));
if (!this->FileName)
{
vtkErrorMacro("Requires valid input file name") ;
return 0;
}
this->Reader->In.open(this->FileName, std::ifstream::binary);
if (!this->Reader->In)
{
std::cerr << "File not found:" << this->FileName << std::endl;
return 0;
}
this->Is3D = this->Reader->Is3DComputeParameters(
this->DataExtent, this->DataOrigin, this->DataSpacing);
const char* outputTypeName = this->Is3D ? "vtkImageData" : "vtkStructuredGrid";
if (!output || !output->IsA(outputTypeName))
{
vtkDataSet* newOutput = nullptr;
if (this->Is3D)
{
newOutput = vtkImageData::New();
}
else
{
newOutput = vtkStructuredGrid::New();
}
info->Set(vtkDataObject::DATA_OBJECT(), newOutput);
newOutput->Delete();
}
std::cout << "RequestDataObject" << std::endl;
return 1;
}
......@@ -16,36 +16,114 @@
#ifndef vtkSegY3DReader_h
#define vtkSegY3DReader_h
#include "vtkImageAlgorithm.h"
#include "vtkSegYReader.h" // For reader implementation
#include "vtkSmartPointer.h" // For smart pointer
#include "vtkDataSetAlgorithm.h"
#include <vtkIOSegYModule.h> // For export macro
// Forward declarations
class vtkImageData;
class vtkSegYReader;
class VTKIOSEGY_EXPORT vtkSegY3DReader : public vtkImageAlgorithm
class VTKIOSEGY_EXPORT vtkSegY3DReader : public vtkDataSetAlgorithm
{
public:
static vtkSegY3DReader* New();
vtkTypeMacro(vtkSegY3DReader, vtkImageAlgorithm);
vtkTypeMacro(vtkSegY3DReader, vtkDataSetAlgorithm);
void PrintSelf(ostream& os, vtkIndent indent) override;
vtkSegY3DReader();
~vtkSegY3DReader();
vtkSetStringMacro(FileName);
vtkGetStringMacro(FileName);
virtual vtkImageData* GetImage();
enum VTKSegYCoordinateModes
{
VTK_SEGY_SOURCE = 0, // default
VTK_SEGY_CDP = 1,
VTK_SEGY_CUSTOM = 2
};
//@{
/**
* Specify whether to use source x/y coordinates or CDP coordinates or custom
* byte positions for data position in the SEG-Y trace header. Defaults to
* source x/y coordinates.
*
* As per SEG-Y rev 2.0 specification,
* Source XY coordinate bytes = (73, 77)
* CDP XY coordinate bytes = (181, 185)
*/
vtkSetClampMacro(XYCoordMode, int, VTK_SEGY_SOURCE, VTK_SEGY_CUSTOM);
vtkGetMacro(XYCoordMode, int);
vtkBooleanMacro(XYCoordMode, int);
void SetXYCoordModeToSource();
void SetXYCoordModeToCDP();
void SetXYCoordModeToCustom();
//@}
//@{
/**
* Specify X and Y byte positions for custom XYCoordinateMode.
* By default, XCoordByte = 73, YCoordByte = 77 i.e. source xy.
*
* \sa SetXYCoordinatesModeToCustom()
*/
vtkSetMacro(XCoordByte, int);
vtkGetMacro(XCoordByte, int);
vtkSetMacro(YCoordByte, int);
vtkGetMacro(YCoordByte, int);
//@}
enum VTKSegYVerticalCRS
{
VTK_SEGY_VERTICAL_HEIGHTS = 0, // default
VTK_SEGY_VERTICAL_DEPTHS
};
//@{
/**
* Specify whether the vertical coordinates in the SEG-Y file are heights
* (positive up) or depths (positive down). By default, the vertical
* coordinates are treated as heights (i.e. positive up). This means that the
* Z-axis of the dataset goes from 0 (surface) to -ve depth (last sample).
* \note As per the SEG-Y rev 2.0 specification, this information is defined
* in the Location Data Stanza of the Extended Textual Header. However, as of
* this revision, vtkSegY2DReader does not support reading the extended
* textual header.
*/
vtkSetMacro(VerticalCRS, int);
vtkGetMacro(VerticalCRS, int);
//@}
protected:
vtkSegY3DReader();
~vtkSegY3DReader();
int RequestData(vtkInformation* request,
vtkInformationVector** inputVector,
vtkInformationVector* outputVector) override;
int RequestInformation(vtkInformation* request,
vtkInformationVector** inputVector,
vtkInformationVector* outputVector) override;
int FillOutputPortInformation(int port, vtkInformation* info) override;
int RequestDataObject(vtkInformation* request,
vtkInformationVector** inputVector,
vtkInformationVector* outputVector) override;
protected:
vtkSegYReader* Reader;
char *FileName;
bool Is3D;
double DataOrigin[3];
double DataSpacing[3];
int DataExtent[6];
int XYCoordMode;
// Custom XY coordinate byte positions
int XCoordByte;
int YCoordByte;
int VerticalCRS;
char* FileName;
vtkSegYReader reader;
vtkSmartPointer<vtkImageData> image;
private:
vtkSegY3DReader(const vtkSegY3DReader&) = delete;
......
......@@ -62,5 +62,6 @@ private:
}
};
#endif // vtkSegYBinaryHeaderBytesPositions_h
// VTK-HeaderTest-Exclude: vtkSegYBinaryHeaderBytesPositions.h
......@@ -22,7 +22,7 @@
//----------------------------------------------------------------------------
vtkSegYIOUtils::vtkSegYIOUtils()
{
isBigEndian = checkIfBigEndian();
this->IsBigEndian = checkIfBigEndian();
}
//----------------------------------------------------------------------------
......@@ -39,7 +39,7 @@ int vtkSegYIOUtils::readShortInteger(int pos, std::ifstream& in)
char buffer[2];
in.read(buffer, sizeof(buffer));
if (!isBigEndian)
if (!this->IsBigEndian)
{
swap(buffer, buffer + 1);
}
......@@ -56,7 +56,7 @@ int vtkSegYIOUtils::readLongInteger(int pos, std::ifstream& in)
char buffer[4];
in.read(buffer, sizeof(buffer));
if (!isBigEndian)
if (!this->IsBigEndian)
{
swap(buffer, buffer + 3);
swap(buffer + 1, buffer + 2);
......@@ -73,7 +73,7 @@ int vtkSegYIOUtils::readLongInteger(std::ifstream& in)
char buffer[4];
in.read(buffer, sizeof(buffer));
if (!isBigEndian)
if (!this->IsBigEndian)
{
swap(buffer, buffer + 3);
swap(buffer + 1, buffer + 2);
......@@ -90,7 +90,7 @@ float vtkSegYIOUtils::readFloat(std::ifstream& in)
char buffer[4];
in.read(buffer, sizeof(buffer));
if (!isBigEndian)
if (!this->IsBigEndian)
{
swap(buffer, buffer + 3);
swap(buffer + 1, buffer + 2);
......@@ -107,7 +107,7 @@ float vtkSegYIOUtils::readIBMFloat(std::ifstream& in)
char buffer[4];
in.read(buffer, sizeof(buffer));
if (!isBigEndian)
if (!this->IsBigEndian)
{
swap(buffer, buffer + 3);
swap(buffer + 1, buffer + 2);
......
......@@ -22,18 +22,19 @@
class vtkSegYIOUtils
{
public:
char readChar(std::ifstream& in);
int readShortInteger(int pos, std::ifstream& in);
int readLongInteger(int pos, std::ifstream& in);
int readLongInteger(std::ifstream& in);
float readFloat(std::ifstream& in);
float readIBMFloat(std::ifstream& in);
char readChar(std::ifstream& in);
unsigned char readUChar(std::ifstream& in);
void swap(char* a, char* b);
bool isBigEndian;
static vtkSegYIOUtils* Instance();
int getFileSize(std::ifstream& in);
bool IsBigEndian;
private:
vtkSegYIOUtils();
bool checkIfBigEndian()
......
This diff is collapsed.
......@@ -35,22 +35,32 @@ public:
~vtkSegYReader();
public:
bool GetImageData(vtkImageData*); // export the data in segy file as 3D volume
bool ExportData3D(vtkImageData*); // export the data in segy file as 3D volume
bool LoadFromFile(std::string path);
bool Is3DComputeParameters(int* extent, double* origin, double* spacing);
void LoadTraces();
void ExportData3D(vtkImageData*, int* extent, double* origin, double* spacing);
void ExportData2D(vtkStructuredGrid*);
void AddScalars(vtkStructuredGrid*);
void SetXYCoordBytePositions(int x, int y);
void SetVerticalCRS(int);
std::ifstream In;
protected:
bool ReadHeader();
private:
bool ReadHeader(std::ifstream& in);
std::vector<vtkSegYTrace*> Traces;
int FormatCode;
vtkSegYBinaryHeaderBytesPositions* BinaryHeaderBytesPos;
vtkSegYTraceReader* TraceReader;
int SampleCountPerTrace;
int VerticalCRS;
// Binary Header
short SampleInterval;
int FormatCode;
int SampleCountPerTrace;
};
#endif
......
......@@ -90,23 +90,15 @@ void vtkSegYTraceReader::PrintTraceHeader(std::ifstream& in, int startPos)
}
//-----------------------------------------------------------------------------
bool vtkSegYTraceReader::ReadTrace(int& startPos,
std::ifstream& in,
int formatCode,
vtkSegYTrace* trace)
void vtkSegYTraceReader::ReadTrace(int& startPos,
std::ifstream& in,
int formatCode,
vtkSegYTrace* trace)
{
int fileSize = vtkSegYIOUtils::Instance()->getFileSize(in);
if (startPos + 240 >= fileSize)
{
return false;
}
// PrintTraceHeader(in, startPos);
trace->crosslineNumber = vtkSegYIOUtils::Instance()->readLongInteger(
startPos + traceHeaderBytesPos.CrosslineNumber, in);
trace->inlineNumber = vtkSegYIOUtils::Instance()->readLongInteger(
startPos + traceHeaderBytesPos.InlineNumber, in);
trace->crosslineNumber = vtkSegYIOUtils::Instance()->readLongInteger(
startPos + traceHeaderBytesPos.CrosslineNumber, in);
int numSamples = vtkSegYIOUtils::Instance()->readShortInteger(
startPos + traceHeaderBytesPos.NumberSamples, in);
trace->xCoordinate = vtkSegYIOUtils::Instance()->readLongInteger(
......@@ -119,29 +111,56 @@ bool vtkSegYTraceReader::ReadTrace(int& startPos,
startPos + traceHeaderBytesPos.SampleInterval, in);
in.seekg(startPos + 240, in.beg);
for (int i = 0; i < numSamples; i++)
float value;
switch (formatCode)
{
float value;
switch (formatCode)
case 1:
for (int i = 0; i < numSamples; i++)
{
value = vtkSegYIOUtils::Instance()->readIBMFloat(in);
trace->data.push_back(value);
}
break;
case 5:
for (int i = 0; i < numSamples; i++)
{
value = vtkSegYIOUtils::Instance()->readFloat(in);
trace->data.push_back(value);
}
break;
case 8:
for (int i = 0; i < numSamples; i++)
{
case 1:
value = vtkSegYIOUtils::Instance()->readIBMFloat(in);
break;
case 5:
value = vtkSegYIOUtils::Instance()->readFloat(in);
break;
default:
std::cerr << "Data sample format code " << formatCode
<< " not supported." << std::endl;
value = 0;
value = vtkSegYIOUtils::Instance()->readChar(in);
trace->data.push_back(value);
}
trace->data.push_back(value);
break;
default:
std::cerr << "Data sample format code " << formatCode
<< " not supported." << std::endl;
value = 0;
}
startPos += 240 + GetTraceSize(numSamples, formatCode);
return true;
startPos += 240 + this->GetTraceSize(numSamples, formatCode);
}
//-----------------------------------------------------------------------------
void vtkSegYTraceReader::ReadInlineCrossline(int& startPos,
std::ifstream& in,
int formatCode,
int* inlineNumber, int* crosslineNumber)
{
*inlineNumber = vtkSegYIOUtils::Instance()->readLongInteger(
startPos + traceHeaderBytesPos.InlineNumber, in);
*crosslineNumber = vtkSegYIOUtils::Instance()->readLongInteger(
startPos + traceHeaderBytesPos.CrosslineNumber, in);
int numSamples = vtkSegYIOUtils::Instance()->readShortInteger(
startPos + traceHeaderBytesPos.NumberSamples, in);
startPos += 240 + this->GetTraceSize(numSamples, formatCode);
}
//-----------------------------------------------------------------------------
int vtkSegYTraceReader::GetTraceSize(int numSamples, int formatCode)
{
......
......@@ -53,10 +53,15 @@ public:
void SetXYCoordBytePositions(int x, int y);
void PrintTraceHeader(std::ifstream& in, int startPos);
bool ReadTrace(int& startPos,
std::ifstream& in,
int formatCode,
vtkSegYTrace* trace);
void ReadTrace(int& startPos,
std::ifstream& in,
int formatCode,
vtkSegYTrace* trace);
void ReadInlineCrossline(int& startPos,
std::ifstream& in,
int formatCode,
int* inlineNumber, int* crosslineNumber);
int GetTraceSize(int numSamples, int formatCode);
};
......
fdbdde34486d2d89163c9d2d5e4ab340
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