BUG: Fix loading of empty segmentation

NRRD file cannot be created without voxel data, so an 1x1x1 array was used as image data.
However, when the empty segmentation was loaded then this single-voxel geometry was treated as a valid geometry.
Fixed by detecting this special case of a single voxel in the segmentation, and considering the geometry invalid in this case.

fixes #7576

Docs/developer_guide/modules/segmentations.md

@@ -10,6 +10,9 @@ The **[slicerio](https://pypi.org/project/slicerio/) Python package** can read .
 
 If segments do not overlap then the file has 3 spatial dimensions. Even if a single-slice image is segmented, the spatial dimension must be still 3 because that is required for specification origin, spacing, and axis directions in 3D space. If segments overlap then the file has one `list` dimension and 3 spatial dimensions. Each 3D volume in the list is referred to as a `layer`.
 
+If the image contains a single voxel then it has a special meaning: it means that no image data is available and the image geometry (origin, spacing, axis directions, extents) will be ignored. This special case is necessary because it is not possible to create a nrrd image file without any voxel data.
+Therefore, a segmentation without segments or with only empty segments is saved to file (e.g., to be used as a template for new segmentations) then it will be saved with an image data containing a single voxel.
+
 ### Metadata
 
 Additional metadata is stored in custom data fields (starting with `Segmentation_` or `SegmentN_` prefixes), which provide hints on how the segments should be displayed or what they contain.

Libs/MRML/Core/Testing/CMakeLists.txt

@@ -193,6 +193,7 @@ simple_test( vtkMRMLSegmentationStorageNodeTest1
   DATA{${INPUT}/ITKSnapSegmentation.nii.gz}
   DATA{${INPUT}/OldSlicerSegmentation.seg.nrrd}
   DATA{${INPUT}/SlicerSegmentation.seg.nrrd}
+  ${TEMP}
   )
 simple_test( vtkMRMLSelectionNodeTest1 )
 simple_test( vtkMRMLSliceCompositeNodeTest1 )

Libs/MRML/Core/Testing/vtkMRMLSegmentationStorageNodeTest1.cxx

@@ -23,6 +23,7 @@ Care Ontario.
 #include "vtkMRMLScene.h"
 #include "vtkMRMLSegmentationNode.h"
 #include "vtkMRMLSegmentationStorageNode.h"
+#include "vtkOrientedImageData.h"
 #include "vtkSegmentationConverterFactory.h"
 
 // Converter rules
@@ -39,10 +40,10 @@ int vtkMRMLSegmentationStorageNodeTest1(int argc, char * argv[] )
   scene->AddNode(node1.GetPointer());
   EXERCISE_ALL_BASIC_MRML_METHODS(node1.GetPointer());
 
-  if (argc != 4)
+  if (argc != 5)
   {
     std::cerr << "Line " << __LINE__
-              << " - Missing parameters !\n"
+              << " - Missing or extra parameters!\n"
               << "Usage: " << argv[0] << " /path/to/ITKSnapSegmentation.nii.gz /path/to/OldSlicerSegmentation.seg.nrrd /path/to/SlicerSegmentation.seg.nrrd"
               << std::endl;
     return EXIT_FAILURE;
@@ -57,6 +58,7 @@ int vtkMRMLSegmentationStorageNodeTest1(int argc, char * argv[] )
   const char* itkSnapSegmentationFilename = argv[1]; // ITKSnapSegmentation.nii.gz
   const char* oldSlicerSegmentationFilename = argv[2]; // OldSlicerSegmentation.seg.nrrd: Segmentation before shared labelmaps implemented.
   const char* slicerSegmentationFilename = argv[3]; // SlicerSegmentation.seg.nrrd: Segmentation with shared labelmaps.
+  const char* tempDir = argv[4]; // Temporary folder where test segmentation files will be created
 
   // Test segmentation exported from ITK-SNAP
   std::cout << "Testing ITK-SNAP segmentation" << std::endl;
@@ -114,5 +116,44 @@ int vtkMRMLSegmentationStorageNodeTest1(int argc, char * argv[] )
     CHECK_INT(numberOfLayers, 2);
   }
 
-  return EXIT_SUCCESS;
+  std::cout << "Testing empty segmentation" << std::endl;
+  {
+    // Create empty segmentation
+    vtkNew<vtkMRMLSegmentationNode> segmentationNode;
+    scene->AddNode(segmentationNode);
+    segmentationNode->GetSegmentation()->AddEmptySegment();
+    segmentationNode->GetSegmentation()->AddEmptySegment();
+    segmentationNode->GetSegmentation()->AddEmptySegment();
+
+    // Write to file
+    vtkNew<vtkMRMLSegmentationStorageNode> segmentationStorageNode;
+    scene->AddNode(segmentationStorageNode);
+    std::string emptySegmentationFilename = std::string(tempDir) + "/EmptySegmentation.seg.nrrd";
+    std::cout << "Write empty segmentation file: " << emptySegmentationFilename;
+    segmentationStorageNode->SetFileName(emptySegmentationFilename.c_str());
+    CHECK_INT(segmentationStorageNode->WriteData(segmentationNode), 1);
+
+    // Read from file
+    vtkNew<vtkMRMLSegmentationNode> segmentationNodeFromFile;
+    scene->AddNode(segmentationNodeFromFile);
+    segmentationStorageNode->ReadData(segmentationNodeFromFile);
+
+    // Check basic content
+    vtkSegmentation* segmentation = segmentationNode->GetSegmentation();
+    CHECK_NOT_NULL(segmentation);
+    int numberOfSegments = segmentation->GetNumberOfSegments();
+    CHECK_INT(numberOfSegments, 3);
+
+    // Check that no valid geometry is found.
+    // The segmentation is stored as a single voxel, which would specify a geometry,
+    // the storage node should ignore that when reading the file (single-voxel volume is a special case).
+    std::string segmentationGeometryString = segmentation->DetermineCommonLabelmapGeometry(
+      vtkSegmentation::EXTENT_UNION_OF_EFFECTIVE_SEGMENTS_AND_REFERENCE_GEOMETRY);
+    CHECK_STD_STRING(segmentationGeometryString, "");
+
+    // Clean up
+    vtksys::SystemTools::RemoveFile(emptySegmentationFilename);
+  }
+
+    return EXIT_SUCCESS;
 }

