From 7af1e17286ab676a9c57b254bc64dbcdcd9fbf49 Mon Sep 17 00:00:00 2001
From: Matt McCormick <matt@mmmccormick.com>
Date: Mon, 23 Sep 2024 22:58:53 -0400
Subject: [PATCH] BUG: Add VTKPolyDataMeshIO FIELD data support

Add reading of VTK legacy FIELD data entries in both POINT_DATA
and CELL_DATA sections. When FIELD data is present and the number
of tuples matches the point/cell count, the data is read as a
VariableLengthVector pixel type.

Only the first field data array per section is used; additional
arrays are silently skipped. Add a test with a gourd.vtk file
that exercises the FIELD data path.

Co-authored-by: Matt McCormick <matt@mmmccormick.com>
---
 .../IO/MeshVTK/src/itkVTKPolyDataMeshIO.cxx   | 88 +++++++++++++++++++
 Modules/IO/MeshVTK/test/CMakeLists.txt        |  8 ++
 Modules/IO/MeshVTK/test/Input/gourd.vtk.cid   |  1 +
 3 files changed, 97 insertions(+)
 create mode 100755 Modules/IO/MeshVTK/test/Input/gourd.vtk.cid

diff --git a/Modules/IO/MeshVTK/src/itkVTKPolyDataMeshIO.cxx b/Modules/IO/MeshVTK/src/itkVTKPolyDataMeshIO.cxx
index 7129af2b97b..89fffb95c1e 100644
--- a/Modules/IO/MeshVTK/src/itkVTKPolyDataMeshIO.cxx
+++ b/Modules/IO/MeshVTK/src/itkVTKPolyDataMeshIO.cxx
@@ -655,6 +655,50 @@ VTKPolyDataMeshIO::ReadMeshInformation()
         this->m_NumberOfPointPixelComponents = this->m_PointDimension * (this->m_PointDimension + 1) / 2;
         this->m_UpdatePointData = true;
       }
+      if (line.find("FIELD") != std::string::npos)
+      {
+        StringStreamType sss;
+        sss << line;
+
+        sss >> item; // should be "FIELD"
+        if (item == "FIELD")
+        {
+          sss >> item; // field data name (e.g. "FieldData")
+          sss >> item; // number of field data arrays
+          const int numArrays = std::stoi(item);
+          if (numArrays >= 1)
+          {
+            // Read the first field data array header.
+            // Additional arrays are silently skipped.
+            std::getline(inputFile, line, '\n');
+
+            StringStreamType fieldStream;
+            fieldStream << line;
+            fieldStream >> item; // array name
+
+            fieldStream >> item; // numComponents
+            const unsigned int numComponents = std::stoi(item);
+
+            fieldStream >> item; // numTuples
+            const unsigned int numTuples = std::stoi(item);
+
+            fieldStream >> item; // data type
+            const IOComponentEnum componentType = this->GetComponentTypeFromString(item);
+
+            if (this->m_NumberOfPoints == numTuples)
+            {
+              this->m_NumberOfPointPixelComponents = numComponents;
+              this->m_PointPixelType = IOPixelEnum::VARIABLELENGTHVECTOR;
+              this->m_PointPixelComponentType = componentType;
+              if (this->m_PointPixelComponentType == IOComponentEnum::UNKNOWNCOMPONENTTYPE)
+              {
+                itkExceptionMacro("Unknown point pixel component type");
+              }
+              this->m_UpdatePointData = true;
+            }
+          }
+        }
+      }
     }
     else if (line.find("CELL_DATA") != std::string::npos)
     {
@@ -788,6 +832,50 @@ VTKPolyDataMeshIO::ReadMeshInformation()
         this->m_NumberOfCellPixelComponents = this->m_PointDimension * (this->m_PointDimension + 1) / 2;
         this->m_UpdateCellData = true;
       }
+      if (line.find("FIELD") != std::string::npos)
+      {
+        StringStreamType sss;
+        sss << line;
+
+        sss >> item; // should be "FIELD"
+        if (item == "FIELD")
+        {
+          sss >> item; // field data name (e.g. "FieldData")
+          sss >> item; // number of field data arrays
+          const int numArrays = std::stoi(item);
+          if (numArrays >= 1)
+          {
+            // Read the first field data array header.
+            // Additional arrays are silently skipped.
+            std::getline(inputFile, line, '\n');
+
+            StringStreamType fieldStream;
+            fieldStream << line;
+            fieldStream >> item; // array name
+
+            fieldStream >> item; // numComponents
+            const unsigned int numComponents = std::stoi(item);
+
+            fieldStream >> item; // numTuples
+            const unsigned int numTuples = std::stoi(item);
+
+            fieldStream >> item; // data type
+            const IOComponentEnum componentType = this->GetComponentTypeFromString(item);
+
+            if (this->m_NumberOfCells == numTuples)
+            {
+              this->m_NumberOfCellPixelComponents = numComponents;
+              this->m_CellPixelType = IOPixelEnum::VARIABLELENGTHVECTOR;
+              this->m_CellPixelComponentType = componentType;
+              if (this->m_CellPixelComponentType == IOComponentEnum::UNKNOWNCOMPONENTTYPE)
+              {
+                itkExceptionMacro("Unknown cell pixel component type");
+              }
+              this->m_UpdateCellData = true;
+            }
+          }
+        }
+      }
     }
     else if (line.find("OFFSETS") != std::string::npos)
     {
diff --git a/Modules/IO/MeshVTK/test/CMakeLists.txt b/Modules/IO/MeshVTK/test/CMakeLists.txt
index 982ec69c556..413f0b239e2 100644
--- a/Modules/IO/MeshVTK/test/CMakeLists.txt
+++ b/Modules/IO/MeshVTK/test/CMakeLists.txt
@@ -176,6 +176,14 @@ itk_add_test(
   2
   1
   1)
+itk_add_test(
+  NAME
+  itkMeshFileReadWriteTestField
+  COMMAND
+  ITKIOMeshVTKTestDriver
+  itkMeshFileReadWriteTest
+  DATA{Input/gourd.vtk}
+  ${ITK_TEST_OUTPUT_DIR}/gourd.vtk)
 
 set(ITKIOMeshVTKGTests itkVTKPolyDataMeshIOGTest.cxx)
 
diff --git a/Modules/IO/MeshVTK/test/Input/gourd.vtk.cid b/Modules/IO/MeshVTK/test/Input/gourd.vtk.cid
new file mode 100755
index 00000000000..fb4ea539551
--- /dev/null
+++ b/Modules/IO/MeshVTK/test/Input/gourd.vtk.cid
@@ -0,0 +1 @@
+bafkreidv66i2s5mleluikmvcpq32yziqt7qo3tvc5imncrx4auoimwtdua