fix(image-sets-normalization): fix sorting of files in series

Author: PaulHax

packages/dicom/gdcm/SortSpatially.h

@@ -0,0 +1,94 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         https://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef SORT_SPATIALLY_H
+#define SORT_SPATIALLY_H
+
+#include <string>
+#include <unordered_set>
+
+#include "gdcmSerieHelper.h"
+
+class CustomSerieHelper : public gdcm::SerieHelper
+{
+public:
+    void AddFileName(std::string const &fileName)
+    {
+        SerieHelper::AddFileName(fileName);
+    }
+};
+
+using FileNamesContainer = std::vector<std::string>;
+
+FileNamesContainer sortSpatially(std::vector<std::string> unsortedSerieFileNames)
+{
+    std::unique_ptr<CustomSerieHelper> serieHelper(new CustomSerieHelper());
+    for (const std::string &fileName : unsortedSerieFileNames)
+    {
+        serieHelper->AddFileName(fileName);
+    }
+    serieHelper->SetUseSeriesDetails(true);
+    // Add the default restrictions to refine the file set into multiple series.
+    serieHelper->CreateDefaultUniqueSeriesIdentifier();
+    using SeriesIdContainer = std::vector<std::string>;
+    SeriesIdContainer seriesUIDs;
+    // Accessing the first serie found (assume there is at least one)
+    gdcm::FileList *flist = serieHelper->GetFirstSingleSerieUIDFileSet();
+    while (flist)
+    {
+        if (!flist->empty()) // make sure we have at least one serie
+        {
+            gdcm::File *file = (*flist)[0]; // for example take the first one
+
+            // Create its unique series ID
+            const std::string id(serieHelper->CreateUniqueSeriesIdentifier(file));
+
+            seriesUIDs.push_back(id);
+        }
+        flist = serieHelper->GetNextSingleSerieUIDFileSet();
+    }
+
+    FileNamesContainer fileNames;
+    flist = serieHelper->GetFirstSingleSerieUIDFileSet();
+    const std::string serie = seriesUIDs[0];
+    bool found = false;
+    while (flist && !found)
+    {
+        if (!flist->empty()) // make sure we have at least one serie
+        {
+            gdcm::File *file = (*flist)[0]; // for example take the first one
+            const std::string id(serieHelper->CreateUniqueSeriesIdentifier(file));
+            if (id == serie)
+            {
+                found = true; // we found a match
+                break;
+            }
+        }
+        flist = serieHelper->GetNextSingleSerieUIDFileSet();
+    }
+    serieHelper->OrderFileList(flist);
+
+    gdcm::FileList::iterator it;
+    for (it = flist->begin(); it != flist->end(); ++it)
+    {
+        gdcm::FileWithName *header = *it;
+        fileNames.push_back(header->filename);
+    }
+    return fileNames;
+}
+
+#endif // SORT_SPATIALLY_H

packages/dicom/gdcm/image-sets-normalization.cxx

@@ -15,25 +15,25 @@
  *  limitations under the License.
  *
  *=========================================================================*/
+#include <unordered_set>
+
+#include "rapidjson/document.h"
+#include "rapidjson/stringbuffer.h"
+#include "rapidjson/writer.h"
+
 #include "itkPipeline.h"
 #include "itkOutputTextStream.h"
 
 #include "gdcmDiscriminateVolume.h"
 #include "DICOMTagReader.h"
 #include "Tags.h"
-
-#include "rapidjson/document.h"
-#include "rapidjson/stringbuffer.h"
-#include "rapidjson/writer.h"
-
-#include <unordered_set>
+#include "SortSpatially.h"
 
 const std::string STUDY_INSTANCE_UID = "0020|000D";
 const std::string SERIES_INSTANCE_UID = "0020|000e";
 
 using File = std::string;
 using TagValues = itk::DICOMTagReader::TagMapType;
