Merge pull request #239 from inviwo/feature/dicom-fix

MedVis: DICOM sanity check fix

medvis/dicom/include/inviwo/dicom/datastructures/dicomdirtypes.h

@@ -102,9 +102,9 @@ struct IVW_MODULE_DICOM_API Series {
      *   * pixel format
      *   * color space (i.e. Photometric Interpretation) containing  no. of samples per pixel
      *   * intercept and slope
-     *   * TODO: pixel spacing
-     *   * TODO: image orientation (patient) with summed squared difference < 1e-4
-     *   * TODO: image position (patient) with summed squared difference < 1e-4
+     *   * pixel spacing
+     *   * image orientation (patient) with summed squared difference < 1e-4
+     *   * image position (patient) with summed squared difference < 1e-4
      *
      * @param series     series containing images with valid paths
      * @param dicompath  used only in case of an error and should refer to the main DICOM dataset
@@ -119,7 +119,6 @@ struct IVW_MODULE_DICOM_API Series {
      * TODO: proper volume sorting with more heuristics according to
      * https://nipy.org/nibabel/dicom/spm_dicom.html#spm-volume-sorting
      */
-
     void sortImages();
 
     bool empty() const;

medvis/dicom/src/datastructures/dicomdirtypes.cpp

@@ -33,6 +33,7 @@
 #include <inviwo/core/io/datareaderexception.h>
 #include <inviwo/core/util/stdextensions.h>
 #include <inviwo/core/util/filesystem.h>
+#include <inviwo/core/util/glmmat.h>
 
 #include <fmt/format.h>
 #include <fmt/std.h>
@@ -50,6 +51,7 @@
 #include <warn/pop>
 
 #include <algorithm>
+#include <array>
 
 namespace inviwo {
 
@@ -147,14 +149,14 @@ struct ImageMetaData {
         slope = interceptSlopeVec[1];
     }
 
-    size2_t dims;
-    gdcm::PixelFormat pixelformat;
-    gdcm::PhotometricInterpretation photometric;
-    double intercept;
-    double slope;
-    dvec3 pixelSpacing;
-    dvec3 orientation[2];
-    dvec3 origin;
+    size2_t dims{0};
+    gdcm::PixelFormat pixelformat = gdcm::PixelFormat::UNKNOWN;
+    gdcm::PhotometricInterpretation photometric = gdcm::PhotometricInterpretation::UNKNOWN;
+    double intercept{0.0};
+    double slope{1.0};
+    dvec3 pixelSpacing{0.0};
+    std::array<dvec3, 2> orientation = {dvec3{0.0}, dvec3{0.0}};
+    dvec3 origin{0.0};
 };
 
 }  // namespace
@@ -167,16 +169,19 @@ void Series::updateImageInformation(const std::filesystem::path& dicompath) {
     // series:
     //   * dimensions
     //   * pixel format
-    //   * color space (i.e. Photometric Interpretation) containing  no. of samples per
-    //   pixel
+    //   * color space (i.e. Photometric Interpretation) containing no. of samples per
+    //     pixel
     //   * intercept and slope
     //   * pixel spacing
     //   * image orientation (patient) with summed squared difference < 1e-4
-    //   * image position (patient) with summed squared difference < 1e-4
+    //   * image position (patient) with summed squared difference < 1e-4 execpt along axis
+    //     orthogonal to image orientation
+
     bool warnSlopeIntercept = false;
     bool warnPixelSpacing = false;
     bool warnOrientation = false;
     bool warnOrigin = false;
+    auto squaredSum = [](auto vec) { return glm::dot(vec, vec); };
     auto sanityCheck = [&](const ImageMetaData& ref, const ImageMetaData& img) {
         if (ref.dims != img.dims) {
             throw DataReaderException(
@@ -197,22 +202,20 @@ void Series::updateImageInformation(const std::filesystem::path& dicompath) {
         }
         const double dicomDelta = 1.0e-4;
         if (std::abs(ref.slope - img.slope) > dicomDelta ||
-            std::abs(ref.intercept - img.intercept) < dicomDelta) {
+            std::abs(ref.intercept - img.intercept) > dicomDelta) {
             warnSlopeIntercept = true;
         }
 
-        auto ssd = [](auto vec) {  // compute squared sums
-            return glm::dot(vec, vec);
-        };
-
-        if (ssd(ref.pixelSpacing - img.pixelSpacing) > dicomDelta) {
+        if (squaredSum(ref.pixelSpacing - img.pixelSpacing) > dicomDelta) {
             warnPixelSpacing = true;
         }
-        if ((ssd(ref.orientation[0] - img.orientation[0]) > dicomDelta) ||
-            (ssd(ref.orientation[1] - img.orientation[1]) > dicomDelta)) {
+        if ((squaredSum(ref.orientation[0] - img.orientation[0]) > dicomDelta) ||
+            (squaredSum(ref.orientation[1] - img.orientation[1]) > dicomDelta)) {
             warnOrientation = true;
         }
-        if (ssd(ref.origin - img.origin) > dicomDelta) {
+        const dmat3 refT{ref.orientation[0], ref.orientation[1], dvec3{0.0}};
+        const dmat3 imgT{img.orientation[0], img.orientation[1], dvec3{0.0}};
+        if (squaredSum(refT * ref.origin - imgT * img.origin) > dicomDelta) {
             warnOrigin = true;
         }
     };

medvis/dicom/src/io/gdcmvolumereader.cpp

@@ -664,7 +664,7 @@ std::shared_ptr<Volume> GdcmVolumeReader::generateVolume(const gdcm::Image& imag
     LogInfo("spacing: " << spacestr);
     LogInfo("volume size: " << size);
     LogInfo("voxel size: " << voxelsz << "(components: " << (format->getComponents())
-                           << ", component size: " << (format->getSize()) << ")");
+                           << ", component size: " << (format->getSizeInBytes()) << ")");
     LogInfo("sample value range is [" << pixelformat.GetMin() << ", " << pixelformat.GetMax()
                                       << "].");
     pixelformat.Print(gdcm::Trace::GetDebugStream());

tensorvis/tensorvisbase/src/datastructures/tensorfield3d.cpp

@@ -677,9 +677,6 @@ void TensorField3D::computeEigenValuesAndEigenVectors() {
         if (tensor == dmat3(0.0)) {
             return {{std::make_pair(0, dvec3{0}), std::make_pair(0, dvec3{0}),
                      std::make_pair(0, dvec3{0})}};
-            return std::array<std::pair<double, dvec3>, 3>{std::pair<double, dvec3>{0, dvec3(0)},
-                                                           std::pair<double, dvec3>{0, dvec3(0)},
-                                                           std::pair<double, dvec3>{0, dvec3(0)}};
         }
 
         Eigen::EigenSolver<Eigen::Matrix3d> solver(util::glm2eigen(tensor));

tensorvis/tensorvisbase/src/processors/tensorfield2danisotropy.cpp

@@ -104,13 +104,12 @@ void TensorField2DAnisotropy::process() {
             default:
                 throw Exception(IVW_CONTEXT, "Unsupported tensorutil::Anisotropy");
         }
-        return nullptr;
     }();
 
     auto [min, max] = util::layerMinMax(layerRam.get(), IgnoreSpecialValues::No);
     dvec2 range{glm::compMin(min), glm::compMax(max)};
 
-    std::shared_ptr<Layer> layer = std::make_shared<Layer>(layerRam);
+    auto layer = std::make_shared<Layer>(layerRam);
     layer->dataMap.dataRange = range;
     layer->dataMap.valueRange = range;
 

tensorvis/tensorvisio/src/processors/tensorfield2dimport.cpp

@@ -71,12 +71,12 @@ void TensorField2DImport::process() {
         return;
     }
 
-    size_t version;
-    size2_t dimensions;
+    size_t version{0};
+    size2_t dimensions{1};
     auto extents = dvec2(1.0);
-    size_t rank;
-    size_t dimensionality;
-    bool hasEigenInfo;
+    size_t rank{0};
+    size_t dimensionality{0};
+    bool hasEigenInfo = false;
 
     std::string versionStr;
     size_t size;

tensorvis/tensorvisio/src/tensorvisiomodule.cpp

@@ -80,7 +80,6 @@ bool TensorVisIOModule::Converter::convert(TxElement* root) {
         default:
             return false;  // No changes
     }
-    return true;
 }
 
 }  // namespace inviwo