Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/ImageGeometry.cpp @ 182:2cbfb08f3a95 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ImageGeometry.cpp
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 15 Mar 2018 12:12:01 +0100
parents
children 465b294a55f0
comparison
equal deleted inserted replaced
181:ff8556874557 182:2cbfb08f3a95
1 /**
2 * Stone of Orthanc
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
4 * Department, University Hospital of Liege, Belgium
5 * Copyright (C) 2017-2018 Osimis S.A., Belgium
6 *
7 * This program is free software: you can redistribute it and/or
8 * modify it under the terms of the GNU Affero General Public License
9 * as published by the Free Software Foundation, either version 3 of
10 * the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Affero General Public License for more details.
16 *
17 * You should have received a copy of the GNU Affero General Public License
18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 **/
20
21
22 #include "ImageGeometry.h"
23
24 #include "Extent2D.h"
25 #include "SubpixelReader.h"
26
27 #include <Core/Images/ImageProcessing.h>
28 #include <Core/Logging.h>
29
30
31 namespace OrthancStone
32 {
33 static void AddTransformedPoint(Extent2D& extent,
34 const Matrix& a,
35 double x,
36 double y)
37 {
38 assert(a.size1() == 3 &&
39 a.size2() == 3);
40
41 Vector p = LinearAlgebra::Product(a, LinearAlgebra::CreateVector(x, y, 1));
42
43 if (!LinearAlgebra::IsCloseToZero(p[2]))
44 {
45 extent.AddPoint(p[0] / p[2], p[1] / p[2]);
46 }
47 }
48
49
50 bool GetPerpectiveTransformExtent(unsigned int& x1,
51 unsigned int& y1,
52 unsigned int& x2,
53 unsigned int& y2,
54 const Matrix& a,
55 unsigned int sourceWidth,
56 unsigned int sourceHeight,
57 unsigned int targetWidth,
58 unsigned int targetHeight)
59 {
60 if (targetWidth == 0 ||
61 targetHeight == 0)
62 {
63 return false;
64 }
65
66 Extent2D extent;
67 AddTransformedPoint(extent, a, 0, 0);
68 AddTransformedPoint(extent, a, sourceWidth, 0);
69 AddTransformedPoint(extent, a, 0, sourceHeight);
70 AddTransformedPoint(extent, a, sourceWidth, sourceHeight);
71
72 if (extent.IsEmpty())
73 {
74 return false;
75 }
76 else
77 {
78 int tmp;
79
80 tmp = std::floor(extent.GetX1());
81 if (tmp < 0)
82 {
83 x1 = 0;
84 }
85 else
86 {
87 x1 = static_cast<unsigned int>(tmp);
88 }
89
90 tmp = std::floor(extent.GetY1());
91 if (tmp < 0)
92 {
93 y1 = 0;
94 }
95 else
96 {
97 y1 = static_cast<unsigned int>(tmp);
98 }
99
100 tmp = std::ceil(extent.GetX2());
101 if (tmp < 0)
102 {
103 return false;
104 }
105 else if (static_cast<unsigned int>(tmp) >= targetWidth)
106 {
107 x2 = targetWidth - 1;
108 }
109 else
110 {
111 x2 = static_cast<unsigned int>(tmp);
112 }
113
114 tmp = std::ceil(extent.GetY2());
115 if (tmp < 0)
116 {
117 return false;
118 }
119 else if (static_cast<unsigned int>(tmp) >= targetHeight)
120 {
121 y2 = targetHeight - 1;
122 }
123 else
124 {
125 y2 = static_cast<unsigned int>(tmp);
126 }
127
128 return (x1 <= x2 &&
129 y1 <= y2);
130 }
131 }
132
133
134 template <typename Reader,
135 bool HasOffsetX,
136 bool HasOffsetY>
137 static void ApplyAffineTransformToRow(typename Reader::PixelType* p,
138 Reader& reader,
139 unsigned int x1,
140 unsigned int x2,
141 float positionX,
142 float positionY,
143 float offsetX,
144 float offsetY)
145 {
146 typename Reader::PixelType value;
147
148 for (unsigned int x = x1; x <= x2; x++, p++)
149 {
150 if (reader.GetValue(value, positionX, positionY))
151 {
152 *p = value;
153 }
154 else
155 {
156 Reader::Traits::SetZero(*p);
157 }
158
159 if (HasOffsetX)
160 {
161 positionX += offsetX;
162 }
163
164 if (HasOffsetY)
165 {
166 positionY += offsetY;
167 }
168 }
169 }
170
171
172 template <Orthanc::PixelFormat Format,
173 ImageInterpolation Interpolation>
174 static void ApplyAffineInternal(Orthanc::ImageAccessor& target,
175 const Orthanc::ImageAccessor& source,
176 const Matrix& a)
177 {
178 assert(target.GetFormat() == Format &&
179 source.GetFormat() == Format);
180
181 typedef SubpixelReader<Format, Interpolation> Reader;
182 typedef typename Reader::PixelType PixelType;
183
184 if (Format == Orthanc::PixelFormat_RGB24)
185 {
186 Orthanc::ImageProcessing::Set(target, 0, 0, 0, 255);
187 }
188 else
189 {
190 Orthanc::ImageProcessing::Set(target, 0);
191 }
192
193 Matrix inva;
194 if (!LinearAlgebra::InvertMatrixUnsafe(inva, a))
195 {
196 // Singular matrix
197 return;
198 }
199
200 Reader reader(source);
201
202 unsigned int x1, y1, x2, y2;
203
204 if (GetPerpectiveTransformExtent(x1, y1, x2, y2, a,
205 source.GetWidth(), source.GetHeight(),
206 target.GetWidth(), target.GetHeight()))
207 {
208 const size_t targetPitch = target.GetPitch();
209 uint8_t *targetRow = reinterpret_cast<uint8_t*>(reinterpret_cast<PixelType*>(target.GetRow(y1)) + x1);
210
211 for (unsigned int y = y1; y <= y2; y++)
212 {
213 Vector start;
214 LinearAlgebra::AssignVector(start, static_cast<double>(x1) + 0.5,
215 static_cast<double>(y) + 0.5, 1);
216 start = boost::numeric::ublas::prod(inva, start);
217 assert(LinearAlgebra::IsNear(1.0, start(2)));
218
219 Vector offset;
220 LinearAlgebra::AssignVector(offset, static_cast<double>(x1) + 1.5,
221 static_cast<double>(y) + 0.5, 1);
222 offset = boost::numeric::ublas::prod(inva, offset) - start;
223 assert(LinearAlgebra::IsNear(0.0, offset(2)));
224
225 float startX = static_cast<float>(start[0]);
226 float startY = static_cast<float>(start[1]);
227 float offsetX = static_cast<float>(offset[0]);
228 float offsetY = static_cast<float>(offset[1]);
229
230 PixelType* pixel = reinterpret_cast<PixelType*>(targetRow);
231 if (LinearAlgebra::IsCloseToZero(offsetX))
232 {
233 ApplyAffineTransformToRow<Reader, false, true>
234 (pixel, reader, x1, x2, startX, startY, offsetX, offsetY);
235 }
236 else if (LinearAlgebra::IsCloseToZero(offsetY))
237 {
238 ApplyAffineTransformToRow<Reader, true, false>
239 (pixel, reader, x1, x2, startX, startY, offsetX, offsetY);
240 }
241 else
242 {
243 ApplyAffineTransformToRow<Reader, true, true>
244 (pixel, reader, x1, x2, startX, startY, offsetX, offsetY);
245 }
246
247 targetRow += targetPitch;
248 }
249 }
250 }
251
252
253 void ApplyAffineTransform(Orthanc::ImageAccessor& target,
254 const Orthanc::ImageAccessor& source,
255 double a11,
256 double a12,
257 double b1,
258 double a21,
259 double a22,
260 double b2,
261 ImageInterpolation interpolation)
262 {
263 if (source.GetFormat() != target.GetFormat())
264 {
265 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
266 }
267
268 if (interpolation != ImageInterpolation_Nearest &&
269 interpolation != ImageInterpolation_Bilinear)
270 {
271 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
272 }
273
274 Matrix a;
275 a.resize(3, 3);
276 a(0, 0) = a11;
277 a(0, 1) = a12;
278 a(0, 2) = b1;
279 a(1, 0) = a21;
280 a(1, 1) = a22;
281 a(1, 2) = b2;
282 a(2, 0) = 0;
283 a(2, 1) = 0;
284 a(2, 2) = 1;
285
286 switch (source.GetFormat())
287 {
288 case Orthanc::PixelFormat_Grayscale8:
289 switch (interpolation)
290 {
291 case ImageInterpolation_Nearest:
292 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale8,
293 ImageInterpolation_Nearest>(target, source, a);
294 break;
295
296 case ImageInterpolation_Bilinear:
297 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale8,
298 ImageInterpolation_Bilinear>(target, source, a);
299 break;
300
301 default:
302 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
303 }
304 break;
305
306 case Orthanc::PixelFormat_Grayscale16:
307 switch (interpolation)
308 {
309 case ImageInterpolation_Nearest:
310 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale16,
311 ImageInterpolation_Nearest>(target, source, a);
312 break;
313
314 case ImageInterpolation_Bilinear:
315 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale16,
316 ImageInterpolation_Bilinear>(target, source, a);
317 break;
318
319 default:
320 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
321 }
322 break;
323
324 case Orthanc::PixelFormat_SignedGrayscale16:
325 switch (interpolation)
326 {
327 case ImageInterpolation_Nearest:
328 ApplyAffineInternal<Orthanc::PixelFormat_SignedGrayscale16,
329 ImageInterpolation_Nearest>(target, source, a);
330 break;
331
332 case ImageInterpolation_Bilinear:
333 ApplyAffineInternal<Orthanc::PixelFormat_SignedGrayscale16,
334 ImageInterpolation_Bilinear>(target, source, a);
335 break;
336
337 default:
338 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
339 }
340 break;
341
342 case Orthanc::PixelFormat_RGB24:
343 switch (interpolation)
344 {
345 case ImageInterpolation_Nearest:
346 ApplyAffineInternal<Orthanc::PixelFormat_RGB24,
347 ImageInterpolation_Nearest>(target, source, a);
348 break;
349
350 default:
351 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
352 }
353 break;
354
355 default:
356 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
357 }
358 }
359
360
361 template <Orthanc::PixelFormat Format,
362 ImageInterpolation Interpolation>
363 static void ApplyPerspectiveInternal(Orthanc::ImageAccessor& target,
364 const Orthanc::ImageAccessor& source,
365 const Matrix& a,
366 const Matrix& inva)
367 {
368 assert(target.GetFormat() == Format &&
369 source.GetFormat() == Format);
370
371 typedef SubpixelReader<Format, Interpolation> Reader;
372 typedef typename Reader::PixelType PixelType;
373
374 Reader reader(source);
375 unsigned int x1, y1, x2, y2;
376
377 const float floatWidth = source.GetWidth();
378 const float floatHeight = source.GetHeight();
379
380 if (GetPerpectiveTransformExtent(x1, y1, x2, y2, a,
381 source.GetWidth(), source.GetHeight(),
382 target.GetWidth(), target.GetHeight()))
383 {
384 const size_t targetPitch = target.GetPitch();
385 uint8_t *targetRow = reinterpret_cast<uint8_t*>(reinterpret_cast<PixelType*>(target.GetRow(y1)) + x1);
386
387 for (unsigned int y = y1; y <= y2; y++)
388 {
389 PixelType *p = reinterpret_cast<PixelType*>(targetRow);
390
391 for (unsigned int x = x1; x <= x2; x++)
392 {
393 Vector v;
394 LinearAlgebra::AssignVector(v, static_cast<double>(x) + 0.5,
395 static_cast<double>(y) + 0.5, 1);
396
397 Vector vv = LinearAlgebra::Product(inva, v);
398
399 assert(!LinearAlgebra::IsCloseToZero(vv[2]));
400 const double w = 1.0 / vv[2];
401 const float sourceX = static_cast<float>(vv[0] * w);
402 const float sourceY = static_cast<float>(vv[1] * w);
403
404 // Make sure no integer overflow will occur after truncation
405 // (the static_cast<unsigned int> could otherwise throw an
406 // exception in WebAssembly if strong perspective effects)
407 if (sourceX < floatWidth &&
408 sourceY < floatHeight)
409 {
410 reader.GetValue(*p, sourceX, sourceY);
411 }
412 else
413 {
414 Reader::Traits::SetZero(*p);
415 }
416
417 p++;
418 }
419
420 targetRow += targetPitch;
421 }
422 }
423 }
424
425
426 void ApplyPerspectiveTransform(Orthanc::ImageAccessor& target,
427 const Orthanc::ImageAccessor& source,
428 const Matrix& a,
429 ImageInterpolation interpolation)
430 {
431 if (source.GetFormat() != target.GetFormat())
432 {
433 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
434 }
435
436 if (a.size1() != 3 ||
437 a.size2() != 3)
438 {
439 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageSize);
440 }
441
442 if (interpolation != ImageInterpolation_Nearest &&
443 interpolation != ImageInterpolation_Bilinear)
444 {
445 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
446 }
447
448 // Check whether we are dealing with an affine transform
449 if (LinearAlgebra::IsCloseToZero(a(2, 0)) &&
450 LinearAlgebra::IsCloseToZero(a(2, 1)))
451 {
452 double w = a(2, 2);
453 if (LinearAlgebra::IsCloseToZero(w))
454 {
455 LOG(ERROR) << "Singular perspective matrix";
456 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
457 }
458 else
459 {
460 ApplyAffineTransform(target, source,
461 a(0, 0) / w, a(0, 1) / w, a(0, 2) / w,
462 a(1, 0) / w, a(1, 1) / w, a(1, 2) / w,
463 interpolation);
464 return;
465 }
466 }
467
468 if (target.GetFormat() == Orthanc::PixelFormat_RGB24)
469 {
470 Orthanc::ImageProcessing::Set(target, 0, 0, 0, 255);
471 }
472 else
473 {
474 Orthanc::ImageProcessing::Set(target, 0);
475 }
476
477 Matrix inva;
478 if (!LinearAlgebra::InvertMatrixUnsafe(inva, a))
479 {
480 return;
481 }
482
483 switch (source.GetFormat())
484 {
485 case Orthanc::PixelFormat_Grayscale8:
486 switch (interpolation)
487 {
488 case ImageInterpolation_Nearest:
489 ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale8,
490 ImageInterpolation_Nearest>(target, source, a, inva);
491 break;
492
493 case ImageInterpolation_Bilinear:
494 ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale8,
495 ImageInterpolation_Bilinear>(target, source, a, inva);
496 break;
497
498 default:
499 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
500 }
501 break;
502
503 case Orthanc::PixelFormat_Grayscale16:
504 switch (interpolation)
505 {
506 case ImageInterpolation_Nearest:
507 ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale16,
508 ImageInterpolation_Nearest>(target, source, a, inva);
509 break;
510
511 case ImageInterpolation_Bilinear:
512 ApplyPerspectiveInternal<Orthanc::PixelFormat_Grayscale16,
513 ImageInterpolation_Bilinear>(target, source, a, inva);
514 break;
515
516 default:
517 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
518 }
519 break;
520
521 case Orthanc::PixelFormat_SignedGrayscale16:
522 switch (interpolation)
523 {
524 case ImageInterpolation_Nearest:
525 ApplyPerspectiveInternal<Orthanc::PixelFormat_SignedGrayscale16,
526 ImageInterpolation_Nearest>(target, source, a, inva);
527 break;
528
529 case ImageInterpolation_Bilinear:
530 ApplyPerspectiveInternal<Orthanc::PixelFormat_SignedGrayscale16,
531 ImageInterpolation_Bilinear>(target, source, a, inva);
532 break;
533
534 default:
535 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
536 }
537 break;
538
539 case Orthanc::PixelFormat_RGB24:
540 switch (interpolation)
541 {
542 case ImageInterpolation_Nearest:
543 ApplyPerspectiveInternal<Orthanc::PixelFormat_RGB24,
544 ImageInterpolation_Nearest>(target, source, a, inva);
545 break;
546
547 default:
548 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
549 }
550 break;
551
552 default:
553 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
554 }
555 }
556 }