-using FileToTags = std::unordered_map<File, TagValues>;
 
 rapidjson::Value mapToJsonObj(const TagValues &tags, rapidjson::Document::AllocatorType &allocator)
 {
@@ -118,24 +118,23 @@ bool compareTags(const TagValues &tags1, const TagValues &tags2, const Tags &tag
   return true;
 }
 
+bool isSameVolume(const TagValues &tags1, const TagValues &tags2)
+{
+  // TODO check cosines  
+  return compareTags(tags1, tags2, {SERIES_INSTANCE_UID});
+}
+
 Volumes groupByVolume(const DicomFiles &dicomFiles)
 {
   Volumes volumes;
   for (const DicomFile &dicomFile : dicomFiles)
   {
-    const auto &[file, tags] = dicomFile;
-    Volume *matchingVolume = nullptr;
-    for (Volume &volume : volumes)
-    {
-      const DicomFile fileInVolume = *volume.begin();
-      const TagValues volumeTags = fileInVolume.second;
-      if (compareTags(volumeTags, tags, {SERIES_INSTANCE_UID}))
-      {
-        matchingVolume = &volume;
-        break;
-      }
-    }
-    if (matchingVolume)
+    const auto tags = dicomFile.second;
+    auto matchingVolume = std::find_if(volumes.begin(), volumes.end(), [&tags](const Volume &volume)
+                                       {
+      return isSameVolume(volume.begin()->second, tags); });
+
+    if (matchingVolume != volumes.end())
     {
       matchingVolume->push_back(dicomFile);
     }
@@ -154,17 +153,11 @@ ImageSets groupByImageSet(const Volumes &volumes)
   for (const Volume &volume : volumes)
   {
     const auto volumeTags = volume.begin()->second;
-    Volumes *matchingImageSet = nullptr;
-    for (Volumes &volumes : imageSets)
-    {
+    auto matchingImageSet = std::find_if(imageSets.begin(), imageSets.end(), [&volumeTags](const Volumes &volumes)
+                                         {
       const TagValues imageSetTags = volumes.begin()->begin()->second;
-      if (compareTags(imageSetTags, volumeTags, {STUDY_INSTANCE_UID}))
-      {
-        matchingImageSet = &volumes;
-        break;
-      }
-    }
-    if (matchingImageSet)
+      return compareTags(imageSetTags, volumeTags, {STUDY_INSTANCE_UID}); });
+    if (matchingImageSet != imageSets.end())
     {
       matchingImageSet->push_back(volume);
     }
@@ -177,6 +170,33 @@ ImageSets groupByImageSet(const Volumes &volumes)
   return imageSets;
 }
 
+Volumes sortSpatially(Volumes &volumes)
+{
+  Volumes sortedVolumes;
+  for (Volume &volume : volumes)
+  {
+    std::vector<std::string> unsortedSerieFileNames;
+    for (const DicomFile &dicomFile : volume)
+    {
+      unsortedSerieFileNames.push_back(dicomFile.first);
+    }
+    std::vector<std::string> sortedFileNames = sortSpatially(unsortedSerieFileNames);
+
+    Volume sorted;
+    for (const auto &fileName : sortedFileNames)
+    {
+      const auto matchingDicomFile = std::find_if(volume.begin(), volume.end(), [&fileName](const DicomFile &dicomFile)
+                                                  { return dicomFile.first == fileName; });
+      if (matchingDicomFile != volume.end())
+      {
+        sorted.push_back(*matchingDicomFile);
+      }
+    }
+    sortedVolumes.push_back(sorted);
+  }
+  return sortedVolumes;
+}
+
 rapidjson::Document toJson(const ImageSets &imageSets)
 {
   rapidjson::Document imageSetsJson(rapidjson::kArrayType);
@@ -271,7 +291,8 @@ int main(int argc, char *argv[])
   ITK_WASM_PARSE(pipeline);
 
   const DicomFiles dicomFiles = readTags(files);
-  const Volumes volumes = groupByVolume(dicomFiles);
+  Volumes volumes = groupByVolume(dicomFiles);
+  volumes = sortSpatially(volumes);
   const ImageSets imageSets = groupByImageSet(volumes);
 
   rapidjson::Document imageSetsJson = toJson(imageSets);

packages/dicom/gdcm/read-image-dicom-file-series.cxx

@@ -20,6 +20,10 @@
 #include <fstream>
 #include <memory>
 
+#include "rapidjson/document.h"
+#include "rapidjson/prettywriter.h"
+#include "rapidjson/ostreamwrapper.h"
+
 #include "itkCommonEnums.h"
 #include "gdcmSerieHelper.h"
 #include "itkImageIOBase.h"
@@ -32,18 +36,7 @@
 #include "itkOutputImage.h"
 #include "itkOutputTextStream.h"
 
-#include "rapidjson/document.h"
-#include "rapidjson/prettywriter.h"
-#include "rapidjson/ostreamwrapper.h"
-
-class CustomSerieHelper: public gdcm::SerieHelper
-{
-public:
-  void AddFileName(std::string const &fileName)
-  {
-    SerieHelper::AddFileName(fileName);
-  }
-};
+#include "SortSpatially.h"
 
 namespace itk
 {
@@ -210,60 +203,7 @@ int runPipeline(itk::wasm::Pipeline & pipeline, std::vector<std::string> & input
 
   if (!singleSortedSeries)
   {
-    std::unique_ptr<CustomSerieHelper> serieHelper(new CustomSerieHelper());
-    for (const std::string & fileName: inputFileNames)
-    {
-      serieHelper->AddFileName(fileName);
-    }
-    serieHelper->SetUseSeriesDetails(true);
-    // Add the default restrictions to refine the file set into multiple series.
-    serieHelper->CreateDefaultUniqueSeriesIdentifier();
-    using SeriesIdContainer = std::vector<std::string>;
-    SeriesIdContainer seriesUIDs;
-    // Accessing the first serie found (assume there is at least one)
-    gdcm::FileList * flist = serieHelper->GetFirstSingleSerieUIDFileSet();
-    while (flist)
-    {
-      if (!flist->empty()) // make sure we have at leat one serie
-      {
-        gdcm::File * file = (*flist)[0]; // for example take the first one
-
-        // Create its unique series ID
-        const std::string id( serieHelper->CreateUniqueSeriesIdentifier(file));
-
-        seriesUIDs.push_back(id);
-      }
-      flist = serieHelper->GetNextSingleSerieUIDFileSet();
-    }
-
-    using FileNamesContainer = std::vector<std::string>;
-    FileNamesContainer fileNames;
-    flist = serieHelper->GetFirstSingleSerieUIDFileSet();
-    const std::string serie = seriesUIDs[0];
-    bool found = false;
-    while (flist && !found)
-    {
-      if (!flist->empty()) // make sure we have at leat one serie
-      {
-        gdcm::File * file = (*flist)[0]; // for example take the first one
-        const std::string id( serieHelper->CreateUniqueSeriesIdentifier(file));
-        if (id == serie)
-        {
-          found = true; // we found a match
-          break;
-        }
-      }
-      flist = serieHelper->GetNextSingleSerieUIDFileSet();
-    }
-    serieHelper->OrderFileList(flist);
-
-    gdcm::FileList::iterator it;
-    for (it = flist->begin(); it != flist->end(); ++it)
-    {
-      gdcm::FileWithName * header = *it;
-      fileNames.push_back(header->filename);
-    }
-
+    std::vector<std::string> fileNames = sortSpatially(inputFileNames);
     reader->SetFileNames(fileNames);
   }
   else