# HG changeset patch # User Sebastien Jodogne # Date 1516282953 -3600 # Node ID 2115530d37032c1e25a34d1073b67d9a364d413e # Parent 22628d37ef5cfa0625b7351a4f79e54b3e2afcce OrientedBoundingBox diff -r 22628d37ef5c -r 2115530d3703 Framework/Toolbox/GeometryToolbox.cpp --- a/Framework/Toolbox/GeometryToolbox.cpp Tue Jan 16 17:44:16 2018 +0100 +++ b/Framework/Toolbox/GeometryToolbox.cpp Thu Jan 18 14:42:33 2018 +0100 @@ -259,40 +259,40 @@ // (2005). This is a direct, non-optimized translation of Algorithm // 2 in the paper. - static uint8_t tab1[16] = { 255 /* none */, - 0, - 0, - 1, - 1, - 255 /* na */, - 0, - 2, - 2, - 0, - 255 /* na */, - 1, - 1, - 0, - 0, - 255 /* none */ }; + static const uint8_t tab1[16] = { 255 /* none */, + 0, + 0, + 1, + 1, + 255 /* na */, + 0, + 2, + 2, + 0, + 255 /* na */, + 1, + 1, + 0, + 0, + 255 /* none */ }; - static uint8_t tab2[16] = { 255 /* none */, - 3, - 1, - 3, - 2, - 255 /* na */, - 2, - 3, - 3, - 2, - 255 /* na */, - 2, - 3, - 1, - 3, - 255 /* none */ }; + static const uint8_t tab2[16] = { 255 /* none */, + 3, + 1, + 3, + 2, + 255 /* na */, + 2, + 3, + 3, + 2, + 255 /* na */, + 2, + 3, + 1, + 3, + 255 /* none */ }; // Create the coordinates of the rectangle Vector x[4]; @@ -378,5 +378,59 @@ spacingY = 1; } } + + + Matrix CreateRotationMatrixAlongX(double a) + { + // Rotate along X axis (R_x) + // https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations + Matrix r(3, 3); + r(0,0) = 1; + r(0,1) = 0; + r(0,2) = 0; + r(1,0) = 0; + r(1,1) = cos(a); + r(1,2) = -sin(a); + r(2,0) = 0; + r(2,1) = sin(a); + r(2,2) = cos(a); + return r; + } + + + Matrix CreateRotationMatrixAlongY(double a) + { + // Rotate along Y axis (R_y) + // https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations + Matrix r(3, 3); + r(0,0) = cos(a); + r(0,1) = 0; + r(0,2) = sin(a); + r(1,0) = 0; + r(1,1) = 1; + r(1,2) = 0; + r(2,0) = -sin(a); + r(2,1) = 0; + r(2,2) = cos(a); + return r; + } + + + Matrix CreateRotationMatrixAlongZ(double a) + { + // Rotate along Z axis (R_z) + // https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations + Matrix r(3, 3); + r(0,0) = cos(a); + r(0,1) = -sin(a); + r(0,2) = 0; + r(1,0) = sin(a); + r(1,1) = cos(a); + r(1,2) = 0; + r(2,0) = 0; + r(2,1) = 0; + r(2,2) = 1; + return r; + } } } diff -r 22628d37ef5c -r 2115530d3703 Framework/Toolbox/GeometryToolbox.h --- a/Framework/Toolbox/GeometryToolbox.h Tue Jan 16 17:44:16 2018 +0100 +++ b/Framework/Toolbox/GeometryToolbox.h Thu Jan 18 14:42:33 2018 +0100 @@ -28,12 +28,14 @@ # include #endif +#include #include #include namespace OrthancStone { + typedef boost::numeric::ublas::matrix Matrix; typedef boost::numeric::ublas::vector Vector; namespace GeometryToolbox @@ -118,5 +120,11 @@ { return boost::numeric::ublas::inner_prod(point, normal); } + + Matrix CreateRotationMatrixAlongX(double a); + + Matrix CreateRotationMatrixAlongY(double a); + + Matrix CreateRotationMatrixAlongZ(double a); }; } diff -r 22628d37ef5c -r 2115530d3703 Framework/Toolbox/OrientedBoundingBox.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/OrientedBoundingBox.cpp Thu Jan 18 14:42:33 2018 +0100 @@ -0,0 +1,267 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2018 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see . + **/ + + +#include "OrientedBoundingBox.h" + +#include + +#include + +namespace OrthancStone +{ + static bool IntersectPlaneAndSegment(Vector& p, + const Vector& normal, + double d, + const Vector& edgeFrom, + const Vector& edgeTo) + { + // http://geomalgorithms.com/a05-_intersect-1.html#Line-Plane-Intersection + + // Check for parallel line and plane + Vector direction = edgeTo - edgeFrom; + double denominator = boost::numeric::ublas::inner_prod(direction, normal); + + if (fabs(denominator) < 100.0 * std::numeric_limits::epsilon()) + { + return false; + } + else + { + // Compute intersection + double t = -(normal[0] * edgeFrom[0] + + normal[1] * edgeFrom[1] + + normal[2] * edgeFrom[2] + d) / denominator; + + if (t >= 0 && t <= 1) + { + // The intersection lies inside edge segment + p = edgeFrom + t * direction; + return true; + } + else + { + return false; + } + } + } + + + OrientedBoundingBox::OrientedBoundingBox(const ImageBuffer3D& image) + { + unsigned int n = image.GetDepth(); + if (n < 1) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageSize); + } + + const CoordinateSystem3D& geometry = image.GetAxialGeometry(); + Vector dim = image.GetVoxelDimensions(VolumeProjection_Axial); + + u_ = geometry.GetAxisX(); + v_ = geometry.GetAxisY(); + w_ = geometry.GetNormal(); + + hu_ = static_cast(image.GetWidth() * dim[0] / 2.0); + hv_ = static_cast(image.GetHeight() * dim[1] / 2.0); + hw_ = static_cast(image.GetDepth() * dim[2] / 2.0); + + c_ = (geometry.GetOrigin() + + (hu_ - dim[0] / 2.0) * u_ + + (hv_ - dim[1] / 2.0) * v_ + + (hw_ - dim[2] / 2.0) * w_); + } + + + bool OrientedBoundingBox::HasIntersectionWithPlane(std::vector& points, + const Vector& normal, + double d) const + { + assert(normal.size() == 3); + + double r = (hu_ * fabs(boost::numeric::ublas::inner_prod(normal, u_)) + + hv_ * fabs(boost::numeric::ublas::inner_prod(normal, v_)) + + hw_ * fabs(boost::numeric::ublas::inner_prod(normal, w_))); + + double s = boost::numeric::ublas::inner_prod(normal, c_) + d; + + if (fabs(s) >= r) + { + // No intersection, or intersection is reduced to a single point + return false; + } + else + { + Vector p; + + // Loop over all the 12 edges (segments) of the oriented + // bounding box, and check whether they intersect the plane + + // X-aligned edges + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ - v_ * hv_ - w_ * hw_, + c_ + u_ * hu_ - v_ * hv_ - w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ + v_ * hv_ - w_ * hw_, + c_ + u_ * hu_ + v_ * hv_ - w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ - v_ * hv_ + w_ * hw_, + c_ + u_ * hu_ - v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ + v_ * hv_ + w_ * hw_, + c_ + u_ * hu_ + v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + // Y-aligned edges + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ - v_ * hv_ - w_ * hw_, + c_ - u_ * hu_ + v_ * hv_ - w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ + u_ * hu_ - v_ * hv_ - w_ * hw_, + c_ + u_ * hu_ + v_ * hv_ - w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ - v_ * hv_ + w_ * hw_, + c_ - u_ * hu_ + v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ + u_ * hu_ - v_ * hv_ + w_ * hw_, + c_ + u_ * hu_ + v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + // Z-aligned edges + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ - v_ * hv_ - w_ * hw_, + c_ - u_ * hu_ - v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ + u_ * hu_ - v_ * hv_ - w_ * hw_, + c_ + u_ * hu_ - v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ - u_ * hu_ + v_ * hv_ - w_ * hw_, + c_ - u_ * hu_ + v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + if (IntersectPlaneAndSegment(p, normal, d, + c_ + u_ * hu_ + v_ * hv_ - w_ * hw_, + c_ + u_ * hu_ + v_ * hv_ + w_ * hw_)) + { + points.push_back(p); + } + + return true; + } + } + + + bool OrientedBoundingBox::HasIntersection(std::vector& points, + const CoordinateSystem3D& plane) const + { + // From the vector equation of a 3D plane (specified by origin + // and normal), to the general equation of a 3D plane (which + // looses information about the origin of the coordinate system) + const Vector& normal = plane.GetNormal(); + const Vector& origin = plane.GetOrigin(); + double d = -(normal[0] * origin[0] + normal[1] * origin[1] + normal[2] * origin[2]); + + return HasIntersectionWithPlane(points, normal, d); + } + + + bool OrientedBoundingBox::Contains(const Vector& p) const + { + assert(p.size() == 3); + + const Vector q = p - c_; + + return (fabs(boost::numeric::ublas::inner_prod(q, u_)) <= hu_ && + fabs(boost::numeric::ublas::inner_prod(q, v_)) <= hv_ && + fabs(boost::numeric::ublas::inner_prod(q, w_)) <= hw_); + } + + + void OrientedBoundingBox::FromInternalCoordinates(Vector& target, + double x, + double y, + double z) const + { + target = (c_ + + u_ * 2.0 * hu_ * (x - 0.5) + + v_ * 2.0 * hv_ * (y - 0.5) + + w_ * 2.0 * hw_ * (z - 0.5)); + } + + + void OrientedBoundingBox::FromInternalCoordinates(Vector& target, + const Vector& source) const + { + assert(source.size() == 3); + FromInternalCoordinates(target, source[0], source[1], source[2]); + } + + + void OrientedBoundingBox::ToInternalCoordinates(Vector& target, + const Vector& source) const + { + assert(source.size() == 3); + const Vector q = source - c_; + + double x = boost::numeric::ublas::inner_prod(q, u_) / (2.0 * hu_) + 0.5; + double y = boost::numeric::ublas::inner_prod(q, v_) / (2.0 * hv_) + 0.5; + double z = boost::numeric::ublas::inner_prod(q, w_) / (2.0 * hw_) + 0.5; + + GeometryToolbox::AssignVector(target, x, y, z); + } +} diff -r 22628d37ef5c -r 2115530d3703 Framework/Toolbox/OrientedBoundingBox.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Framework/Toolbox/OrientedBoundingBox.h Thu Jan 18 14:42:33 2018 +0100 @@ -0,0 +1,69 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2018 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see . + **/ + + +#pragma once + +#include "GeometryToolbox.h" +#include "../Volumes/ImageBuffer3D.h" + +namespace OrthancStone +{ + class OrientedBoundingBox : public boost::noncopyable + { + private: + Vector c_; // center + Vector u_; // normalized width vector + Vector v_; // normalized height vector + Vector w_; // normalized depth vector + double hu_; // half width + double hv_; // half height + double hw_; // half depth + + public: + OrientedBoundingBox(const ImageBuffer3D& image); + + bool HasIntersectionWithPlane(std::vector& points, + const Vector& normal, + double d) const; + + bool HasIntersection(std::vector& points, + const CoordinateSystem3D& plane) const; + + bool Contains(const Vector& p) const; + + void FromInternalCoordinates(Vector& target, + double x, + double y, + double z) const; + + void FromInternalCoordinates(Vector& target, + const Vector& source) const; + + void ToInternalCoordinates(Vector& target, + const Vector& source) const; + + const Vector& GetCenter() const + { + return c_; + } + }; +} + diff -r 22628d37ef5c -r 2115530d3703 Resources/CMake/OrthancStoneConfiguration.cmake --- a/Resources/CMake/OrthancStoneConfiguration.cmake Tue Jan 16 17:44:16 2018 +0100 +++ b/Resources/CMake/OrthancStoneConfiguration.cmake Thu Jan 18 14:42:33 2018 +0100 @@ -182,6 +182,7 @@ ${ORTHANC_STONE_DIR}/Framework/Toolbox/Extent2D.cpp ${ORTHANC_STONE_DIR}/Framework/Toolbox/GeometryToolbox.cpp ${ORTHANC_STONE_DIR}/Framework/Toolbox/MessagingToolbox.cpp + ${ORTHANC_STONE_DIR}/Framework/Toolbox/OrientedBoundingBox.cpp ${ORTHANC_STONE_DIR}/Framework/Toolbox/OrthancSlicesLoader.cpp ${ORTHANC_STONE_DIR}/Framework/Toolbox/ParallelSlices.cpp ${ORTHANC_STONE_DIR}/Framework/Toolbox/ParallelSlicesCursor.cpp