# HG changeset patch # User Sebastien Jodogne # Date 1574181005 -3600 # Node ID 19b1c8caade4f49d0bd1823774af7da0d0ea655f # Parent 59906485896f79d5ef87e52cdfee8e1443a6eb96 fix sagittal geometry diff -r 59906485896f -r 19b1c8caade4 Framework/Toolbox/CoordinateSystem3D.cpp --- a/Framework/Toolbox/CoordinateSystem3D.cpp Tue Nov 19 12:52:07 2019 +0100 +++ b/Framework/Toolbox/CoordinateSystem3D.cpp Tue Nov 19 17:30:05 2019 +0100 @@ -241,4 +241,14 @@ return s; } + + CoordinateSystem3D CoordinateSystem3D::NormalizeCuttingPlane(const CoordinateSystem3D& plane) + { + double ox, oy; + plane.ProjectPoint(ox, oy, LinearAlgebra::CreateVector(0, 0, 0)); + + CoordinateSystem3D normalized(plane); + normalized.SetOrigin(plane.MapSliceToWorldCoordinates(ox, oy)); + return normalized; + } } diff -r 59906485896f -r 19b1c8caade4 Framework/Toolbox/CoordinateSystem3D.h --- a/Framework/Toolbox/CoordinateSystem3D.h Tue Nov 19 12:52:07 2019 +0100 +++ b/Framework/Toolbox/CoordinateSystem3D.h Tue Nov 19 17:30:05 2019 +0100 @@ -22,6 +22,7 @@ #pragma once #include "LinearAlgebra.h" +#include "../Scene2D/ScenePoint2D.h" #include @@ -95,12 +96,24 @@ Vector MapSliceToWorldCoordinates(double x, double y) const; + Vector MapSliceToWorldCoordinates(const ScenePoint2D& p) const + { + return MapSliceToWorldCoordinates(p.GetX(), p.GetY()); + } + double ProjectAlongNormal(const Vector& point) const; void ProjectPoint(double& offsetX, double& offsetY, const Vector& point) const; + ScenePoint2D ProjectPoint(const Vector& point) const + { + double x, y; + ProjectPoint(x, y, point); + return ScenePoint2D(x, y); + } + /* Alternated faster implementation (untested yet) */ @@ -120,5 +133,9 @@ static bool ComputeDistance(double& distance, const CoordinateSystem3D& a, const CoordinateSystem3D& b); + + // Normalize a cutting plane so that the origin (0,0,0) of the 3D + // world is mapped to the origin of its (x,y) coordinate system + static CoordinateSystem3D NormalizeCuttingPlane(const CoordinateSystem3D& plane); }; } diff -r 59906485896f -r 19b1c8caade4 Framework/Toolbox/DicomInstanceParameters.cpp --- a/Framework/Toolbox/DicomInstanceParameters.cpp Tue Nov 19 12:52:07 2019 +0100 +++ b/Framework/Toolbox/DicomInstanceParameters.cpp Tue Nov 19 17:30:05 2019 +0100 @@ -470,4 +470,33 @@ return (value * factor + offset); } + + + bool DicomInstanceParameters::Data::ComputeRegularSpacing(double& spacing) const + { + if (frameOffsets_.size() == 0) // Not a RT-DOSE + { + return false; + } + else if (frameOffsets_.size() == 1) + { + spacing = 1; // Edge case: RT-DOSE with one single frame + return true; + } + else + { + spacing = std::abs(frameOffsets_[1] - frameOffsets_[0]); + + for (size_t i = 1; i + 1 < frameOffsets_.size(); i++) + { + double s = frameOffsets_[i + 1] - frameOffsets_[i]; + if (!LinearAlgebra::IsNear(spacing, s, 0.001)) + { + return false; + } + } + + return true; + } + } } diff -r 59906485896f -r 19b1c8caade4 Framework/Toolbox/DicomInstanceParameters.h --- a/Framework/Toolbox/DicomInstanceParameters.h Tue Nov 19 12:52:07 2019 +0100 +++ b/Framework/Toolbox/DicomInstanceParameters.h Tue Nov 19 17:30:05 2019 +0100 @@ -75,6 +75,8 @@ bool useDouble) const; double ApplyRescale(double value) const; + + bool ComputeRegularSpacing(double& target) const; }; @@ -209,6 +211,11 @@ return data_.doseUnits_; } + void SetDoseGridScaling(double value) + { + data_.doseGridScaling_ = value; + } + double GetDoseGridScaling() const { return data_.doseGridScaling_; @@ -218,5 +225,11 @@ { return data_.ApplyRescale(value); } + + // Required for RT-DOSE + bool ComputeRegularSpacing(double& target) const + { + return data_.ComputeRegularSpacing(target); + } }; } diff -r 59906485896f -r 19b1c8caade4 Framework/Volumes/DicomVolumeImageMPRSlicer.cpp --- a/Framework/Volumes/DicomVolumeImageMPRSlicer.cpp Tue Nov 19 12:52:07 2019 +0100 +++ b/Framework/Volumes/DicomVolumeImageMPRSlicer.cpp Tue Nov 19 17:30:05 2019 +0100 @@ -45,8 +45,7 @@ revision_(volume_.GetRevision()) { valid_ = (volume_.HasDicomParameters() && - volume_.GetGeometry().DetectSlice(projection_, sliceIndex_, - cuttingPlane)); + volume_.GetGeometry().DetectSlice(projection_, sliceIndex_, cuttingPlane)); } @@ -82,21 +81,6 @@ ImageBuffer3D::SliceReader reader(volume_.GetPixelData(), projection_, sliceIndex_); texture.reset(dynamic_cast (configurator->CreateTextureFromDicom(reader.GetAccessor(), parameters))); - - // -#if 0 - Orthanc::JpegWriter writer; - writer.SetQuality(60); - static int index = 0; - std::string filePath = "C:\\temp\\sliceReader_P"; - filePath += boost::lexical_cast(projection_); - filePath += "_I"; - filePath += boost::lexical_cast(index); - filePath += ".jpg"; - index++; - writer.WriteToFile(filePath, reader.GetAccessor()); -#endif - // } const CoordinateSystem3D& system = volume_.GetGeometry().GetProjectionGeometry(projection_); @@ -105,53 +89,11 @@ cuttingPlane.ProjectPoint(x0, y0, system.GetOrigin()); cuttingPlane.ProjectPoint(x1, y1, system.GetOrigin() + system.GetAxisX()); - // -#if 0 { - LOG(ERROR) << "+----------------------------------------------------+"; - LOG(ERROR) << "| DicomVolumeImageMPRSlicer::Slice::CreateSceneLayer |"; - LOG(ERROR) << "+----------------------------------------------------+"; - std::string projectionString; - switch (projection_) - { - case VolumeProjection_Coronal: - projectionString = "CORONAL"; - break; - case VolumeProjection_Axial: - projectionString = "CORONAL"; - break; - case VolumeProjection_Sagittal: - projectionString = "SAGITTAL"; - break; - default: - ORTHANC_ASSERT(false); - } - if(volume_.GetGeometry().GetDepth() == 200) - LOG(ERROR) << "| CT IMAGE 512x512 with projection " << projectionString; - else - LOG(ERROR) << "| RTDOSE IMAGE NNNxNNN with projection " << projectionString; - LOG(ERROR) << "+----------------------------------------------------+"; - LOG(ERROR) << "| cuttingPlane = " << cuttingPlane; - LOG(ERROR) << "| point to project = " << system.GetOrigin(); - LOG(ERROR) << "| result = x0: " << x0 << " y0: " << y0; - LOG(ERROR) << "+----------------------- END ------------------------+"; + double xz, yz; + cuttingPlane.ProjectPoint(xz, yz, LinearAlgebra::CreateVector(0, 0, 0)); + texture->SetOrigin(x0 - xz, y0 - yz); } -#endif - // - -#if 1 // BGO 2019-08-13 - // The sagittal coordinate system has a Y vector going down. The displayed - // image (scene coords) has a Y vector pointing upwards (towards the patient - // coord Z index) - // we need to flip the Y local coordinates to get the scene-coord offset. - // TODO: this is quite ugly. Isn't there a better way? - if(projection_ == VolumeProjection_Sagittal) - texture->SetOrigin(x0, -y0); - else - texture->SetOrigin(x0, y0); -#else - texture->SetOrigin(x0, y0); -#endif double dx = x1 - x0; double dy = y1 - y0; @@ -164,19 +106,6 @@ Vector tmp = volume_.GetGeometry().GetVoxelDimensions(projection_); texture->SetPixelSpacing(tmp[0], tmp[1]); - // - { - //using std::endl; - //std::stringstream ss; - //ss << "DicomVolumeImageMPRSlicer::Slice::CreateSceneLayer | cuttingPlane = " << cuttingPlane << " | projection_ = " << projection_ << endl; - //ss << "volume_.GetGeometry().GetProjectionGeometry(projection_) = " << system << endl; - //ss << "cuttingPlane.ProjectPoint(x0, y0, system.GetOrigin()); --> | x0 = " << x0 << " | y0 = " << y0 << "| x1 = " << x1 << " | y1 = " << y1 << endl; - //ss << "volume_.GetGeometry() = " << volume_.GetGeometry() << endl; - //ss << "volume_.GetGeometry() = " << volume_.GetGeometry() << endl; - //LOG(ERROR) << ss.str(); - } - // - return texture.release(); } diff -r 59906485896f -r 19b1c8caade4 Framework/Volumes/VolumeImageGeometry.cpp --- a/Framework/Volumes/VolumeImageGeometry.cpp Tue Nov 19 12:52:07 2019 +0100 +++ b/Framework/Volumes/VolumeImageGeometry.cpp Tue Nov 19 17:30:05 2019 +0100 @@ -39,7 +39,7 @@ sagittalGeometry_ = CoordinateSystem3D(p, axialGeometry_.GetAxisY(), - axialGeometry_.GetNormal()); + -axialGeometry_.GetNormal()); Vector origin = ( axialGeometry_.MapSliceToWorldCoordinates(-0.5 * voxelDimensions_[0],