GeneralTransform that has a GeoProjection does not correctly perform inverse

Development
1

Heyy,
when trying to take the inverse of a GeneralTransform that contains a GeoProjection it seems like the GeoProjection is “skipped” in the inverse transform. I believe this is a bug. Below is my code to reproduce. I think the comments quite clearly explain the problem. I tried posting this issue on GitLab but it told me “The form contains the following error: Your issue has been recognized as spam. Please, change the content to proceed.”

// Triangle.cxx (Apologies for the illogical filenames I edited a VTK template for this)
#include <vtkNew.h>
#include <vtkGeoTransform.h>
#include <vtkGeoProjection.h>
#include <vtkGeneralTransform.h>
#include <vtkTransform.h>

int main(int, char *[]) {
  // Create a mercator projection transform
  vtkNew<vtkGeoProjection> proj;
  proj->SetName("merc");
  vtkNew<vtkGeoTransform> geoTransform;
  geoTransform->SetDestinationProjection(proj);

  // Create simplest linear transform
  vtkNew<vtkTransform> linearTransform;
  linearTransform->Identity();

  // Compose them with totalProjection = linearTransform * geoTransform
  // i.e totalProjection(x) = linearTransform(geoTransform(x))
  vtkNew<vtkGeneralTransform> totalProjection;
  totalProjection->PostMultiply();
  totalProjection->Identity();
  totalProjection->Concatenate(geoTransform);
  totalProjection->Concatenate(linearTransform);

  // Choose random point transform and then take the inverse SHOULD get back original point.
  double in[3] = {4.846871030623073, 52.364810061968335, 0};
  totalProjection->TransformPoint(in, in);

  // Prints "in[3] = {539557,6.8322e+06,0}" - no problem here
  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;

  totalProjection->Inverse();
  totalProjection->TransformPoint(in, in);

  // Prints "in[3] = {539557,6.8322e+06,0}" - the same value as above?!
  // EXPECTED: "in[3] = {4.846871030623073, 52.364810061968335,0}" (possibly with some different rounding)
  cout << "in[3] = {" << in[0] << "," << in[1] << "," << in[2] << "}" << endl;

  return EXIT_SUCCESS;
}

# CMakeLists.txt
cmake_minimum_required(VERSION 3.12 FATAL_ERROR)

project(Triangle)

find_package(VTK COMPONENTS
  GeovisCore
)

if (NOT VTK_FOUND)
  message(FATAL_ERROR "Triangle: Unable to find the VTK build folder.")
endif()

# Prevent a "command line is too long" failure in Windows.
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
add_executable(Triangle MACOSX_BUNDLE Triangle.cxx )
  target_link_libraries(Triangle PRIVATE ${VTK_LIBRARIES}
)
# vtk_module_autoinit is needed
vtk_module_autoinit(
  TARGETS Triangle
  MODULES ${VTK_LIBRARIES}
)
2

It looks like the vtkGeoTransform is incomplete, it doesn’t implement the InternalDeepCopy() copy method. The vtkGeneralTransform needs to be able to copy the internal transforms before it inverts them.

To work around this bug, the easiest option is to create two copies of the GeoTransform:

  vtkNew<vtkGeoTransform> geoTransform;
  geoTransform->SetDestinationProjection(proj);

  vtkNew<vtkGeoTransform> geoTransformInverse;
  geoTransformInverse->SetSourceProjection(proj);

and two copies of the combined transform:

  vtkNew<vtkGeneralTransform> totalProjection;
  totalProjection->PostMultiply();
  totalProjection->Identity();
  totalProjection->Concatenate(geoTransform);
  totalProjection->Concatenate(linearTransform);

  vtkNew<vtkGeneralTransform> totalProjectionInverse;
  totalProjectionInverse->PostMultiply();
  totalProjectionInverse->Identity();
  totalProjectionInverse->Concatenate(linearTransform->GetInverse());
  totalProjectionInverse->Concatenate(geoTransformInverse);

A proper bug fix would require adding an InternalDeepCopy() method to vtkGeoTransform within the VTK code itself. Something like this:

void vtkGeoTransform::InternalDeepCopy(vtkAbstractTransform* transform)
{
  vtkGeoTransform* t = static_cast<vtkGeoTransform*>(transform);
  this->SetSourceProjection(t->GetSourceProjection());
  this->SetDestinationProjection(t->GetDestinationProjection());
  this->SetTransformZCoordinate(t->GetTransformZCoordinate());
}