Sort DICOMs with discriminate volume

packages/dicom/gdcm/discriminate-volume.cxx

@@ -0,0 +1,264 @@
+#include "gdcmScanner.h"
+#include "gdcmTesting.h"
+#include "gdcmIPPSorter.h"
+#include "gdcmDirectionCosines.h"
+#include <cmath>
+
+/*
+ * The following example is a basic sorted which should work in generic cases.
+ * It sort files based on:
+ * Study Instance UID
+ *   Series Instance UID
+ *     Frame of Reference UID
+ *       Image Orientation (Patient)
+ *         Image Position (Patient) (Sorting based on IPP + IOP)
+ */
+
+namespace gdcm
+{
+    const Tag t1(0x0020, 0x000d); // Study Instance UID
+    const Tag t2(0x0020, 0x000e); // Series Instance UID
+    const Tag t3(0x0020, 0x0052); // Frame of Reference UID
+    const Tag t4(0x0020, 0x0037); // Image Orientation (Patient)
+
+    class DiscriminateVolume
+    {
+    private:
+        std::vector<Directory::FilenamesType> SortedFiles;
+        std::vector<Directory::FilenamesType> UnsortedFiles;
+
+        Directory::FilenamesType GetAllFilenamesFromTagToValue(
+            Scanner const &s, Directory::FilenamesType const &filesubset, Tag const &t, const char *valueref)
+        {
+            Directory::FilenamesType theReturn;
+            if (valueref)
+            {
+                size_t len = strlen(valueref);
+                Directory::FilenamesType::const_iterator file = filesubset.begin();
+                for (; file != filesubset.end(); ++file)
+                {
+                    const char *filename = file->c_str();
+                    const char *value = s.GetValue(filename, t);
+                    if (value && strncmp(value, valueref, len) == 0)
+                    {
+                        theReturn.push_back(filename);
+                    }
+                }
+            }
+            return theReturn;
+        }
+
+        void ProcessAIOP(Scanner const &, Directory::FilenamesType const &subset, const char *iopval)
+        {
+            std::cout << "IOP: " << iopval << std::endl;
+            IPPSorter ipp;
+            ipp.SetComputeZSpacing(true);
+            ipp.SetZSpacingTolerance(1e-3); // ??
+            bool b = ipp.Sort(subset);
+            if (!b)
+            {
+                // If you reach here this means you need one more parameter to discriminiat this
+                // series. Eg. T1 / T2 intertwinted. Multiple Echo (0018,0081)
+                std::cerr << "Failed to sort: " << subset.begin()->c_str() << std::endl;
+                for (
+                    Directory::FilenamesType::const_iterator file = subset.begin();
+                    file != subset.end(); ++file)
+                {
+                    std::cerr << *file << std::endl;
+                }
+                UnsortedFiles.push_back(subset);
+                return;
+            }
+            ipp.Print(std::cout);
+            SortedFiles.push_back(ipp.GetFilenames());
+        }
+
+        void ProcessAFrameOfRef(Scanner const &s, Directory::FilenamesType const &subset, const char *frameuid)
+        {
+            // In this subset of files (belonging to same series), let's find those
+            // belonging to the same Frame ref UID:
+            Directory::FilenamesType files = GetAllFilenamesFromTagToValue(
+                s, subset, t3, frameuid);
+
+            std::set<std::string> iopset;
+
+            for (
+                Directory::FilenamesType::const_iterator file = files.begin();
+                file != files.end(); ++file)
+            {
+                // std::cout << *file << std::endl;
+                const char *value = s.GetValue(file->c_str(), gdcm::t4);
+                assert(value);
+                iopset.insert(value);
+            }
+            size_t n = iopset.size();
+            if (n == 0)
+            {
+                assert(files.empty());
+                return;
+            }
+
+            std::cout << "Frame of Ref: " << frameuid << std::endl;
+            if (n == 1)
+            {
+                ProcessAIOP(s, files, iopset.begin()->c_str());
+            }
+            else
+            {
+                const char *f = files.begin()->c_str();
+                std::cerr << "More than one IOP: " << f << std::endl;
+                // Make sure that there is actually 'n' different IOP
+                gdcm::DirectionCosines ref;
+                gdcm::DirectionCosines dc;
+                for (
+                    std::set<std::string>::const_iterator it = iopset.begin();
+                    it != iopset.end(); ++it)
+                {
+                    ref.SetFromString(it->c_str());
+                    for (
+                        Directory::FilenamesType::const_iterator file = files.begin();
+                        file != files.end(); ++file)
+                    {
+                        std::string value = s.GetValue(file->c_str(), gdcm::t4);
+                        if (value != it->c_str())
+                        {
+                            dc.SetFromString(value.c_str());
+                            const double crossdot = ref.CrossDot(dc);
+                            const double eps = std::fabs(1. - crossdot);
+                            if (eps < 1e-6)
+                            {
+                                std::cerr << "Problem with IOP discrimination: " << file->c_str()
+                                          << " " << it->c_str() << std::endl;
+                                return;
+                            }
+                        }
+                    }
+                }
+                // If we reach here this means there is actually 'n' different IOP
+                for (
+                    std::set<std::string>::const_iterator it = iopset.begin();
+                    it != iopset.end(); ++it)
+                {
+                    const char *iopvalue = it->c_str();
+                    Directory::FilenamesType iopfiles = GetAllFilenamesFromTagToValue(
+                        s, files, t4, iopvalue);
+                    ProcessAIOP(s, iopfiles, iopvalue);
+                }
+            }
+        }
+
+        void ProcessASeries(Scanner const &s, const char *seriesuid)
+        {
+            std::cout << "Series: " << seriesuid << std::endl;
+            // let's find all files belonging to this series:
+            Directory::FilenamesType seriesfiles = GetAllFilenamesFromTagToValue(
+                s, s.GetFilenames(), t2, seriesuid);
+
+            gdcm::Scanner::ValuesType vt3 = s.GetValues(t3);
+            for (
+                gdcm::Scanner::ValuesType::const_iterator it = vt3.begin(); it != vt3.end(); ++it)
+            {
+                ProcessAFrameOfRef(s, seriesfiles, it->c_str());
+            }
+        }
+
+        void ProcessAStudy(Scanner const &s, const char *studyuid)
+        {
+            std::cout << "Study: " << studyuid << std::endl;
+            gdcm::Scanner::ValuesType vt2 = s.GetValues(t2);
+            for (
+                gdcm::Scanner::ValuesType::const_iterator it = vt2.begin(); it != vt2.end(); ++it)
+            {
+                ProcessASeries(s, it->c_str());
+            }
+        }
+
+    public:
+        void Print(std::ostream &os)
+        {
+            os << "Sorted Files: " << std::endl;
+            for (
+                std::vector<Directory::FilenamesType>::const_iterator it = SortedFiles.begin();
+                it != SortedFiles.end(); ++it)
+            {
+                os << "Group: " << std::endl;
+                for (
+                    Directory::FilenamesType::const_iterator file = it->begin();
+                    file != it->end(); ++file)
+                {
+                    os << *file << std::endl;
+                }
+            }
+            os << "Unsorted Files: " << std::endl;
+            for (
+                std::vector<Directory::FilenamesType>::const_iterator it = UnsortedFiles.begin();
+                it != UnsortedFiles.end(); ++it)
+            {
+                os << "Group: " << std::endl;
+                for (
+                    Directory::FilenamesType::const_iterator file = it->begin();
+                    file != it->end(); ++file)
+                {
+                    os << *file << std::endl;
+                }
+            }
+        }
+
+        std::vector<Directory::FilenamesType> const &GetSortedFiles() const { return SortedFiles; }
+        std::vector<Directory::FilenamesType> const &GetUnsortedFiles() const { return UnsortedFiles; }
+
+        void ProcessIntoVolume(Scanner const &s)
+        {
+            gdcm::Scanner::ValuesType vt1 = s.GetValues(gdcm::t1);
+            for (
+                gdcm::Scanner::ValuesType::const_iterator it = vt1.begin(); it != vt1.end(); ++it)
+            {
+                ProcessAStudy(s, it->c_str());
+            }
+        }
+    };
+
+} // namespace gdcm
+
+// int main(int argc, char *argv[])
+// {
+//     std::string dir1;
+//     if (argc < 2)
+//     {
+//         const char *extradataroot = nullptr;
+// #ifdef GDCM_BUILD_TESTING
+//         extradataroot = gdcm::Testing::GetDataExtraRoot();
+// #endif
+//         if (!extradataroot)
+//         {
+//             return 1;
+//         }
+//         dir1 = extradataroot;
+//         dir1 += "/gdcmSampleData/ForSeriesTesting/VariousIncidences/ST1";
+//     }
+//     else
+//     {
+//         dir1 = argv[1];
+//     }
+
+//     gdcm::Directory d;
+//     d.Load(dir1, true); // recursive !
+
+//     gdcm::Scanner s;
+//     s.AddTag(gdcm::t1);
+//     s.AddTag(gdcm::t2);
+//     s.AddTag(gdcm::t3);
+//     s.AddTag(gdcm::t4);
+//     bool b = s.Scan(d.GetFilenames());
+//     if (!b)
+//     {
+//         std::cerr << "Scanner failed" << std::endl;
+//         return 1;
+//     }
+
+//     gdcm::DiscriminateVolume dv;
+//     dv.ProcessIntoVolume(s);
+//     dv.Print(std::cout);
+
+//     return 0;
+// }

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

