# HG changeset patch # User Sebastien Jodogne # Date 1620842270 -7200 # Node ID fca942f4b4a7a49ab3f30304312b133106443787 # Parent 95ece40bb29845536165cbece66597e27dee8829 fix conversion from voxel centers to texture borders diff -r 95ece40bb298 -r fca942f4b4a7 OrthancStone/Resources/Documentation/Coordinates.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/OrthancStone/Resources/Documentation/Coordinates.txt Wed May 12 19:57:50 2021 +0200 @@ -0,0 +1,24 @@ + +Some notes about the coordinate systems +======================================= + +* The 3D coordinates are expressed in the PATIENT system, *not* the + "gantry" coordinates. As a consequence, the tags "Image Position + Patient" (0020,0032) and "Image Orientation Patient" (0020,0037) + describe the 3D geometry of a slice. + +* The tag "Patient Position" (0018,5100) could be used to convert from + patient coordinates to gantry coordinates, but this would be useful + to device manufacturers, whereas Stone primarly deals with users of + modalities. + +* The "Image Position Patient" gives the 3D coordinates of the CENTER + of the voxel of slice. + +* In 2D compositors, the origin of a texture corresponds to the CORNER + of the texture (*not* to the center of the first pixel). Roughly + speaking, the operation "floor()" must be applied to move from + floating coordinates to pixel coordinates. + +* The classes deriving from "IVolumeSlicer" must pay to attention to + convert from centers of 3D voxels to bounding boxes of 2D pixels. diff -r 95ece40bb298 -r fca942f4b4a7 OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.cpp --- a/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.cpp Wed May 12 17:43:51 2021 +0200 +++ b/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.cpp Wed May 12 19:57:50 2021 +0200 @@ -165,10 +165,17 @@ const Vector& pixelOffsetX, const Vector& pixelOffsetY) { + /** + * Shift from the center of the voxel (DICOM convention for 3D + * slices) to the corner of the voxel, because 2D textures are + * expressed relatively to their borders. (*) + **/ + Vector p = origin + cuttingPlane.GetOrigin() - 0.5 * pixelOffsetX - 0.5 * pixelOffsetY; + double x0, y0, x1, y1, x2, y2; - cuttingPlane.ProjectPoint(x0, y0, origin + cuttingPlane.GetOrigin()); - cuttingPlane.ProjectPoint(x1, y1, origin + cuttingPlane.GetOrigin() + pixelOffsetX); - cuttingPlane.ProjectPoint(x2, y2, origin + cuttingPlane.GetOrigin() + pixelOffsetY); + cuttingPlane.ProjectPoint(x0, y0, p); + cuttingPlane.ProjectPoint(x1, y1, p + pixelOffsetX); + cuttingPlane.ProjectPoint(x2, y2, p + pixelOffsetY); /** @@ -215,7 +222,7 @@ AffineTransform2D::CreateOffset(originX_, originY_), AffineTransform2D::CreateRotation(angle_), AffineTransform2D::CreateScaling(pixelSpacingX_, pixelSpacingY_), - AffineTransform2D::CreateOffset(-0.5, -0.5), + AffineTransform2D::CreateOffset(-0.5, -0.5), // (*) AffineTransform2D::CreateFlip(flipX_, flipY_, width, height)); } else diff -r 95ece40bb298 -r fca942f4b4a7 OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.h --- a/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.h Wed May 12 17:43:51 2021 +0200 +++ b/OrthancStone/Sources/Scene2D/TextureBaseSceneLayer.h Wed May 12 19:57:50 2021 +0200 @@ -59,6 +59,12 @@ public: TextureBaseSceneLayer(); + private: + /** + * TODO - The methods below could be removed, as well as the + * corresponding members + **/ + // Center of the top-left pixel void SetOrigin(double x, double y); @@ -69,8 +75,6 @@ // In radians void SetAngle(double angle); - void SetLinearInterpolation(bool isLinearInterpolation); - void SetFlipX(bool flip); void SetFlipY(bool flip); @@ -100,16 +104,6 @@ return angle_; } - bool IsLinearInterpolation() const - { - return isLinearInterpolation_; - } - - bool HasTexture() const - { - return (texture_.get() != NULL); - } - bool IsFlipX() const { return flipX_; @@ -120,15 +114,32 @@ return flipY_; } + public: + bool IsLinearInterpolation() const + { + return isLinearInterpolation_; + } + + void SetLinearInterpolation(bool isLinearInterpolation); + + bool HasTexture() const + { + return (texture_.get() != NULL); + } + const Orthanc::ImageAccessor& GetTexture() const; void SetTransform(const AffineTransform2D& transform); void ClearTransform(); - // Initialize a transform that maps a texture slice in 3D, to a - // cutting plane (the cutting plane should be parallel to the 3D - // slice). The "pixelOffsetX/Y" must take pixel spacing into account. + /** + * Initialize a transform that maps a texture slice in 3D, to a + * cutting plane (the cutting plane should be parallel to the 3D + * slice). The "pixelOffsetX/Y" must take pixel spacing into + * account. This method automatically converts from voxel centers + * (3D) to pixel corners (2D). + **/ void SetCuttingPlaneTransform(const CoordinateSystem3D& cuttingPlane, const Vector& origin, // coordinates of the center of the voxel const Vector& pixelOffsetX, // 3D offset from (0,0) voxel to (1,0) voxel diff -r 95ece40bb298 -r fca942f4b4a7 OrthancStone/Sources/Volumes/DicomVolumeImageMPRSlicer.cpp --- a/OrthancStone/Sources/Volumes/DicomVolumeImageMPRSlicer.cpp Wed May 12 17:43:51 2021 +0200 +++ b/OrthancStone/Sources/Volumes/DicomVolumeImageMPRSlicer.cpp Wed May 12 19:57:50 2021 +0200 @@ -95,11 +95,6 @@ const CoordinateSystem3D& system = volume_.GetGeometry().GetProjectionGeometry(projection_); - /** - * TODO => There was a shift of (0.5, 0.5) introduced by - * TextureBaseSceneLayer::GetTransform(). Is it an error? - **/ - Vector pixelSpacing = volume_.GetGeometry().GetVoxelDimensions(projection_); texture->SetCuttingPlaneTransform(cuttingPlane, system.GetOrigin(), diff -r 95ece40bb298 -r fca942f4b4a7 OrthancStone/Sources/Volumes/DicomVolumeImageReslicer.cpp --- a/OrthancStone/Sources/Volumes/DicomVolumeImageReslicer.cpp Wed May 12 17:43:51 2021 +0200 +++ b/OrthancStone/Sources/Volumes/DicomVolumeImageReslicer.cpp Wed May 12 19:57:50 2021 +0200 @@ -87,11 +87,14 @@ const Vector p1 = cuttingPlane.MapSliceToWorldCoordinates(x1, y1); const Vector p2 = cuttingPlane.MapSliceToWorldCoordinates(x1, y2); + // The "0.5" shift is to move from the corner of voxel to the center of the voxel + if (1) { - texture->SetCuttingPlaneTransform(cuttingPlane, p1, - s * cuttingPlane.GetAxisX(), - s * cuttingPlane.GetAxisY()); + texture->SetCuttingPlaneTransform( + cuttingPlane, p1 + 0.5 * s * cuttingPlane.GetAxisX() + 0.5 * s * cuttingPlane.GetAxisY(), + s * cuttingPlane.GetAxisX(), + s * cuttingPlane.GetAxisY()); } else { @@ -100,9 +103,10 @@ * possible for the X axis? **/ - texture->SetCuttingPlaneTransform(cuttingPlane, p2, - s * cuttingPlane.GetAxisX(), - -s * cuttingPlane.GetAxisY()); + texture->SetCuttingPlaneTransform( + cuttingPlane, p2 + 0.5 * s * cuttingPlane.GetAxisX() + 0.5 * s * cuttingPlane.GetAxisY(), + s * cuttingPlane.GetAxisX(), + -s * cuttingPlane.GetAxisY()); } #else diff -r 95ece40bb298 -r fca942f4b4a7 UnitTestsSources/VolumeRenderingTests.cpp --- a/UnitTestsSources/VolumeRenderingTests.cpp Wed May 12 17:43:51 2021 +0200 +++ b/UnitTestsSources/VolumeRenderingTests.cpp Wed May 12 19:57:50 2021 +0200 @@ -21,6 +21,7 @@ #include "../OrthancStone/Sources/Scene2D/CairoCompositor.h" #include "../OrthancStone/Sources/Scene2D/CopyStyleConfigurator.h" +#include "../OrthancStone/Sources/Scene2D/ColorTextureSceneLayer.h" #include "../OrthancStone/Sources/Volumes/DicomVolumeImageMPRSlicer.h" #include "../OrthancStone/Sources/Volumes/DicomVolumeImageReslicer.h" @@ -29,7 +30,7 @@ #include -TEST(VolumeRendering, Basic) +TEST(VolumeRendering, Axial) { Orthanc::DicomMap dicom; dicom.SetValue(Orthanc::DICOM_TAG_STUDY_INSTANCE_UID, "study", false); @@ -43,7 +44,7 @@ OrthancStone::VolumeImageGeometry geometry; geometry.SetSizeInVoxels(3, 3, 1); geometry.SetAxialGeometry(axial); - + boost::shared_ptr volume(new OrthancStone::DicomVolumeImage); volume->Initialize(geometry, Orthanc::PixelFormat_Grayscale8, false); volume->SetDicomParameters(OrthancStone::DicomInstanceParameters(dicom)); @@ -62,6 +63,11 @@ } } + OrthancStone::Vector v = volume->GetGeometry().GetVoxelDimensions(OrthancStone::VolumeProjection_Axial); + ASSERT_FLOAT_EQ(1, v[0]); + ASSERT_FLOAT_EQ(1, v[1]); + ASSERT_FLOAT_EQ(1, v[2]); + OrthancStone::CoordinateSystem3D viewpoint; for (unsigned int mode = 0; mode < 2; mode++) @@ -83,6 +89,15 @@ OrthancStone::CopyStyleConfigurator configurator; std::unique_ptr layer(slice->CreateSceneLayer(&configurator, viewpoint)); + ASSERT_EQ(OrthancStone::ISceneLayer::Type_FloatTexture, layer->GetType()); + + OrthancStone::Extent2D box; + layer->GetBoundingBox(box); + ASSERT_FLOAT_EQ(-1.0f, box.GetX1()); + ASSERT_FLOAT_EQ(-1.0f, box.GetY1()); + ASSERT_FLOAT_EQ(2.0f, box.GetX2()); + ASSERT_FLOAT_EQ(2.0f, box.GetY2()); + { const Orthanc::ImageAccessor& a = dynamic_cast(*layer).GetTexture(); Orthanc::Image i(Orthanc::PixelFormat_Grayscale8, a.GetWidth(), a.GetHeight(), false); @@ -144,4 +159,3 @@ Orthanc::ImageTraits::GetPixel(pixel, j, 4, 4); ASSERT_EQ(200, pixel.red_); } } -