Mercurial > hg > orthanc

comparison OrthancServer/SliceOrdering.cpp @ 2804:d88970f1ffbf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
2803:579acc5e5412 2804:d88970f1ffbf
93 return TokenizeVector(result, value->GetContent(), expectedSize); 93 return TokenizeVector(result, value->GetContent(), expectedSize);
94 } 94 }
95 } 95 }
96 96
97 97
98 static bool IsCloseToZero(double x)
99 {
100 return fabs(x) < 10.0 * std::numeric_limits<float>::epsilon();
101 }
102
103
104 bool SliceOrdering::ComputeNormal(Vector& normal,
105 const DicomMap& dicom)
106 {
107 std::vector<float> cosines;
108
109 if (TokenizeVector(cosines, dicom, DICOM_TAG_IMAGE_ORIENTATION_PATIENT, 6))
110 {
111 assert(cosines.size() == 6);
112 normal[0] = cosines[1] * cosines[5] - cosines[2] * cosines[4];
113 normal[1] = cosines[2] * cosines[3] - cosines[0] * cosines[5];
114 normal[2] = cosines[0] * cosines[4] - cosines[1] * cosines[3];
115 return true;
116 }
117 else
118 {
119 return false;
120 }
121 }
122
123
124 bool SliceOrdering::IsParallelOrOpposite(const Vector& u,
125 const Vector& v)
126 {
127 // Check out "GeometryToolbox::IsParallelOrOpposite()" in Stone of
128 // Orthanc for explanations
129 const double u1 = u[0];
130 const double u2 = u[1];
131 const double u3 = u[2];
132 const double normU = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
133
134 const double v1 = v[0];
135 const double v2 = v[1];
136 const double v3 = v[2];
137 const double normV = sqrt(v1 * v1 + v2 * v2 + v3 * v3);
138
139 if (IsCloseToZero(normU * normV))
140 {
141 return false;
142 }
143 else
144 {
145 const double cosAngle = (u1 * v1 + u2 * v2 + u3 * v3) / (normU * normV);
146
147 return (IsCloseToZero(cosAngle - 1.0) || // Close to +1: Parallel, non-opposite
148 IsCloseToZero(fabs(cosAngle) - 1.0)); // Close to -1: Parallel, opposite
149 }
150 }
151
152
98 struct SliceOrdering::Instance : public boost::noncopyable 153 struct SliceOrdering::Instance : public boost::noncopyable
99 { 154 {
100 private: 155 private:
101 std::string instanceId_; 156 std::string instanceId_;
102 bool hasPosition_; 157 bool hasPosition_;
103 Vector position_; 158 Vector position_;
159 bool hasNormal_;
160 Vector normal_;
104 bool hasIndexInSeries_; 161 bool hasIndexInSeries_;
105 size_t indexInSeries_; 162 size_t indexInSeries_;
106 unsigned int framesCount_; 163 unsigned int framesCount_;
107 164
108 public: 165 public:
139 position_[0] = tmp[0]; 196 position_[0] = tmp[0];
140 position_[1] = tmp[1]; 197 position_[1] = tmp[1];
141 position_[2] = tmp[2]; 198 position_[2] = tmp[2];
142 } 199 }
143 200
201 hasNormal_ = ComputeNormal(normal_, instance);
202
144 std::string s; 203 std::string s;
145 hasIndexInSeries_ = false; 204 hasIndexInSeries_ = false;
146 205
147 try 206 try
148 { 207 {
188 247
189 unsigned int GetFramesCount() const 248 unsigned int GetFramesCount() const
190 { 249 {
191 return framesCount_; 250 return framesCount_;
192 } 251 }
252
253 bool HasNormal() const
254 {
255 return hasNormal_;
256 }
257
258 const Vector& GetNormal() const
259 {
260 assert(hasNormal_);
261 return normal_;
262 }
193 }; 263 };
194 264
195 265
196 class SliceOrdering::PositionComparator 266 class SliceOrdering::PositionComparator
197 { 267 {
224 if (!index_.GetMainDicomTags(series, seriesId_, ResourceType_Series, ResourceType_Series)) 294 if (!index_.GetMainDicomTags(series, seriesId_, ResourceType_Series, ResourceType_Series))
225 { 295 {
226 throw OrthancException(ErrorCode_UnknownResource); 296 throw OrthancException(ErrorCode_UnknownResource);
227 } 297 }
228 298
229 std::vector<float> cosines; 299 hasNormal_ = ComputeNormal(normal_, series);
230 hasNormal_ = TokenizeVector(cosines, series, DICOM_TAG_IMAGE_ORIENTATION_PATIENT, 6);
231
232 if (hasNormal_)
233 {
234 normal_[0] = cosines[1] * cosines[5] - cosines[2] * cosines[4];
235 normal_[1] = cosines[2] * cosines[3] - cosines[0] * cosines[5];
236 normal_[2] = cosines[0] * cosines[4] - cosines[1] * cosines[3];
237 }
238 } 300 }
239 301
240 302
241 void SliceOrdering::CreateInstances() 303 void SliceOrdering::CreateInstances()
242 { 304 {
266 } 328 }
267 329
268 for (size_t i = 0; i < instances_.size(); i++) 330 for (size_t i = 0; i < instances_.size(); i++)
269 { 331 {
270 assert(instances_[i] != NULL); 332 assert(instances_[i] != NULL);
271 if (!instances_[i]->HasPosition()) 333
334 if (!instances_[i]->HasPosition() ||
335 (instances_[i]->HasNormal() &&
336 !IsParallelOrOpposite(instances_[i]->GetNormal(), normal_)))
272 { 337 {
273 return false; 338 return false;
274 } 339 }
275 } 340 }
276 341