Mercurial > hg > orthanc
diff OrthancServer/SliceOrdering.cpp @ 2804:d88970f1ffbf
fix ordering of non-parallel slices + /tools/reconstruct
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Wed, 29 Aug 2018 13:24:28 +0200 |
parents | 878b59270859 |
children | bbfd95a0c429 |
line wrap: on
line diff
--- a/OrthancServer/SliceOrdering.cpp Tue Aug 28 15:14:33 2018 +0200 +++ b/OrthancServer/SliceOrdering.cpp Wed Aug 29 13:24:28 2018 +0200 @@ -95,12 +95,69 @@ } + static bool IsCloseToZero(double x) + { + return fabs(x) < 10.0 * std::numeric_limits<float>::epsilon(); + } + + + bool SliceOrdering::ComputeNormal(Vector& normal, + const DicomMap& dicom) + { + std::vector<float> cosines; + + if (TokenizeVector(cosines, dicom, DICOM_TAG_IMAGE_ORIENTATION_PATIENT, 6)) + { + assert(cosines.size() == 6); + normal[0] = cosines[1] * cosines[5] - cosines[2] * cosines[4]; + normal[1] = cosines[2] * cosines[3] - cosines[0] * cosines[5]; + normal[2] = cosines[0] * cosines[4] - cosines[1] * cosines[3]; + return true; + } + else + { + return false; + } + } + + + bool SliceOrdering::IsParallelOrOpposite(const Vector& u, + const Vector& v) + { + // Check out "GeometryToolbox::IsParallelOrOpposite()" in Stone of + // Orthanc for explanations + const double u1 = u[0]; + const double u2 = u[1]; + const double u3 = u[2]; + const double normU = sqrt(u1 * u1 + u2 * u2 + u3 * u3); + + const double v1 = v[0]; + const double v2 = v[1]; + const double v3 = v[2]; + const double normV = sqrt(v1 * v1 + v2 * v2 + v3 * v3); + + if (IsCloseToZero(normU * normV)) + { + return false; + } + else + { + const double cosAngle = (u1 * v1 + u2 * v2 + u3 * v3) / (normU * normV); + + return (IsCloseToZero(cosAngle - 1.0) || // Close to +1: Parallel, non-opposite + IsCloseToZero(fabs(cosAngle) - 1.0)); // Close to -1: Parallel, opposite + } + } + + struct SliceOrdering::Instance : public boost::noncopyable { private: std::string instanceId_; bool hasPosition_; Vector position_; + bool hasNormal_; + Vector normal_; bool hasIndexInSeries_; size_t indexInSeries_; unsigned int framesCount_; @@ -141,6 +198,8 @@ position_[2] = tmp[2]; } + hasNormal_ = ComputeNormal(normal_, instance); + std::string s; hasIndexInSeries_ = false; @@ -190,6 +249,17 @@ { return framesCount_; } + + bool HasNormal() const + { + return hasNormal_; + } + + const Vector& GetNormal() const + { + assert(hasNormal_); + return normal_; + } }; @@ -226,15 +296,7 @@ throw OrthancException(ErrorCode_UnknownResource); } - std::vector<float> cosines; - hasNormal_ = TokenizeVector(cosines, series, DICOM_TAG_IMAGE_ORIENTATION_PATIENT, 6); - - if (hasNormal_) - { - normal_[0] = cosines[1] * cosines[5] - cosines[2] * cosines[4]; - normal_[1] = cosines[2] * cosines[3] - cosines[0] * cosines[5]; - normal_[2] = cosines[0] * cosines[4] - cosines[1] * cosines[3]; - } + hasNormal_ = ComputeNormal(normal_, series); } @@ -268,7 +330,10 @@ for (size_t i = 0; i < instances_.size(); i++) { assert(instances_[i] != NULL); - if (!instances_[i]->HasPosition()) + + if (!instances_[i]->HasPosition() || + (instances_[i]->HasNormal() && + !IsParallelOrOpposite(instances_[i]->GetNormal(), normal_))) { return false; }