Libs/MRML/Core/vtkMRMLSegmentationStorageNode.cxx

@@ -574,6 +574,7 @@ int vtkMRMLSegmentationStorageNode::ReadBinaryLabelmapRepresentation(vtkMRMLSegm
   vtkMatrix4x4* rasToFileIjk = nullptr;
   int imageExtentInFile[6] = { 0, -1, 0, -1, 0, -1 };
   int commonGeometryExtent[6] = { 0, -1, 0, -1, 0, -1 };
+  bool isExtentValid = false;
   int referenceImageExtentOffset[3] = { 0, 0, 0 };
 
   if (archetypeImageReader->CanReadFile(path.c_str()))
@@ -630,6 +631,29 @@ int vtkMRMLSegmentationStorageNode::ReadBinaryLabelmapRepresentation(vtkMRMLSegm
       imageData->GetExtent(commonGeometryExtent);
     }
 
+    // Special case: extent = [0, 0, 0, 0, 0, 0] means there is no image data
+    isExtentValid = true;
+    if (imageExtentInFile[0] == 0
+      && imageExtentInFile[1] == 0
+      && imageExtentInFile[2] == 0
+      && imageExtentInFile[3] == 0
+      && imageExtentInFile[4] == 0
+      && imageExtentInFile[5] == 0)
+    {
+      imageExtentInFile[1] = -1;
+      imageExtentInFile[3] = -1;
+      imageExtentInFile[5] = -1;
+
+      commonGeometryExtent[0] = 0;
+      commonGeometryExtent[1] = -1;
+      commonGeometryExtent[2] = 0;
+      commonGeometryExtent[3] = -1;
+      commonGeometryExtent[4] = 0;
+      commonGeometryExtent[5] = -1;
+
+      isExtentValid = false;
+    }
+
     // Read conversion parameters
     std::string conversionParameters;
     if (this->GetSegmentationMetaDataFromDicitionary(conversionParameters, dictionary, KEY_SEGMENTATION_CONVERSION_PARAMETERS))
@@ -713,7 +737,7 @@ int vtkMRMLSegmentationStorageNode::ReadBinaryLabelmapRepresentation(vtkMRMLSegm
     // Create binary labelmap volume
     vtkSmartPointer<vtkOrientedImageData> currentBinaryLabelmap = nullptr;
 
-    if (numberOfSegments == 0)
+    if (numberOfSegments == 0 && isExtentValid)
     {
       // No segment metadata. We are loading from a plain volume (not seg.nrrd).
 
@@ -1248,8 +1272,11 @@ int vtkMRMLSegmentationStorageNode::WriteBinaryLabelmapRepresentation(vtkMRMLSeg
     std::string commonGeometryString = segmentation->DetermineCommonLabelmapGeometry(
       this->CropToMinimumExtent ?
         vtkSegmentation::EXTENT_UNION_OF_EFFECTIVE_SEGMENTS : vtkSegmentation::EXTENT_UNION_OF_EFFECTIVE_SEGMENTS_AND_REFERENCE_GEOMETRY);
-    vtkSegmentationConverter::DeserializeImageGeometry(commonGeometryString, commonGeometryImage, true, scalarType, 1);
-    commonGeometryImage->GetExtent(commonGeometryExtent);
+    if (!commonGeometryString.empty())
+    {
+      vtkSegmentationConverter::DeserializeImageGeometry(commonGeometryString, commonGeometryImage, true, scalarType, 1);
+      commonGeometryImage->GetExtent(commonGeometryExtent);
+    }
   }
   if (commonGeometryExtent[0] > commonGeometryExtent[1]
     || commonGeometryExtent[2] > commonGeometryExtent[3]