@@ -23,6 +23,8 @@
 #include "rapidjson/stringbuffer.h"
 #include "rapidjson/writer.h"
 
+#include "discriminate-volume.cxx"
+
 int main(int argc, char *argv[])
 {
   itk::wasm::Pipeline pipeline("image-sets-normalization", "Group DICOM files into image sets", argc, argv);
@@ -35,13 +37,57 @@ int main(int argc, char *argv[])
 
   ITK_WASM_PARSE(pipeline);
 
+  gdcm::Scanner s;
+
+  const gdcm::Tag t1(0x0020, 0x000d); // Study Instance UID
+  const gdcm::Tag t2(0x0020, 0x000e); // Series Instance UID
+  const gdcm::Tag t3(0x0020, 0x0052); // Frame of Reference UID
+  const gdcm::Tag t4(0x0020, 0x0037); // Image Orientation (Patient)
+
+  s.AddTag(t1);
+  s.AddTag(t2);
+  s.AddTag(t3);
+  s.AddTag(t4);
+
+  bool b = s.Scan(files);
+  if (!b)
+  {
+    std::cerr << "Scanner failed" << std::endl;
+    return 1;
+  }
+
+  gdcm::DiscriminateVolume dv;
+  dv.ProcessIntoVolume(s);
+
+  std::vector<gdcm::Directory::FilenamesType> sortedFiles = dv.GetSortedFiles();
+
   rapidjson::Document imageSetsJson;
-  imageSetsJson.SetObject();
   rapidjson::Document::AllocatorType &allocator = imageSetsJson.GetAllocator();
+  imageSetsJson.SetObject();
+
+  int groupId = 0;
+  for (
+      std::vector<gdcm::Directory::FilenamesType>::const_iterator fileNames = sortedFiles.begin();
+      fileNames != sortedFiles.end(); ++fileNames)
+  {
+    groupId++;
+    rapidjson::Value sortedFileNameArray(rapidjson::kArrayType);
+
+    for (
+        gdcm::Directory::FilenamesType::const_iterator fileName = fileNames->begin();
+        fileName != fileNames->end(); ++fileName)
+    {
+      rapidjson::Value fileNameValue;
+      fileNameValue.SetString(fileName->c_str(), fileName->size(), allocator);
+      sortedFileNameArray.PushBack(fileNameValue, allocator);
+    }
+
+    std::string fileName = fileNames->front();
+    const char* studyInstanceUID = s.GetValue(fileName.c_str(), t1);
 
-  rapidjson::Value almostEqualValue;
-  almostEqualValue.SetBool(false);
-  imageSetsJson.AddMember("", almostEqualValue, allocator);
+    rapidjson::Value groupIdKey(studyInstanceUID, allocator);
+    imageSetsJson.AddMember(groupIdKey, sortedFileNameArray, allocator);
+  }
 
   rapidjson::StringBuffer stringBuffer;
   rapidjson::Writer<rapidjson::StringBuffer> writer(stringBuffer);

packages/dicom/python/itkwasm-dicom-wasi/tests/test_image_sets_normalization.py

@@ -3,14 +3,35 @@
 from .common import test_input_path, test_output_path
 
 
-def test_image_sets_normalization():
+def test_one_series():
     files = [
+        test_input_path / "DicomImageOrientationTest" / "ImageOrientation.3.dcm",
         test_input_path / "DicomImageOrientationTest" / "ImageOrientation.1.dcm",
         test_input_path / "DicomImageOrientationTest" / "ImageOrientation.2.dcm",
-        test_input_path / "DicomImageOrientationTest" / "ImageOrientation.3.dcm",
     ]
 
     assert files[0].exists()
 
-    output_text = image_sets_normalization(files)
-    assert output_text
+    image_sets = image_sets_normalization(files)
+    print(image_sets)
+    # assert image_sets
+
+    assert image_sets == [
+        str(files[1]),
+        str(files[2]),
+        str(files[0]),
+    ]
+
+# def test_two_series():
+#     files = [
+#         test_input_path / "DicomImageOrientationTest" / "ImageOrientation.3.dcm",
+#         test_input_path / "DicomImageOrientationTest" / "ImageOrientation.1.dcm",
+#         test_input_path / "DicomImageOrientationTest" / "ImageOrientation.2.dcm",
+#         test_input_path / "dicom-images" / "CT" / "1-1.dcm",
+#         test_input_path / "dicom-images" / "CT" / "1-2.dcm",
+#     ]
+
+#     assert files[0].exists()
+
+#     image_sets = image_sets_normalization(files)
+#     assert image_sets