diff --git a/README.md b/README.md index 06b3dfe88e..32d131e2e9 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,8 @@ # AGAVE : Advanced GPU Accelerated Volume Explorer -AGAVE is a desktop application for viewing multichannel volume data. Several formats are supported, including OME-ZARR 0.4 and 0.5, OME-TIFF and Zeiss .czi files. +AGAVE is a desktop application for viewing multichannel volume data. Several formats are supported, including OME-ZARR 0.4 and 0.5, OME-TIFF, Zeiss .czi, and grayscale DICOM files (.dcm, .dicom, .ima). + +The DICOM reader uses DCMTK 3.7.0. It opens single 8-bit or 16-bit grayscale DICOM objects and folders of DICOM slices, including nested study/series folders. Each `SeriesInstanceUID` is a selectable scene. Slices are ordered by image position and orientation, then by instance number or filename when spatial metadata is absent; irregular stacks are rejected. Multi-frame objects become Z stacks. The reader decodes uncompressed, RLE, JPEG, and JPEG-LS transfer syntaxes supported by DCMTK and maps displayed intensity into Agave's unsigned 16-bit volume format. Color images and JPEG 2000 are not yet supported. ![screenshot](https://github.com/user-attachments/assets/b96618f2-7020-4b93-936e-9b32b795ea83) diff --git a/agave_app/agaveGui.cpp b/agave_app/agaveGui.cpp index b0ed4b1cbe..162d3122f1 100644 --- a/agave_app/agaveGui.cpp +++ b/agave_app/agaveGui.cpp @@ -223,7 +223,7 @@ agaveGui::createActions() m_openUrlAction->setStatusTip(tr("Open an existing volume in the cloud")); connect(m_openUrlAction, SIGNAL(triggered()), this, SLOT(openUrl())); - m_openDirectoryAction = new QAction(tr("&Open directory (.zarr)"), this); + m_openDirectoryAction = new QAction(tr("&Open directory (.zarr, DICOM)"), this); m_openDirectoryAction->setShortcuts(QKeySequence::Open); m_openDirectoryAction->setStatusTip(tr("Open an existing volume from local directory")); connect(m_openDirectoryAction, SIGNAL(triggered()), this, SLOT(openDirectory())); diff --git a/renderlib/io/CMakeLists.txt b/renderlib/io/CMakeLists.txt index 25485464ac..3996af847c 100644 --- a/renderlib/io/CMakeLists.txt +++ b/renderlib/io/CMakeLists.txt @@ -15,6 +15,8 @@ target_sources(renderlib PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/FileReaderCCP4.h" "${CMAKE_CURRENT_SOURCE_DIR}/FileReaderCzi.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/FileReaderCzi.h" +"${CMAKE_CURRENT_SOURCE_DIR}/FileReaderDicom.cpp" +"${CMAKE_CURRENT_SOURCE_DIR}/FileReaderDicom.h" "${CMAKE_CURRENT_SOURCE_DIR}/FileReaderImageSequence.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/FileReaderImageSequence.h" "${CMAKE_CURRENT_SOURCE_DIR}/FileReaderTIFF.cpp" @@ -41,6 +43,26 @@ add_subdirectory(libCZI) # end libczi dependency +# Build only the DCMTK modules needed for grayscale image decoding. The JPEG +# and JPEG-LS modules also link to dcmimage, even for grayscale images. +set(DCMTK_MODULES "ofstd;oflog;oficonv;dcmdata;dcmimgle;dcmimage;dcmjpeg;dcmjpls" CACHE STRING "DCMTK modules" FORCE) +set(DCMTK_DEFAULT_DICT builtin CACHE STRING "Embed the DCMTK dictionary for standalone Agave builds" FORCE) +set(BUILD_APPS OFF CACHE BOOL "Build DCMTK command line applications" FORCE) +set(DCMTK_WITH_OPENSSL OFF CACHE BOOL "DCMTK TLS support" FORCE) +set(DCMTK_WITH_ICONV OFF CACHE BOOL "DCMTK iconv support" FORCE) +set(DCMTK_WITH_ICU OFF CACHE BOOL "DCMTK ICU support" FORCE) +set(DCMTK_WITH_PNG OFF CACHE BOOL "DCMTK PNG support" FORCE) +set(DCMTK_WITH_TIFF OFF CACHE BOOL "DCMTK TIFF support" FORCE) +set(DCMTK_WITH_XML OFF CACHE BOOL "DCMTK XML support" FORCE) +set(DCMTK_WITH_OPENJPEG OFF CACHE BOOL "DCMTK JPEG 2000 support" FORCE) +set(DCMTK_WITH_DOXYGEN OFF CACHE BOOL "DCMTK API documentation" FORCE) +FetchContent_Declare( + dcmtk + URL "https://github.com/DCMTK/dcmtk/archive/refs/tags/DCMTK-3.7.0.tar.gz" + URL_HASH SHA256=5bb3ec8317dc465788bed2ca789e76d03ae5848c9381cce3b14c1a3f8b6aca56 +) +FetchContent_MakeAvailable(dcmtk) + # Use these system/vcpkg-provided dependencies instead of tensorstore's bundled ones. set(TENSORSTORE_USE_SYSTEM_TIFF ON) set(TENSORSTORE_USE_SYSTEM_ZLIB ON) @@ -56,7 +78,7 @@ FetchContent_MakeAvailable(tensorstore) target_link_libraries(renderlib libCZI + dcmdata dcmimgle dcmimage dcmjpeg dcmjpls tensorstore::tensorstore tensorstore::all_drivers ) - diff --git a/renderlib/io/FileReader.cpp b/renderlib/io/FileReader.cpp index 8c2e14df45..81a09d3543 100644 --- a/renderlib/io/FileReader.cpp +++ b/renderlib/io/FileReader.cpp @@ -2,6 +2,7 @@ #include "FileReaderCCP4.h" #include "FileReaderCzi.h" +#include "FileReaderDicom.h" #include "FileReaderImageSequence.h" #include "FileReaderTIFF.h" #include "FileReaderTestVolume.h" @@ -54,16 +55,22 @@ FileReader::getReader(const std::string& filepath, bool isImageSequence) return new FileReaderTIFF(filepath); } else if (extstr == ".czi") { return new FileReaderCzi(filepath); + } else if (extstr == ".dcm" || extstr == ".dicom" || extstr == ".ima") { + return new FileReaderDicom(filepath); } else if (extstr == ".map" || extstr == ".mrc") { return new FileReaderCCP4(filepath); } else if (extstr == ".zarr") { return new FileReaderZarr(filepath); } - // if it's a directory, and contains the string zarr anywhere, we assume it's a zarr + // A directory with Zarr metadata is a Zarr store; other directories can + // contain DICOM series, including files without a DICOM extension. else if (std::filesystem::is_directory(filepath)) { - if (filepath.find("zarr") != std::string::npos) { + const std::filesystem::path directory(filepath); + if (std::filesystem::exists(directory / "zarr.json") || std::filesystem::exists(directory / ".zgroup") || + std::filesystem::exists(directory / ".zarray") || std::filesystem::exists(directory / ".zattrs")) { return new FileReaderZarr(filepath); } + return new FileReaderDicom(filepath); } return nullptr; diff --git a/renderlib/io/FileReaderDicom.cpp b/renderlib/io/FileReaderDicom.cpp new file mode 100644 index 0000000000..b24cc396fe --- /dev/null +++ b/renderlib/io/FileReaderDicom.cpp @@ -0,0 +1,468 @@ +#include "FileReaderDicom.h" + +#include "ImageXYZC.h" +#include "Logging.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +struct DicomFolderIndex +{ + struct Slice + { + std::string path; + uint32_t width = 0; + uint32_t height = 0; + uint32_t frames = 1; + uint16_t bits = 0; + std::string photometric; + float spacingX = 1.0f; + float spacingY = 1.0f; + float spacingZ = 1.0f; + std::array position{}; + std::array orientation{}; + bool hasGeometry = false; + int instance = 0; + bool hasInstance = false; + double projection = 0.0; + }; + + struct Series + { + std::vector slices; + VolumeDimensions dimensions; + }; + + std::vector series; +}; + +namespace { + +void +registerDecoders() +{ + static std::once_flag once; + std::call_once(once, [] { + DcmRLEDecoderRegistration::registerCodecs(); + DJDecoderRegistration::registerCodecs(); + DJLSDecoderRegistration::registerCodecs(); + }); +} + +float +positiveSpacing(DcmDataset* dataset, const DcmTagKey& tag, unsigned long index = 0) +{ + Float64 value = 0.0; + return dataset->findAndGetFloat64(tag, value, index).good() && std::isfinite(value) && value > 0.0 + ? static_cast(value) + : 1.0f; +} + +float +sliceSpacing(DcmDataset* dataset) +{ + Float64 value = 0.0; + if (dataset->findAndGetFloat64(DCM_SpacingBetweenSlices, value).good() && std::isfinite(value) && value > 0.0) { + return static_cast(value); + } + return positiveSpacing(dataset, DCM_SliceThickness); +} + +bool +isDirectory(const std::string& filepath) +{ + std::error_code error; + return std::filesystem::is_directory(filepath, error); +} + +bool +readSliceMetadata(const std::filesystem::path& path, std::string& seriesUid, DicomFolderIndex::Slice& slice) +{ + DcmFileFormat file; + if (file.loadFile(path.string().c_str()).bad()) { + return false; + } + DcmDataset* dataset = file.getDataset(); + OFString uid; + OFString photometric; + Uint16 rows = 0; + Uint16 columns = 0; + Uint16 samples = 0; + Uint16 bits = 0; + if (dataset->findAndGetOFString(DCM_SeriesInstanceUID, uid).bad() || uid.empty() || + dataset->findAndGetUint16(DCM_Rows, rows).bad() || rows == 0 || + dataset->findAndGetUint16(DCM_Columns, columns).bad() || columns == 0 || + dataset->findAndGetUint16(DCM_SamplesPerPixel, samples).bad() || samples != 1 || + dataset->findAndGetUint16(DCM_BitsAllocated, bits).bad() || (bits != 8 && bits != 16) || + dataset->findAndGetOFString(DCM_PhotometricInterpretation, photometric).bad() || + (photometric != "MONOCHROME1" && photometric != "MONOCHROME2") || !dataset->tagExists(DCM_PixelData)) { + return false; + } + + slice.path = path.string(); + slice.width = columns; + slice.height = rows; + slice.bits = bits; + slice.photometric = photometric.c_str(); + slice.spacingY = positiveSpacing(dataset, DCM_PixelSpacing, 0); + slice.spacingX = positiveSpacing(dataset, DCM_PixelSpacing, 1); + slice.spacingZ = sliceSpacing(dataset); + OFString frames; + if (dataset->findAndGetOFString(DCM_NumberOfFrames, frames).good()) { + try { + const unsigned long count = std::stoul(frames.c_str()); + if (count > std::numeric_limits::max()) { + return false; + } + slice.frames = static_cast(count); + } catch (const std::exception&) { + return false; + } + if (slice.frames == 0) { + return false; + } + } + slice.hasGeometry = true; + for (unsigned long i = 0; i < 3; ++i) { + slice.hasGeometry &= dataset->findAndGetFloat64(DCM_ImagePositionPatient, slice.position[i], i).good(); + slice.hasGeometry &= std::isfinite(slice.position[i]); + } + for (unsigned long i = 0; i < 6; ++i) { + slice.hasGeometry &= dataset->findAndGetFloat64(DCM_ImageOrientationPatient, slice.orientation[i], i).good(); + slice.hasGeometry &= std::isfinite(slice.orientation[i]); + } + OFString instance; + if (dataset->findAndGetOFString(DCM_InstanceNumber, instance).good()) { + try { + slice.instance = std::stoi(instance.c_str()); + slice.hasInstance = true; + } catch (const std::exception&) { + } + } + seriesUid = uid.c_str(); + return true; +} + +std::shared_ptr +scanFolder(const std::string& filepath) +{ + auto index = std::make_shared(); + std::map> grouped; + std::error_code error; + std::filesystem::recursive_directory_iterator it( + filepath, std::filesystem::directory_options::skip_permission_denied, error); + const std::filesystem::recursive_directory_iterator end; + for (; !error && it != end; it.increment(error)) { + if (!it->is_regular_file(error)) { + error.clear(); + continue; + } + std::string uid; + DicomFolderIndex::Slice slice; + if (readSliceMetadata(it->path(), uid, slice)) { + grouped[uid].push_back(std::move(slice)); + } + } + + for (auto& group : grouped) { + const auto& uid = group.first; + auto& slices = group.second; + if (slices.empty() || slices.size() > std::numeric_limits::max()) { + continue; + } + const auto first = slices.front(); + bool consistent = std::all_of(slices.begin(), slices.end(), [&](const auto& slice) { + return slice.width == first.width && slice.height == first.height && slice.bits == first.bits && + slice.photometric == first.photometric && std::abs(slice.spacingX - first.spacingX) < 0.001f && + std::abs(slice.spacingY - first.spacingY) < 0.001f && (slices.size() == 1 || slice.frames == 1); + }); + if (!consistent) { + LOG_WARNING << "Skipping DICOM series with incompatible slices: " << uid; + continue; + } + + float spacingZ = first.spacingZ; + const bool spatial = std::all_of(slices.begin(), slices.end(), [](const auto& slice) { return slice.hasGeometry; }); + if (spatial && slices.size() > 1) { + const auto& o = first.orientation; + std::array normal = { o[1] * o[5] - o[2] * o[4], + o[2] * o[3] - o[0] * o[5], + o[0] * o[4] - o[1] * o[3] }; + const double norm = std::sqrt(normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]); + if (norm <= 0.9) { + LOG_WARNING << "Skipping DICOM series with invalid orientation: " << uid; + continue; + } + for (auto& slice : slices) { + for (size_t i = 0; i < 6; ++i) { + consistent &= std::abs(slice.orientation[i] - o[i]) < 0.001; + } + slice.projection = + (slice.position[0] * normal[0] + slice.position[1] * normal[1] + slice.position[2] * normal[2]) / norm; + } + if (!consistent) { + LOG_WARNING << "Skipping DICOM series with inconsistent orientation: " << uid; + continue; + } + std::sort(slices.begin(), slices.end(), [](const auto& a, const auto& b) { return a.projection < b.projection; }); + const double step = slices[1].projection - slices[0].projection; + consistent = step > 0.001; + for (size_t i = 2; i < slices.size(); ++i) { + const double delta = slices[i].projection - slices[i - 1].projection; + consistent &= delta > 0.001 && std::abs(delta - step) <= std::max(0.1, step * 0.05); + } + if (!consistent) { + LOG_WARNING << "Skipping DICOM series with irregular slice positions: " << uid; + continue; + } + spacingZ = static_cast(step); + } else { + const bool numbered = + std::all_of(slices.begin(), slices.end(), [](const auto& slice) { return slice.hasInstance; }); + if (numbered) { + std::sort(slices.begin(), slices.end(), [](const auto& a, const auto& b) { + return a.instance == b.instance ? a.path < b.path : a.instance < b.instance; + }); + if (std::adjacent_find(slices.begin(), slices.end(), [](const auto& a, const auto& b) { + return a.instance == b.instance; + }) != slices.end()) { + LOG_WARNING << "Skipping DICOM series with duplicate instance numbers and no slice positions: " << uid; + continue; + } + } else { + std::sort(slices.begin(), slices.end(), [](const auto& a, const auto& b) { return a.path < b.path; }); + if (slices.size() > 1) { + LOG_WARNING << "DICOM series lacks slice positions and instance numbers; sorting by filename: " << uid; + } + } + } + + DicomFolderIndex::Series series; + series.dimensions.sizeX = first.width; + series.dimensions.sizeY = first.height; + series.dimensions.sizeZ = slices.size() == 1 ? first.frames : static_cast(slices.size()); + series.dimensions.sizeC = 1; + series.dimensions.sizeT = 1; + series.dimensions.bitsPerPixel = ImageXYZC::IN_MEMORY_BPP; + series.dimensions.channelNames = { "Intensity" }; + series.dimensions.spatialUnits = "mm"; + series.dimensions.physicalSizeX = first.spacingX; + series.dimensions.physicalSizeY = first.spacingY; + series.dimensions.physicalSizeZ = spacingZ; + series.slices = std::move(slices); + index->series.push_back(std::move(series)); + } + return index; +} + +bool +readDimensions(const std::string& filepath, VolumeDimensions& dims) +{ + registerDecoders(); + DcmFileFormat file; + if (file.loadFile(filepath.c_str()).bad()) { + LOG_ERROR << "Could not open DICOM file: " << filepath; + return false; + } + DcmDataset* dataset = file.getDataset(); + Uint16 samples = 0; + Uint16 bits = 0; + OFString photometric; + if (dataset->findAndGetUint16(DCM_SamplesPerPixel, samples).bad() || samples != 1 || + dataset->findAndGetUint16(DCM_BitsAllocated, bits).bad() || (bits != 8 && bits != 16) || + dataset->findAndGetOFString(DCM_PhotometricInterpretation, photometric).bad() || + (photometric != "MONOCHROME1" && photometric != "MONOCHROME2")) { + LOG_ERROR << "DICOM reader supports 8 or 16 bit grayscale images: " << filepath; + return false; + } + + DicomImage image(filepath.c_str()); + if (image.getStatus() != EIS_Normal || image.getWidth() == 0 || image.getHeight() == 0 || + image.getFrameCount() == 0 || image.getWidth() > std::numeric_limits::max() || + image.getHeight() > std::numeric_limits::max() || + image.getFrameCount() > std::numeric_limits::max()) { + LOG_ERROR << "Could not decode DICOM image: " << filepath; + return false; + } + + dims.sizeX = static_cast(image.getWidth()); + dims.sizeY = static_cast(image.getHeight()); + dims.sizeZ = static_cast(image.getFrameCount()); + dims.sizeC = 1; + dims.sizeT = 1; + dims.bitsPerPixel = ImageXYZC::IN_MEMORY_BPP; + dims.channelNames = { "Intensity" }; + dims.spatialUnits = "mm"; + // Pixel Spacing stores row spacing before column spacing. + dims.physicalSizeY = positiveSpacing(dataset, DCM_PixelSpacing, 0); + dims.physicalSizeX = positiveSpacing(dataset, DCM_PixelSpacing, 1); + dims.physicalSizeZ = sliceSpacing(dataset); + return true; +} + +} // namespace + +FileReaderDicom::FileReaderDicom(const std::string& filepath) +{ + (void)filepath; +} + +FileReaderDicom::~FileReaderDicom() = default; + +std::shared_ptr +FileReaderDicom::folderIndex(const std::string& filepath) +{ + std::lock_guard lock(m_indexMutex); + if (!m_index || m_indexedPath != filepath) { + m_index = scanFolder(filepath); + m_indexedPath = filepath; + } + return m_index; +} + +uint32_t +FileReaderDicom::loadNumScenes(const std::string& filepath) +{ + if (isDirectory(filepath)) { + return static_cast(folderIndex(filepath)->series.size()); + } + VolumeDimensions dims; + return readDimensions(filepath, dims) ? 1 : 0; +} + +VolumeDimensions +FileReaderDicom::loadDimensions(const std::string& filepath, uint32_t scene) +{ + if (isDirectory(filepath)) { + const auto index = folderIndex(filepath); + return scene < index->series.size() ? index->series[scene].dimensions : VolumeDimensions{}; + } + VolumeDimensions dims; + if (scene != 0 || !readDimensions(filepath, dims)) { + return {}; + } + return dims; +} + +std::vector +FileReaderDicom::loadMultiscaleDims(const std::string& filepath, uint32_t scene) +{ + if (scene >= loadNumScenes(filepath)) { + return {}; + } + const VolumeDimensions dims = loadDimensions(filepath, scene); + MultiscaleDims scale; + scale.shape = { 1, 1, dims.sizeZ, dims.sizeY, dims.sizeX }; + scale.scale = { 1.0f, 1.0f, dims.physicalSizeZ, dims.physicalSizeY, dims.physicalSizeX }; + scale.dimensionOrder = { "T", "C", "Z", "Y", "X" }; + scale.dtype = "uint16"; + scale.channelNames = dims.channelNames; + scale.spatialUnits = dims.spatialUnits; + return { scale }; +} + +std::shared_ptr +FileReaderDicom::loadVolumeBlocking(const LoadSpec& loadSpec, LoadProgress& progress) +{ + if (loadSpec.time != 0 || + (!loadSpec.channels.empty() && (loadSpec.channels.size() != 1 || loadSpec.channels[0] != 0))) { + return {}; + } + VolumeDimensions dims; + std::vector files; + if (isDirectory(loadSpec.filepath)) { + const auto index = folderIndex(loadSpec.filepath); + if (loadSpec.scene >= index->series.size()) { + return {}; + } + const auto& series = index->series[loadSpec.scene]; + dims = series.dimensions; + for (const auto& slice : series.slices) { + files.push_back(slice.path); + } + } else { + if (loadSpec.scene != 0 || !readDimensions(loadSpec.filepath, dims)) { + return {}; + } + files.push_back(loadSpec.filepath); + } + const size_t planePixels = static_cast(dims.sizeX) * dims.sizeY; + if (planePixels > std::numeric_limits::max() / sizeof(uint16_t) / dims.sizeZ) { + LOG_ERROR << "DICOM image is too large: " << loadSpec.filepath; + return {}; + } + double globalMin = std::numeric_limits::infinity(); + double globalMax = -std::numeric_limits::infinity(); + for (const auto& path : files) { + if (progress.isCancelled()) { + return {}; + } + DicomImage image(path.c_str()); + double imageMin = 0.0; + double imageMax = 0.0; + if (image.getStatus() != EIS_Normal || !image.getMinMaxValues(imageMin, imageMax)) { + LOG_ERROR << "Could not decode DICOM image: " << path; + return {}; + } + globalMin = std::min(globalMin, imageMin); + globalMax = std::max(globalMax, imageMax); + } + auto data = std::make_unique(planePixels * dims.sizeZ * sizeof(uint16_t)); + const double width = std::max(1.0, globalMax - globalMin + 1.0); + const double center = (globalMin + globalMax + 1.0) / 2.0; + uint32_t z = 0; + for (const auto& path : files) { + DicomImage image(path.c_str()); + if (image.getStatus() != EIS_Normal || !image.setWindow(center, width) || image.getWidth() != dims.sizeX || + image.getHeight() != dims.sizeY || image.getFrameCount() > dims.sizeZ - z) { + LOG_ERROR << "DICOM series changed while loading: " << path; + return {}; + } + for (unsigned long frameIndex = 0; frameIndex < image.getFrameCount(); ++frameIndex, ++z) { + if (progress.isCancelled()) { + return {}; + } + const void* frame = image.getOutputData(16, frameIndex); + if (!frame) { + LOG_ERROR << "Could not decode DICOM frame " << frameIndex << " in " << path; + return {}; + } + std::memcpy( + data.get() + static_cast(z) * planePixels * sizeof(uint16_t), frame, planePixels * sizeof(uint16_t)); + progress.setProgress(z + 1, dims.sizeZ); + } + } + if (z != dims.sizeZ) { + LOG_ERROR << "DICOM series changed while loading: " << loadSpec.filepath; + return {}; + } + auto result = std::make_shared(dims.sizeX, + dims.sizeY, + dims.sizeZ, + 1, + ImageXYZC::IN_MEMORY_BPP, + data.release(), + dims.physicalSizeX, + dims.physicalSizeY, + dims.physicalSizeZ, + dims.spatialUnits); + result->setChannelNames(dims.channelNames); + return result; +} diff --git a/renderlib/io/FileReaderDicom.h b/renderlib/io/FileReaderDicom.h new file mode 100644 index 0000000000..1d533112c6 --- /dev/null +++ b/renderlib/io/FileReaderDicom.h @@ -0,0 +1,31 @@ +#pragma once + +#include "BlockingFileReader.h" +#include "VolumeDimensions.h" + +#include +#include +#include + +struct DicomFolderIndex; + +// Reads a grayscale DICOM object or a directory containing one or more series. +class FileReaderDicom : public BlockingFileReader +{ +public: + explicit FileReaderDicom(const std::string& filepath); + ~FileReaderDicom() override; + + bool supportChunkedLoading() const override { return false; } + uint32_t loadNumScenes(const std::string& filepath) override; + VolumeDimensions loadDimensions(const std::string& filepath, uint32_t scene = 0) override; + std::vector loadMultiscaleDims(const std::string& filepath, uint32_t scene = 0) override; + std::shared_ptr loadVolumeBlocking(const LoadSpec& loadSpec, LoadProgress& progress) override; + +private: + std::shared_ptr folderIndex(const std::string& filepath); + + std::mutex m_indexMutex; + std::string m_indexedPath; + std::shared_ptr m_index; +}; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7d7e97642c..e8dafd8244 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -18,6 +18,7 @@ target_sources(agave_test PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/test_cacheManager.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/test_commands.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/test_fileReaderTestVolume.cpp" + "${CMAKE_CURRENT_SOURCE_DIR}/test_fileReaderDicom.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/test_histogram.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/test_main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/test_mathUtil.cpp" diff --git a/test/test_fileReaderDicom.cpp b/test/test_fileReaderDicom.cpp new file mode 100644 index 0000000000..d940cee143 --- /dev/null +++ b/test/test_fileReaderDicom.cpp @@ -0,0 +1,136 @@ +#include "renderlib/ImageXYZC.h" +#include "renderlib/io/FileReader.h" +#include "renderlib/io/FileReaderDicom.h" +#include "renderlib/io/LoadRequest.h" + +#include + +#include +#include +#include + +#include +#include +#include + +namespace { +void +writeSlice(const std::filesystem::path& path, const char* seriesUid, const char* instanceUid, int z, Uint16 base) +{ + DcmFileFormat file; + DcmDataset* dataset = file.getDataset(); + REQUIRE(dataset->putAndInsertString(DCM_SOPClassUID, UID_SecondaryCaptureImageStorage).good()); + REQUIRE(dataset->putAndInsertString(DCM_SOPInstanceUID, instanceUid).good()); + REQUIRE(dataset->putAndInsertString(DCM_SeriesInstanceUID, seriesUid).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_Rows, 2).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_Columns, 2).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_SamplesPerPixel, 1).good()); + REQUIRE(dataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2").good()); + REQUIRE(dataset->putAndInsertUint16(DCM_BitsAllocated, 16).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_BitsStored, 16).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_HighBit, 15).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_PixelRepresentation, 0).good()); + REQUIRE(dataset->putAndInsertString(DCM_PixelSpacing, "0.5\\0.25").good()); + REQUIRE(dataset->putAndInsertString(DCM_ImageOrientationPatient, "1\\0\\0\\0\\1\\0").good()); + const std::string position = "0\\0\\" + std::to_string(z); + REQUIRE(dataset->putAndInsertString(DCM_ImagePositionPatient, position.c_str()).good()); + const std::string instance = std::to_string(20 - z); + REQUIRE(dataset->putAndInsertString(DCM_InstanceNumber, instance.c_str()).good()); + const Uint16 pixels[] = { + base, static_cast(base + 1), static_cast(base + 2), static_cast(base + 3) + }; + REQUIRE(dataset->putAndInsertUint16Array(DCM_PixelData, pixels, 4).good()); + REQUIRE(file.saveFile(path.string().c_str(), EXS_LittleEndianExplicit).good()); +} +} // namespace + +TEST_CASE("DICOM reader loads a multiframe grayscale volume") +{ + const auto path = + std::filesystem::temp_directory_path() / + ("agave-dicom-" + std::to_string(std::chrono::steady_clock::now().time_since_epoch().count()) + ".dcm"); + DcmFileFormat file; + DcmDataset* dataset = file.getDataset(); + REQUIRE(dataset->putAndInsertString(DCM_SOPClassUID, UID_SecondaryCaptureImageStorage).good()); + REQUIRE(dataset->putAndInsertString(DCM_SOPInstanceUID, "1.2.826.0.1.3680043.2.1125.1").good()); + REQUIRE(dataset->putAndInsertUint16(DCM_Rows, 2).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_Columns, 2).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_SamplesPerPixel, 1).good()); + REQUIRE(dataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2").good()); + REQUIRE(dataset->putAndInsertUint16(DCM_BitsAllocated, 16).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_BitsStored, 16).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_HighBit, 15).good()); + REQUIRE(dataset->putAndInsertUint16(DCM_PixelRepresentation, 0).good()); + REQUIRE(dataset->putAndInsertString(DCM_NumberOfFrames, "2").good()); + REQUIRE(dataset->putAndInsertString(DCM_PixelSpacing, "0.5\\0.25").good()); + REQUIRE(dataset->putAndInsertString(DCM_SpacingBetweenSlices, "1.5").good()); + const Uint16 pixels[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 65535 }; + REQUIRE(dataset->putAndInsertUint16Array(DCM_PixelData, pixels, 8).good()); + REQUIRE(file.saveFile(path.string().c_str(), EXS_LittleEndianExplicit).good()); + + std::unique_ptr reader(FileReader::getReader(path.string())); + REQUIRE(dynamic_cast(reader.get()) != nullptr); + CHECK(reader->loadNumScenes(path.string()) == 1); + const auto dims = reader->loadDimensions(path.string()); + CHECK(dims.sizeX == 2); + CHECK(dims.sizeY == 2); + CHECK(dims.sizeZ == 2); + CHECK(dims.sizeC == 1); + CHECK(dims.physicalSizeX == 0.25f); + CHECK(dims.physicalSizeY == 0.5f); + CHECK(dims.physicalSizeZ == 1.5f); + CHECK(dims.spatialUnits == "mm"); + const auto scales = reader->loadMultiscaleDims(path.string()); + REQUIRE(scales.size() == 1); + CHECK(scales[0].shape == std::vector{ 1, 1, 2, 2, 2 }); + + LoadSpec spec; + spec.filepath = path.string(); + LoadProgress progress; + auto image = dynamic_cast(reader.get())->loadVolumeBlocking(spec, progress); + REQUIRE(image != nullptr); + CHECK(image->sizeZ() == 2); + CHECK(image->sizeC() == 1); + CHECK(image->physicalSizeX() == 0.25f); + const auto* first = reinterpret_cast(image->ptr(0, 0)); + const auto* second = reinterpret_cast(image->ptr(0, 1)); + CHECK(first[0] < first[3]); + CHECK(first[3] < second[3]); + CHECK(progress.progress() == 1.0f); + std::filesystem::remove(path); +} + +TEST_CASE("DICOM folder groups series and orders slices by position") +{ + const auto root = + std::filesystem::temp_directory_path() / + ("agave-dicom-folder-" + std::to_string(std::chrono::steady_clock::now().time_since_epoch().count())); + std::filesystem::create_directories(root / "nested"); + // Filenames and instance numbers run opposite to the spatial order. + writeSlice(root / "a.dcm", "1.2.826.0.1.3680043.2.1125.100", "1.2.826.0.1.3680043.2.1125.101", 2, 3000); + writeSlice(root / "nested" / "slice0", "1.2.826.0.1.3680043.2.1125.100", "1.2.826.0.1.3680043.2.1125.102", 0, 1000); + writeSlice(root / "c.dcm", "1.2.826.0.1.3680043.2.1125.100", "1.2.826.0.1.3680043.2.1125.103", 1, 2000); + writeSlice( + root / "nested" / "other.dcm", "1.2.826.0.1.3680043.2.1125.200", "1.2.826.0.1.3680043.2.1125.201", 0, 4000); + + std::unique_ptr reader(FileReader::getReader(root.string())); + REQUIRE(dynamic_cast(reader.get()) != nullptr); + REQUIRE(reader->loadNumScenes(root.string()) == 2); + const auto dims = reader->loadDimensions(root.string(), 0); + CHECK(dims.sizeZ == 3); + CHECK(dims.physicalSizeZ == 1.0f); + CHECK(reader->loadDimensions(root.string(), 1).sizeZ == 1); + LoadSpec spec; + spec.filepath = root.string(); + LoadProgress progress; + auto image = dynamic_cast(reader.get())->loadVolumeBlocking(spec, progress); + REQUIRE(image != nullptr); + REQUIRE(image->sizeZ() == 3); + const auto first = reinterpret_cast(image->ptr(0, 0))[0]; + const auto second = reinterpret_cast(image->ptr(0, 1))[0]; + const auto third = reinterpret_cast(image->ptr(0, 2))[0]; + CHECK(first < second); + CHECK(second < third); + CHECK(progress.progress() == 1.0f); + std::filesystem::remove_all(root); +}