Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Transform VTK files in accordance with their format and floating point precision #1253

Merged
merged 2 commits into from
Oct 23, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 78 additions & 46 deletions Core/ComponentBaseClasses/elxTransformBase.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -916,63 +916,95 @@ template <typename TElastix>
void
TransformBase<TElastix>::TransformPointsSomePointsVTK(const std::string & filename) const
{
/** Typedef's. \todo test DummyIPPPixelType=bool. */
using DummyIPPPixelType = float;
using MeshTraitsType =
itk::DefaultStaticMeshTraits<DummyIPPPixelType, FixedImageDimension, FixedImageDimension, CoordRepType>;
using MeshType = itk::Mesh<DummyIPPPixelType, FixedImageDimension, MeshTraitsType>;

/** Read the input points. */
const auto meshReader = itk::MeshFileReader<MeshType>::New();
meshReader->SetFileName(filename);
log::info(std::ostringstream{} << " Reading input point file: " << filename);
try
{
meshReader->Update();
}
catch (const itk::ExceptionObject & err)
{
log::error(std::ostringstream{} << " Error while opening input point file.\n" << err);
}

const auto & inputMesh = *(meshReader->GetOutput());

/** Some user-feedback. */
log::info(" Input points are specified in world coordinates.");
const unsigned long nrofpoints = inputMesh.GetNumberOfPoints();
log::info(std::ostringstream{} << " Number of specified input points: " << nrofpoints);

/** Apply the transform. */
log::info(" The input points are transformed.");
const itk::CommonEnums::IOComponent pointComponentType = [&filename] {
const itk::SmartPointer<itk::MeshIOBase> meshIO =
itk::MeshIOFactory::CreateMeshIO(filename.c_str(), itk::IOFileModeEnum::ReadMode);
meshIO->SetFileName(filename);
meshIO->ReadMeshInformation();
return meshIO->GetPointComponentType();
}();

// Reads a mesh from a VTK file, transforms the mesh, and writes the transformed mesh.
const auto ReadAndTransformAndWriteMesh = [&filename, this](auto coordinate) {
// The `coordinate` parameter just specifies the requested coordinate type.
(void)coordinate;

using DummyIPPPixelType = unsigned char;
using MeshType = itk::Mesh<
DummyIPPPixelType,
FixedImageDimension,
itk::DefaultStaticMeshTraits<DummyIPPPixelType, FixedImageDimension, FixedImageDimension, decltype(coordinate)>>;

/** Read the input points. */
const auto meshReader = itk::MeshFileReader<MeshType>::New();
meshReader->SetFileName(filename);
log::info(std::ostringstream{} << " Reading input point file: " << filename);
try
{
meshReader->Update();
}
catch (const itk::ExceptionObject & err)
{
log::error(std::ostringstream{} << " Error while opening input point file.\n" << err);
}

typename MeshType::ConstPointer outputMesh;
const auto & inputMesh = *(meshReader->GetOutput());

try
{
outputMesh = Self::TransformMesh(inputMesh);
}
catch (const itk::ExceptionObject & err)
{
log::error(std::ostringstream{} << " Error while transforming points.\n" << err);
}
/** Some user-feedback. */
log::info(" Input points are specified in world coordinates.");
const unsigned long nrofpoints = inputMesh.GetNumberOfPoints();
log::info(std::ostringstream{} << " Number of specified input points: " << nrofpoints);

const Configuration & configuration = itk::Deref(Superclass::GetConfiguration());
/** Apply the transform. */
log::info(" The input points are transformed.");

if (const std::string outputDirectoryPath = configuration.GetCommandLineArgument("-out");
!outputDirectoryPath.empty())
{
/** Create filename and file stream. */
const std::string outputPointsFileName = configuration.GetCommandLineArgument("-out") + "outputpoints.vtk";
log::info(std::ostringstream{} << " The transformed points are saved in: " << outputPointsFileName);
typename MeshType::ConstPointer outputMesh;

try
{
itk::WriteMesh(outputMesh, outputPointsFileName);
outputMesh = Self::TransformMesh(inputMesh);
}
catch (const itk::ExceptionObject & err)
{
log::error(std::ostringstream{} << " Error while saving points.\n" << err);
log::error(std::ostringstream{} << " Error while transforming points.\n" << err);
}

const Configuration & configuration = itk::Deref(Superclass::GetConfiguration());

if (const std::string outputDirectoryPath = configuration.GetCommandLineArgument("-out");
!outputDirectoryPath.empty())
{
/** Create filename and file stream. */
const std::string outputPointsFileName = configuration.GetCommandLineArgument("-out") + "outputpoints.vtk";
log::info(std::ostringstream{} << " The transformed points are saved in: " << outputPointsFileName);

try
{
const auto writer = itk::MeshFileWriter<MeshType>::New();
writer->SetInput(outputMesh);
writer->SetFileName(outputPointsFileName);

if (itk::Deref(meshReader->GetModifiableMeshIO()).GetFileType() == itk::IOFileEnum::Binary)
{
writer->SetFileTypeAsBINARY();
}

writer->Update();
}
catch (const itk::ExceptionObject & err)
{
log::error(std::ostringstream{} << " Error while saving points.\n" << err);
}
}
};

if (pointComponentType == itk::CommonEnums::IOComponent::FLOAT)
{
ReadAndTransformAndWriteMesh(float());
}
else
{
ReadAndTransformAndWriteMesh(double());
}

} // end TransformPointsSomePointsVTK()
Expand Down
44 changes: 44 additions & 0 deletions Testing/PythonTests/transformix_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -541,6 +541,50 @@ def test_zero_translation_of_vtk_2d_points(self) -> None:
output_mesh = reader.GetOutput()
self.assert_equal_mesh(output_mesh, input_mesh)

def test_zero_translation_of_vtk_2d_points_binary(self) -> None:
"""Tests zero-translation of VTK points in 2D"""

source_directory_path = pathlib.Path(__file__).resolve().parent
output_directory_path = self.create_test_function_output_directory()

parameter_directory_path = source_directory_path / "TransformParameters"

input_mesh = itk.Mesh[itk.D, 2].New()
for i in range(4):
input_mesh.SetPoint(
i, (self.random_finite_float32(), self.random_finite_float32())
)

writer = itk.MeshFileWriter[itk.Mesh[itk.D, 2]].New()
writer.SetInput(input_mesh)
writer.SetFileTypeAsBINARY()
writer.SetFileName(str(output_directory_path / "inputpoints.vtk"))
writer.Write()

subprocess.run(
[
str(self.transformix_exe_file_path),
"-def",
str(output_directory_path / "inputpoints.vtk"),
"-tp",
str(parameter_directory_path / "Translation(0,0).txt"),
"-out",
str(output_directory_path),
],
capture_output=True,
check=True,
)

# Note that itk.meshread does not work, as the following produces a 3D mesh, instead of a 2D
# mesh.
#
# output_mesh = itk.meshread(str(output_directory_path / "outputpoints.vtk"))
reader = itk.MeshFileReader[itk.Mesh[itk.D, 2]].New()
reader.SetFileName(str(output_directory_path / "outputpoints.vtk"))
reader.Update()
output_mesh = reader.GetOutput()
self.assert_equal_mesh(output_mesh, input_mesh)

def test_translation_deformation_field(self) -> None:
"""Tests zero-translation of VTK points in 2D"""

Expand Down
Loading