diff --git a/packages/dicom/gdcm/discriminate-volume.cxx b/packages/dicom/gdcm/discriminate-volume.cxx new file mode 100644 index 000000000..d0a1a5db1 --- /dev/null +++ b/packages/dicom/gdcm/discriminate-volume.cxx @@ -0,0 +1,264 @@ +#include "gdcmScanner.h" +#include "gdcmTesting.h" +#include "gdcmIPPSorter.h" +#include "gdcmDirectionCosines.h" +#include + +/* + * 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 SortedFiles; + std::vector 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 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::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::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::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::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 const &GetSortedFiles() const { return SortedFiles; } + std::vector 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; +// } \ No newline at end of file diff --git a/packages/dicom/gdcm/image-sets-normalization.cxx b/packages/dicom/gdcm/image-sets-normalization.cxx index 185639b97..137b6c58a 100644 --- a/packages/dicom/gdcm/image-sets-normalization.cxx +++ b/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 sortedFiles = dv.GetSortedFiles(); + rapidjson::Document imageSetsJson; - imageSetsJson.SetObject(); rapidjson::Document::AllocatorType &allocator = imageSetsJson.GetAllocator(); + imageSetsJson.SetObject(); + + int groupId = 0; + for ( + std::vector::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 writer(stringBuffer); diff --git a/packages/dicom/python/itkwasm-dicom-wasi/itkwasm_dicom_wasi/wasm_modules/image-sets-normalization.wasi.wasm b/packages/dicom/python/itkwasm-dicom-wasi/itkwasm_dicom_wasi/wasm_modules/image-sets-normalization.wasi.wasm index 5d98c69e4..66b15e5fc 100755 Binary files a/packages/dicom/python/itkwasm-dicom-wasi/itkwasm_dicom_wasi/wasm_modules/image-sets-normalization.wasi.wasm and b/packages/dicom/python/itkwasm-dicom-wasi/itkwasm_dicom_wasi/wasm_modules/image-sets-normalization.wasi.wasm differ diff --git a/packages/dicom/python/itkwasm-dicom-wasi/tests/test_image_sets_normalization.py b/packages/dicom/python/itkwasm-dicom-wasi/tests/test_image_sets_normalization.py index 9c535f69b..087720c9b 100644 --- a/packages/dicom/python/itkwasm-dicom-wasi/tests/test_image_sets_normalization.py +++ b/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 \ No newline at end of file