Mercurial > hg > orthanc-stone

comparison OrthancStone/Sources/Toolbox/ImageGeometry.cpp @ 1512:244ad1e4e76a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
reorganization of folders
author Sebastien Jodogne <s.jodogne@gmail.com>
date Tue, 07 Jul 2020 16:21:02 +0200
parents Framework/Toolbox/ImageGeometry.cpp@30deba7bc8e2
children 4fb8fdf03314
comparison
equal deleted inserted replaced
1511:9dfeee74c1e6 1512:244ad1e4e76a
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-2020 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 <Images/ImageProcessing.h>
28 #include <Logging.h>
29 #include <OrthancException.h>
30
31
32 namespace OrthancStone
33 {
34 static void AddTransformedPoint(Extent2D& extent,
35 const Matrix& a,
36 double x,
37 double y)
38 {
39 assert(a.size1() == 3 &&
40 a.size2() == 3);
41
42 Vector p = LinearAlgebra::Product(a, LinearAlgebra::CreateVector(x, y, 1));
43
44 if (!LinearAlgebra::IsCloseToZero(p[2]))
45 {
46 extent.AddPoint(p[0] / p[2], p[1] / p[2]);
47 }
48 }
49
50
51 bool GetProjectiveTransformExtent(unsigned int& x1,
52 unsigned int& y1,
53 unsigned int& x2,
54 unsigned int& y2,
55 const Matrix& a,
56 unsigned int sourceWidth,
57 unsigned int sourceHeight,
58 unsigned int targetWidth,
59 unsigned int targetHeight)
60 {
61 if (targetWidth == 0 ||
62 targetHeight == 0)
63 {
64 return false;
65 }
66
67 Extent2D extent;
68 AddTransformedPoint(extent, a, 0, 0);
69 AddTransformedPoint(extent, a, sourceWidth, 0);
70 AddTransformedPoint(extent, a, 0, sourceHeight);
71 AddTransformedPoint(extent, a, sourceWidth, sourceHeight);
72
73 if (extent.IsEmpty())
74 {
75 return false;
76 }
77 else
78 {
79 int tmp = static_cast<int>(std::floor(extent.GetX1()));
80 if (tmp < 0)
81 {
82 x1 = 0;
83 }
84 else
85 {
86 x1 = static_cast<unsigned int>(tmp);
87 }
88
89 tmp = static_cast<int>(std::floor(extent.GetY1()));
90 if (tmp < 0)
91 {
92 y1 = 0;
93 }
94 else
95 {
96 y1 = static_cast<unsigned int>(tmp);
97 }
98
99 tmp = static_cast<int>(std::ceil(extent.GetX2()));
100 if (tmp < 0)
101 {
102 return false;
103 }
104 else if (static_cast<unsigned int>(tmp) >= targetWidth)
105 {
106 x2 = targetWidth - 1;
107 }
108 else
109 {
110 x2 = static_cast<unsigned int>(tmp);
111 }
112
113 tmp = static_cast<int>(std::ceil(extent.GetY2()));
114 if (tmp < 0)
115 {
116 return false;
117 }
118 else if (static_cast<unsigned int>(tmp) >= targetHeight)
119 {
120 y2 = targetHeight - 1;
121 }
122 else
123 {
124 y2 = static_cast<unsigned int>(tmp);
125 }
126
127 return (x1 <= x2 &&
128 y1 <= y2);
129 }
130 }
131
132
133 template <typename Reader,
134 bool HasOffsetX,
135 bool HasOffsetY>
136 static void ApplyAffineTransformToRow(typename Reader::PixelType* p,
137 Reader& reader,
138 unsigned int x1,
139 unsigned int x2,
140 float positionX,
141 float positionY,
142 float offsetX,
143 float offsetY)
144 {
145 typename Reader::PixelType value;
146
147 for (unsigned int x = x1; x <= x2; x++, p++)
148 {
149 if (reader.GetValue(value, positionX, positionY))
150 {
151 *p = value;
152 }
153
154 if (HasOffsetX)
155 {
156 positionX += offsetX;
157 }
158
159 if (HasOffsetY)
160 {
161 positionY += offsetY;
162 }
163 }
164 }
165
166
167 template <Orthanc::PixelFormat Format,
168 ImageInterpolation Interpolation>
169 static void ApplyAffineInternal(Orthanc::ImageAccessor& target,
170 const Orthanc::ImageAccessor& source,
171 const Matrix& a,
172 bool clear)
173 {
174 assert(target.GetFormat() == Format &&
175 source.GetFormat() == Format);
176
177 typedef SubpixelReader<Format, Interpolation> Reader;
178 typedef typename Reader::PixelType PixelType;
179
180 if (clear)
181 {
182 if (Format == Orthanc::PixelFormat_RGB24)
183 {
184 Orthanc::ImageProcessing::Set(target, 0, 0, 0, 255);
185 }
186 else
187 {
188 Orthanc::ImageProcessing::Set(target, 0);
189 }
190 }
191
192 Matrix inva;
193 if (!LinearAlgebra::InvertMatrixUnsafe(inva, a))
194 {
195 // Singular matrix
196 return;
197 }
198
199 Reader reader(source);
200
201 unsigned int x1, y1, x2, y2;
202
203 if (GetProjectiveTransformExtent(x1, y1, x2, y2, a,
204 source.GetWidth(), source.GetHeight(),
205 target.GetWidth(), target.GetHeight()))
206 {
207 const size_t targetPitch = target.GetPitch();
208 uint8_t *targetRow = reinterpret_cast<uint8_t*>
209 (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 bool clear)
263 {
264 if (source.GetFormat() != target.GetFormat())
265 {
266 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
267 }
268
269 if (interpolation != ImageInterpolation_Nearest &&
270 interpolation != ImageInterpolation_Bilinear)
271 {
272 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
273 }
274
275 Matrix a;
276 a.resize(3, 3);
277 a(0, 0) = a11;
278 a(0, 1) = a12;
279 a(0, 2) = b1;
280 a(1, 0) = a21;
281 a(1, 1) = a22;
282 a(1, 2) = b2;
283 a(2, 0) = 0;
284 a(2, 1) = 0;
285 a(2, 2) = 1;
286
287 switch (source.GetFormat())
288 {
289 case Orthanc::PixelFormat_Grayscale8:
290 switch (interpolation)
291 {
292 case ImageInterpolation_Nearest:
293 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale8,
294 ImageInterpolation_Nearest>(target, source, a, clear);
295 break;
296
297 case ImageInterpolation_Bilinear:
298 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale8,
299 ImageInterpolation_Bilinear>(target, source, a, clear);
300 break;
301
302 default:
303 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
304 }
305 break;
306
307 case Orthanc::PixelFormat_Grayscale16:
308 switch (interpolation)
309 {
310 case ImageInterpolation_Nearest:
311 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale16,
312 ImageInterpolation_Nearest>(target, source, a, clear);
313 break;
314
315 case ImageInterpolation_Bilinear:
316 ApplyAffineInternal<Orthanc::PixelFormat_Grayscale16,
317 ImageInterpolation_Bilinear>(target, source, a, clear);
318 break;
319
320 default:
321 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
322 }
323 break;
324
325 case Orthanc::PixelFormat_SignedGrayscale16:
326 switch (interpolation)
327 {
328 case ImageInterpolation_Nearest:
329 ApplyAffineInternal<Orthanc::PixelFormat_SignedGrayscale16,
330 ImageInterpolation_Nearest>(target, source, a, clear);
331 break;
332
333 case ImageInterpolation_Bilinear:
334 ApplyAffineInternal<Orthanc::PixelFormat_SignedGrayscale16,
335 ImageInterpolation_Bilinear>(target, source, a, clear);
336 break;
337
338 default:
339 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
340 }
341 break;
342
343 case Orthanc::PixelFormat_Float32:
344 switch (interpolation)
345 {
346 case ImageInterpolation_Nearest:
347 ApplyAffineInternal<Orthanc::PixelFormat_Float32,
348 ImageInterpolation_Nearest>(target, source, a, clear);
349 break;
350
351 case ImageInterpolation_Bilinear:
352 ApplyAffineInternal<Orthanc::PixelFormat_Float32,
353 ImageInterpolation_Bilinear>(target, source, a, clear);
354 break;
355
356 default:
357 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
358 }
359 break;
360
361 case Orthanc::PixelFormat_RGB24:
362 switch (interpolation)
363 {
364 case ImageInterpolation_Nearest:
365 ApplyAffineInternal<Orthanc::PixelFormat_RGB24,
366 ImageInterpolation_Nearest>(target, source, a, clear);
367 break;
368
369 default:
370 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
371 }
372 break;
373
374 default:
375 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
376 }
377 }
378
379
380 template <Orthanc::PixelFormat Format,
381 ImageInterpolation Interpolation>
382 static void ApplyProjectiveInternal(Orthanc::ImageAccessor& target,
383 const Orthanc::ImageAccessor& source,
384 const Matrix& a,
385 const Matrix& inva)
386 {
387 assert(target.GetFormat() == Format &&
388 source.GetFormat() == Format);
389
390 typedef SubpixelReader<Format, Interpolation> Reader;
391 typedef typename Reader::PixelType PixelType;
392
393 Reader reader(source);
394 unsigned int x1, y1, x2, y2;
395
396 const float floatWidth = static_cast<float>(source.GetWidth());
397 const float floatHeight = static_cast<float>(source.GetHeight());
398
399 if (GetProjectiveTransformExtent(x1, y1, x2, y2, a,
400 source.GetWidth(), source.GetHeight(),
401 target.GetWidth(), target.GetHeight()))
402 {
403 const size_t targetPitch = target.GetPitch();
404 uint8_t *targetRow = reinterpret_cast<uint8_t*>
405 (reinterpret_cast<PixelType*>(target.GetRow(y1)) + x1);
406
407 for (unsigned int y = y1; y <= y2; y++)
408 {
409 PixelType *p = reinterpret_cast<PixelType*>(targetRow);
410
411 for (unsigned int x = x1; x <= x2; x++)
412 {
413 Vector v;
414 LinearAlgebra::AssignVector(v, static_cast<double>(x) + 0.5,
415 static_cast<double>(y) + 0.5, 1);
416
417 Vector vv = LinearAlgebra::Product(inva, v);
418
419 assert(!LinearAlgebra::IsCloseToZero(vv[2]));
420 const double w = 1.0 / vv[2];
421 const float sourceX = static_cast<float>(vv[0] * w);
422 const float sourceY = static_cast<float>(vv[1] * w);
423
424 // Make sure no integer overflow will occur after truncation
425 // (the static_cast<unsigned int> could otherwise throw an
426 // exception in WebAssembly if strong projective effects)
427 if (sourceX < floatWidth &&
428 sourceY < floatHeight)
429 {
430 reader.GetValue(*p, sourceX, sourceY);
431 }
432
433 p++;
434 }
435
436 targetRow += targetPitch;
437 }
438 }
439 }
440
441
442 void ApplyProjectiveTransform(Orthanc::ImageAccessor& target,
443 const Orthanc::ImageAccessor& source,
444 const Matrix& a,
445 ImageInterpolation interpolation,
446 bool clear)
447 {
448 if (source.GetFormat() != target.GetFormat())
449 {
450 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat);
451 }
452
453 if (a.size1() != 3 ||
454 a.size2() != 3)
455 {
456 throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageSize);
457 }
458
459 if (interpolation != ImageInterpolation_Nearest &&
460 interpolation != ImageInterpolation_Bilinear)
461 {
462 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
463 }
464
465 // Check whether we are dealing with an affine transform
466 if (LinearAlgebra::IsCloseToZero(a(2, 0)) &&
467 LinearAlgebra::IsCloseToZero(a(2, 1)))
468 {
469 double w = a(2, 2);
470 if (LinearAlgebra::IsCloseToZero(w))
471 {
472 LOG(ERROR) << "Singular projective matrix";
473 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
474 }
475 else
476 {
477 ApplyAffineTransform(target, source,
478 a(0, 0) / w, a(0, 1) / w, a(0, 2) / w,
479 a(1, 0) / w, a(1, 1) / w, a(1, 2) / w,
480 interpolation, clear);
481 return;
482 }
483 }
484
485 if (clear)
486 {
487 if (target.GetFormat() == Orthanc::PixelFormat_RGB24)
488 {
489 Orthanc::ImageProcessing::Set(target, 0, 0, 0, 255);
490 }
491 else
492 {
493 Orthanc::ImageProcessing::Set(target, 0);
494 }
495 }
496
497 Matrix inva;
498 if (!LinearAlgebra::InvertMatrixUnsafe(inva, a))
499 {
500 return;
501 }
502
503 switch (source.GetFormat())
504 {
505 case Orthanc::PixelFormat_Grayscale8:
506 switch (interpolation)
507 {
508 case ImageInterpolation_Nearest:
509 ApplyProjectiveInternal<Orthanc::PixelFormat_Grayscale8,
510 ImageInterpolation_Nearest>(target, source, a, inva);
511 break;
512
513 case ImageInterpolation_Bilinear:
514 ApplyProjectiveInternal<Orthanc::PixelFormat_Grayscale8,
515 ImageInterpolation_Bilinear>(target, source, a, inva);
516 break;
517
518 default:
519 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
520 }
521 break;
522
523 case Orthanc::PixelFormat_Grayscale16:
524 switch (interpolation)
525 {
526 case ImageInterpolation_Nearest:
527 ApplyProjectiveInternal<Orthanc::PixelFormat_Grayscale16,
528 ImageInterpolation_Nearest>(target, source, a, inva);
529 break;
530
531 case ImageInterpolation_Bilinear:
532 ApplyProjectiveInternal<Orthanc::PixelFormat_Grayscale16,
533 ImageInterpolation_Bilinear>(target, source, a, inva);
534 break;
535
536 default:
537 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
538 }
539 break;
540
541 case Orthanc::PixelFormat_SignedGrayscale16:
542 switch (interpolation)
543 {
544 case ImageInterpolation_Nearest:
545 ApplyProjectiveInternal<Orthanc::PixelFormat_SignedGrayscale16,
546 ImageInterpolation_Nearest>(target, source, a, inva);
547 break;
548
549 case ImageInterpolation_Bilinear:
550 ApplyProjectiveInternal<Orthanc::PixelFormat_SignedGrayscale16,
551 ImageInterpolation_Bilinear>(target, source, a, inva);
552 break;
553
554 default:
555 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
556 }
557 break;
558
559 case Orthanc::PixelFormat_Float32:
560 switch (interpolation)
561 {
562 case ImageInterpolation_Nearest:
563 ApplyProjectiveInternal<Orthanc::PixelFormat_Float32,
564 ImageInterpolation_Nearest>(target, source, a, inva);
565 break;
566
567 case ImageInterpolation_Bilinear:
568 ApplyProjectiveInternal<Orthanc::PixelFormat_Float32,
569 ImageInterpolation_Bilinear>(target, source, a, inva);
570 break;
571
572 default:
573 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
574 }
575 break;
576
577 case Orthanc::PixelFormat_RGB24:
578 switch (interpolation)
579 {
580 case ImageInterpolation_Nearest:
581 ApplyProjectiveInternal<Orthanc::PixelFormat_RGB24,
582 ImageInterpolation_Nearest>(target, source, a, inva);
583 break;
584
585 default:
586 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
587 }
588 break;
589
590 default:
591 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
592 }
593 }
594 }