1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2008-2015 Imperial College London
5 * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel
6 * Copyright 2013-2015 Andreas Schuh
7 *
8 * Licensed under the Apache License, Version 2.0 (the "License");
9 * you may not use this file except in compliance with the License.
10 * You may obtain a copy of the License at
11 *
12 * http://www.apache.org/licenses/LICENSE-2.0
13 *
14 * Unless required by applicable law or agreed to in writing, software
15 * distributed under the License is distributed on an "AS IS" BASIS,
16 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 * See the License for the specific language governing permissions and
18 * limitations under the License.
19 */
20
21 #include "mirtk/ImageAttributes.h"
22
23 #include "mirtk/Math.h"
24 #include "mirtk/Array.h"
25 #include "mirtk/Memory.h"
26 #include "mirtk/Stream.h"
27 #include "mirtk/Indent.h"
28 #include "mirtk/Vector.h"
29 #include "mirtk/Vector3.h"
30 #include "mirtk/Matrix.h"
31 #include "mirtk/Matrix3x3.h"
32
33
34 namespace mirtk {
35
36
37 // -----------------------------------------------------------------------------
ImageAttributes()38 ImageAttributes::ImageAttributes()
39 {
40 // Default image size
41 _x = 0;
42 _y = 0;
43 _z = 1;
44 _t = 1;
45
46 // Default voxel size
47 _dx = 1;
48 _dy = 1;
49 _dz = 1;
50 _dt = 1;
51
52 // Default origin
53 _xorigin = 0;
54 _yorigin = 0;
55 _zorigin = 0;
56 _torigin = 0;
57
58 // Default x-axis
59 _xaxis[0] = 1;
60 _xaxis[1] = 0;
61 _xaxis[2] = 0;
62
63 // Default y-axis
64 _yaxis[0] = 0;
65 _yaxis[1] = 1;
66 _yaxis[2] = 0;
67
68 // Default z-axis
69 _zaxis[0] = 0;
70 _zaxis[1] = 0;
71 _zaxis[2] = 1;
72
73 // Default affine transformation
74 _smat.Initialize(4, 4);
75 _smat.Ident();
76
77 _i2w = nullptr;
78 _w2i = nullptr;
79 }
80
81 // -----------------------------------------------------------------------------
ImageAttributes(int x,int y,double dx,double dy)82 ImageAttributes::ImageAttributes(int x, int y, double dx, double dy)
83 {
84 // Default image size
85 _x = x;
86 _y = y;
87 _z = 1;
88 _t = 1;
89
90 // Default voxel size
91 _dx = dx;
92 _dy = dy;
93 _dz = 1;
94 _dt = 1;
95
96 // Default origin
97 _xorigin = 0;
98 _yorigin = 0;
99 _zorigin = 0;
100 _torigin = 0;
101
102 // Default x-axis
103 _xaxis[0] = 1;
104 _xaxis[1] = 0;
105 _xaxis[2] = 0;
106
107 // Default y-axis
108 _yaxis[0] = 0;
109 _yaxis[1] = 1;
110 _yaxis[2] = 0;
111
112 // Default z-axis
113 _zaxis[0] = 0;
114 _zaxis[1] = 0;
115 _zaxis[2] = 1;
116
117 // Default affine transformation
118 _smat.Initialize(4, 4);
119 _smat.Ident();
120
121 _i2w = nullptr;
122 _w2i = nullptr;
123 }
124
125 // -----------------------------------------------------------------------------
ImageAttributes(int x,int y,int z,double dx,double dy,double dz)126 ImageAttributes::ImageAttributes(int x, int y, int z, double dx, double dy, double dz)
127 {
128 _x = x;
129 _y = y;
130 _z = z;
131 _t = 1;
132
133 _dx = dx;
134 _dy = dy;
135 _dz = dz;
136 _dt = 1;
137
138 _xorigin = 0;
139 _yorigin = 0;
140 _zorigin = 0;
141 _torigin = 0;
142
143 _xaxis[0] = 1;
144 _xaxis[1] = 0;
145 _xaxis[2] = 0;
146
147 _yaxis[0] = 0;
148 _yaxis[1] = 1;
149 _yaxis[2] = 0;
150
151 _zaxis[0] = 0;
152 _zaxis[1] = 0;
153 _zaxis[2] = 1;
154
155 _smat.Initialize(4, 4);
156 _smat.Ident();
157
158 _i2w = nullptr;
159 _w2i = nullptr;
160 }
161
162 // -----------------------------------------------------------------------------
ImageAttributes(int x,int y,int z,int t,double dx,double dy,double dz,double dt)163 ImageAttributes::ImageAttributes(int x, int y, int z, int t, double dx, double dy, double dz, double dt)
164 {
165 _x = x;
166 _y = y;
167 _z = z;
168 _t = t;
169
170 _dx = dx;
171 _dy = dy;
172 _dz = dz;
173 _dt = dt;
174
175 _xorigin = 0;
176 _yorigin = 0;
177 _zorigin = 0;
178 _torigin = 0;
179
180 _xaxis[0] = 1;
181 _xaxis[1] = 0;
182 _xaxis[2] = 0;
183
184 _yaxis[0] = 0;
185 _yaxis[1] = 1;
186 _yaxis[2] = 0;
187
188 _zaxis[0] = 0;
189 _zaxis[1] = 0;
190 _zaxis[2] = 1;
191
192 _smat.Initialize(4, 4);
193 _smat.Ident();
194
195 _i2w = nullptr;
196 _w2i = nullptr;
197 }
198
199 // -----------------------------------------------------------------------------
ImageAttributes(const ImageAttributes & attr)200 ImageAttributes::ImageAttributes(const ImageAttributes &attr)
201 {
202 _x = attr._x;
203 _y = attr._y;
204 _z = attr._z;
205 _t = attr._t;
206
207 _dx = attr._dx;
208 _dy = attr._dy;
209 _dz = attr._dz;
210 _dt = attr._dt;
211
212 _xorigin = attr._xorigin;
213 _yorigin = attr._yorigin;
214 _zorigin = attr._zorigin;
215 _torigin = attr._torigin;
216
217 _xaxis[0] = attr._xaxis[0];
218 _xaxis[1] = attr._xaxis[1];
219 _xaxis[2] = attr._xaxis[2];
220
221 _yaxis[0] = attr._yaxis[0];
222 _yaxis[1] = attr._yaxis[1];
223 _yaxis[2] = attr._yaxis[2];
224
225 _zaxis[0] = attr._zaxis[0];
226 _zaxis[1] = attr._zaxis[1];
227 _zaxis[2] = attr._zaxis[2];
228
229 _smat = attr._smat;
230
231 // Do not copy pointers to pre-computed matrices as we don't know how
232 // long the objects pointed to by the other attributes will be valid
233 _i2w = nullptr;
234 _w2i = nullptr;
235 }
236
237 // -----------------------------------------------------------------------------
operator =(const ImageAttributes & attr)238 ImageAttributes &ImageAttributes::operator =(const ImageAttributes &attr)
239 {
240 if (this != &attr) {
241 _x = attr._x;
242 _y = attr._y;
243 _z = attr._z;
244 _t = attr._t;
245
246 _dx = attr._dx;
247 _dy = attr._dy;
248 _dz = attr._dz;
249 _dt = attr._dt;
250
251 _xorigin = attr._xorigin;
252 _yorigin = attr._yorigin;
253 _zorigin = attr._zorigin;
254 _torigin = attr._torigin;
255
256 _xaxis[0] = attr._xaxis[0];
257 _xaxis[1] = attr._xaxis[1];
258 _xaxis[2] = attr._xaxis[2];
259
260 _yaxis[0] = attr._yaxis[0];
261 _yaxis[1] = attr._yaxis[1];
262 _yaxis[2] = attr._yaxis[2];
263
264 _zaxis[0] = attr._zaxis[0];
265 _zaxis[1] = attr._zaxis[1];
266 _zaxis[2] = attr._zaxis[2];
267
268 _smat = attr._smat;
269 }
270
271 // Do not copy pointers to pre-computed matrices as we don't know how
272 // long the objects pointed to by the other attributes will be valid
273 //
274 // Reset pointers even when assigning to itself such that behavior is
275 // the same as when assigning from another instance
276 _i2w = nullptr;
277 _w2i = nullptr;
278
279 return *this;
280 }
281
282 // -----------------------------------------------------------------------------
LatticeToWorld(double * x,double * y,double * z) const283 void ImageAttributes::LatticeToWorld(double *x, double *y, double *z) const
284 {
285 int idx = 0;
286 for (int k = 0; k < _z; ++k)
287 for (int j = 0; j < _y; ++j)
288 for (int i = 0; i < _x; ++i, ++idx) {
289 x[idx] = i, y[idx] = j, z[idx] = k;
290 LatticeToWorld(x[idx], y[idx], z[idx]);
291 }
292 }
293
294 // -----------------------------------------------------------------------------
LatticeToWorld(double * x,double * y,double * z,double * t) const295 void ImageAttributes::LatticeToWorld(double *x, double *y, double *z, double *t) const
296 {
297 const int xyz = _x * _y * _z;
298
299 LatticeToWorld(x, y, z);
300
301 if (_dt == .0) {
302 double tl = LatticeToTime(0);
303 for (int idx = 0; idx < xyz; ++idx) {
304 (*t++) = tl;
305 }
306 } else {
307 int idx = xyz;
308 for (int l = 1; l < _t; ++l) {
309 memcpy(x + idx, x, xyz * sizeof(double));
310 memcpy(y + idx, y, xyz * sizeof(double));
311 memcpy(z + idx, z, xyz * sizeof(double));
312 idx += xyz;
313 }
314 for (int l = 0; l < _t; ++l) {
315 double tl = LatticeToTime(l);
316 for (int idx = 0; idx < xyz; ++idx) {
317 (*t++) = tl;
318 }
319 }
320 }
321 }
322
323 // -----------------------------------------------------------------------------
ContainsInSpace(const ImageAttributes & attr) const324 bool ImageAttributes::ContainsInSpace(const ImageAttributes &attr) const
325 {
326 double x = 0, y = 0, z = 0;
327 attr.LatticeToWorld(x, y, z);
328 this->WorldToLattice(x, y, z);
329 if (x <= -.5 || x >= _x - .5 || y <= -.5 || y >= _y - .5 || z <= -.5 || z >= _z - .5) return false;
330 x = attr._x - 1, y = attr._y - 1, z = attr._z - 1;
331 attr.LatticeToWorld(x, y, z);
332 this->WorldToLattice(x, y, z);
333 if (x <= -.5 || x >= _x - .5 || y <= -.5 || y >= _y - .5 || z <= -.5 || z >= _z - .5) return false;
334 return true;
335 }
336
337 // -----------------------------------------------------------------------------
EqualInSpace(const ImageAttributes & attr) const338 bool ImageAttributes::EqualInSpace(const ImageAttributes &attr) const
339 {
340 // Spatial dimensions
341 if (_x != attr._x || _y != attr._y || _z != attr._z) return false;
342 // Spatial resolution
343 if (!fequal(_dx, attr._dx) || !fequal(_dy, attr._dy) || !fequal(_dz, attr._dz)) return false;
344 // Spatial orientation
345 for (int i = 0; i < 3; ++i) if (!fequal(_xaxis[i], attr._xaxis[i])) return false;
346 for (int i = 0; i < 3; ++i) if (!fequal(_yaxis[i], attr._yaxis[i])) return false;
347 for (int i = 0; i < 3; ++i) if (!fequal(_zaxis[i], attr._zaxis[i])) return false;
348 // Spatial origin
349 if (!fequal(_xorigin, attr._xorigin) || !fequal(_yorigin, attr._yorigin) || !fequal(_zorigin, attr._zorigin)) return false;
350 // Spatial transformation
351 for (int r = 0; r < 4; ++r) {
352 for (int c = 0; c < 4; ++c) {
353 if (!fequal(_smat(r, c), attr._smat(r, c))) return false;
354 }
355 }
356 return true;
357 }
358
359 // -----------------------------------------------------------------------------
EqualInTime(const ImageAttributes & attr) const360 bool ImageAttributes::EqualInTime(const ImageAttributes &attr) const
361 {
362 return (_t == attr._t) && fequal(_dt, attr._dt) && fequal(_torigin, attr._torigin);
363 }
364
365 // -----------------------------------------------------------------------------
Area() const366 double ImageAttributes::Area() const
367 {
368 return _x * _dx * _y * _dy;
369 }
370
371 // -----------------------------------------------------------------------------
Volume() const372 double ImageAttributes::Volume() const
373 {
374 return _x * _dx * _y * _dy * _z * _dz;
375 }
376
377 // -----------------------------------------------------------------------------
Space() const378 double ImageAttributes::Space() const
379 {
380 if (static_cast<bool>(*this)) {
381 if (_dz > 0.) return Volume();
382 else return Area();
383 } else {
384 return 0.;
385 }
386 }
387
388 // -----------------------------------------------------------------------------
Print(ostream & os,Indent indent) const389 void ImageAttributes::Print(ostream &os, Indent indent) const
390 {
391 // Change output stream settings
392 const streamsize w = os.width (0);
393 const streamsize p = os.precision(5);
394 const ios::fmtflags f = os.flags ();
395 // Print attributes of lattice
396 bool sz = (_z > 1 && _dz);
397 bool st = (_t > 1 && _dt);
398 bool dz = (_dz && (_dz != 1.0 || _z > 1));
399 bool dt = (_dt && (_dt != 1.0 || _t > 1));
400 if (_t > 1 && !_dt) os << indent << "Channels: " << setw(10) << _t << "\n";
401 os << indent << "Size: " << setw(10) << _x
402 << " x " << setw(10) << _y;
403 if (sz || st) os << " x " << setw(10) << _z;
404 if ( st) os << " x " << setw(10) << _t;
405 os << "\n";
406 os << indent << "Spacing: " << fixed << setw(10) << _dx
407 << " x " << fixed << setw(10) << _dy;
408 if (dz || dt) os << " x " << fixed << setw(10) << _dz;
409 if ( dt) os << " x " << fixed << setw(10) << _dt;
410 os << "\n";
411 os << indent << "Origin: [" << fixed << setw(10) << _xorigin
412 << " , " << fixed << setw(10) << _yorigin
413 << " , " << fixed << setw(10) << _zorigin
414 << " , " << fixed << setw(10) << _torigin
415 << "]\n";
416 os << indent << "X-axis: [" << fixed << setw(10) << _xaxis[0]
417 << " , " << fixed << setw(10) << _xaxis[1]
418 << " , " << fixed << setw(10) << _xaxis[2]
419 << "]\n";
420 os << indent << "Y-axis: [" << fixed << setw(10) << _yaxis[0]
421 << " , " << fixed << setw(10) << _yaxis[1]
422 << " , " << fixed << setw(10) << _yaxis[2]
423 << "]\n";
424 os << indent << "Z-axis: [" << fixed << setw(10) << _zaxis[0]
425 << " , " << fixed << setw(10) << _zaxis[1]
426 << " , " << fixed << setw(10) << _zaxis[2]
427 << "]\n";
428 if (!_smat.IsIdentity()) {
429 os << indent << "Affine";
430 _smat.Print(os, indent);
431 }
432 // Restore output stream settings
433 os.width (w);
434 os.precision(p);
435 os.flags (f);
436 }
437
438 // -----------------------------------------------------------------------------
Print(Indent indent) const439 void ImageAttributes::Print(Indent indent) const
440 {
441 Print(cout, indent);
442 }
443
444 // -----------------------------------------------------------------------------
GetLatticeToWorldMatrix() const445 Matrix ImageAttributes::GetLatticeToWorldMatrix() const
446 {
447 // Final mat = A * T * R * S * T0
448 Matrix mat(4, 4), tmp(4, 4);
449
450 // T0: Translate image origin
451 mat.Ident();
452 mat(0, 3) = - (_x - 1) / 2.0;
453 mat(1, 3) = - (_y - 1) / 2.0;
454 mat(2, 3) = - (_z - 1) / 2.0;
455
456 // S: Convert to world units
457 tmp.Ident();
458 tmp(0, 0) = _dx;
459 tmp(1, 1) = _dy;
460 tmp(2, 2) = _dz;
461 mat = tmp * mat;
462
463 // R: Orientation
464 tmp(0, 0) = _xaxis[0];
465 tmp(1, 0) = _xaxis[1];
466 tmp(2, 0) = _xaxis[2];
467 tmp(0, 1) = _yaxis[0];
468 tmp(1, 1) = _yaxis[1];
469 tmp(2, 1) = _yaxis[2];
470 tmp(0, 2) = _zaxis[0];
471 tmp(1, 2) = _zaxis[1];
472 tmp(2, 2) = _zaxis[2];
473 mat = tmp * mat;
474
475 // T: Translate world origin
476 tmp.Ident();
477 tmp(0, 3) = _xorigin;
478 tmp(1, 3) = _yorigin;
479 tmp(2, 3) = _zorigin;
480 mat = tmp * mat;
481
482 // A: Transform world coordinate by optional affine transformation
483 // to align (scanner) coordinates with another anatomy
484 mat = _smat * mat;
485
486 return mat;
487 }
488
489 // -----------------------------------------------------------------------------
GetWorldToLatticeMatrix() const490 Matrix ImageAttributes::GetWorldToLatticeMatrix() const
491 {
492 // Final mat = inv(A * T * R * S * T0)
493 // = inv(T0) * inv(S) * inv(R) * inv(T) * inv(A),
494 Matrix mat(4, 4), tmp(4, 4);
495
496 // inv(A): Transform world coordinate by optional inverse of the affine
497 // transformation which aligns (scanner) coordinates with another anatomy
498 if (_smat.IsIdentity()) {
499 mat.Ident(); // avoid numerical imprecisions caused by Matrix::Invert
500 } else {
501 mat = _smat.Inverse();
502 }
503
504 // inv(T): Translate world origin
505 tmp.Ident();
506 tmp(0, 3) = - _xorigin;
507 tmp(1, 3) = - _yorigin;
508 tmp(2, 3) = - _zorigin;
509 mat = tmp * mat;
510
511 // inv(R): Orientation
512 tmp(0, 0) = _xaxis[0];
513 tmp(0, 1) = _xaxis[1];
514 tmp(0, 2) = _xaxis[2];
515 tmp(0, 3) = .0;
516 tmp(1, 0) = _yaxis[0];
517 tmp(1, 1) = _yaxis[1];
518 tmp(1, 2) = _yaxis[2];
519 tmp(1, 3) = .0;
520 tmp(2, 0) = _zaxis[0];
521 tmp(2, 1) = _zaxis[1];
522 tmp(2, 2) = _zaxis[2];
523 tmp(2, 3) = .0;
524 mat = tmp * mat;
525
526 // inv(S): Convert to voxel units
527 tmp.Ident();
528 if (_dx) tmp(0, 0) = 1.0 / _dx;
529 if (_dy) tmp(1, 1) = 1.0 / _dy;
530 if (_dz) tmp(2, 2) = 1.0 / _dz;
531 mat = tmp * mat;
532
533 // inv(T0): Translate image origin
534 tmp(0, 0) = tmp(1, 1) = tmp(2, 2) = 1.0;
535 tmp(0, 3) = (_x - 1) / 2.0;
536 tmp(1, 3) = (_y - 1) / 2.0;
537 tmp(2, 3) = (_z - 1) / 2.0;
538 mat = tmp * mat;
539
540 return mat;
541 }
542
543 // -----------------------------------------------------------------------------
PutAffineMatrix(const Matrix & m,bool apply)544 void ImageAttributes::PutAffineMatrix(const Matrix &m, bool apply)
545 {
546 if (m.IsIdentity()) {
547 _smat.Ident();
548 } else if (apply) {
549
550 // Lattice to world matrix is: A * T * R * S * T0
551 const Matrix B = GetLatticeToWorldMatrix();
552
553 // T0: Translate image origin
554 Matrix T0(4, 4);
555 T0.Ident();
556 T0(0, 3) = - (_x - 1) / 2.;
557 T0(1, 3) = - (_y - 1) / 2.;
558 T0(2, 3) = - (_z - 1) / 2.;
559
560 // Get new lattice to world matrix, excluding T0
561 Matrix C = m * B * T0.Inverse();
562
563 // Decompose C into T * R * S * K components
564 // (K: shear, S: scale, R: rotation, T: translation)
565 double tx, ty, tz, rx, ry, rz, sx, sy, sz, sxy, sxz, syz;
566 MatrixToAffineParameters(C, tx, ty, tz, rx, ry, rz, sx, sy, sz, sxy, sxz, syz);
567
568 // Remove shearing part
569 const Matrix K = AffineParametersToMatrix(0., 0., 0., 0., 0., 0., 1., 1., 1, sxy, sxz, syz);
570 const bool has_shearing = !K.IsIdentity();
571 if (has_shearing) {
572 C = C * K.Inverse();
573 }
574
575 // Decompose again, this time without shearing
576 MatrixToAffineParameters(C, tx, ty, tz, rx, ry, rz, sx, sy, sz, sxy, sxz, syz);
577
578 // Voxel size must be positive
579 Matrix flip(4, 4);
580 flip.Ident();
581 if (sx < 0.) {
582 flip(0, 0) = -1.;
583 sx = -sx;
584 }
585 if (sy < 0.) {
586 flip(1, 1) = -1.;
587 sy = -sy;
588 }
589 if (sz < 0.) {
590 flip(2, 2) = -1.;
591 sz = -sz;
592 }
593
594 // Set new lattice attributes
595 _dx = sx;
596 _dy = sy;
597 _dz = sz;
598
599 // R: Orientation
600 const Matrix R = RigidParametersToMatrix(0., 0., 0., rx, ry, rz) * flip;
601 _xaxis[0] = R(0, 0);
602 _xaxis[1] = R(1, 0);
603 _xaxis[2] = R(2, 0);
604 _yaxis[0] = R(0, 1);
605 _yaxis[1] = R(1, 1);
606 _yaxis[2] = R(2, 1);
607 _zaxis[0] = R(0, 2);
608 _zaxis[1] = R(1, 2);
609 _zaxis[2] = R(2, 2);
610
611 // T: Translate world origin
612 _xorigin = tx;
613 _yorigin = ty;
614 _zorigin = tz;
615
616 // Set remaining difference due to shearing as post-image to world mapping
617 if (has_shearing) {
618 _smat = C * K * T0 * GetWorldToLatticeMatrix();
619 }
620 } else {
621 _smat = m;
622 }
623 }
624
625 // =============================================================================
626 // Image domain helpers
627 // =============================================================================
628
629 // -----------------------------------------------------------------------------
630 /// Adjust size of output domain if necessary such that all corners of the
631 /// input image domain are within the bounds of the output image domain
632 ///
633 /// \param[in] in Image grid that should fit into \p out domain.
634 /// \param[in,out] out Image grid which should fully contain \p in without
635 /// changing its orientation nor voxel size.
AdjustFieldOfView(const ImageAttributes & in,ImageAttributes & out)636 void AdjustFieldOfView(const ImageAttributes &in, ImageAttributes &out)
637 {
638 const Matrix i2w = in .GetImageToWorldMatrix();
639 Matrix w2i = out.GetWorldToImageMatrix();
640 const int di = (in._x > 1 ? in._x - 1 : 1);
641 const int dj = (in._y > 1 ? in._y - 1 : 1);
642 const int dk = (in._z > 1 ? in._z - 1 : 1);
643 for (int k1 = 0; k1 < in._z; k1 += dk)
644 for (int j1 = 0; j1 < in._y; j1 += dj)
645 for (int i1 = 0; i1 < in._x; i1 += di) {
646 // Convert corner voxel indices of input domain to world coordinates
647 double x = i2w(0, 0) * i1 + i2w(0, 1) * j1 + i2w(0, 2) * k1 + i2w(0, 3);
648 double y = i2w(1, 0) * i1 + i2w(1, 1) * j1 + i2w(1, 2) * k1 + i2w(1, 3);
649 double z = i2w(2, 0) * i1 + i2w(2, 1) * j1 + i2w(2, 2) * k1 + i2w(2, 3);
650 // Convert world coordinates to voxel index w.r.t output domain
651 double i2 = w2i(0, 0) * x + w2i(0, 1) * y + w2i(0, 2) * z + w2i(0, 3);
652 double j2 = w2i(1, 0) * x + w2i(1, 1) * y + w2i(1, 2) * z + w2i(1, 3);
653 double k2 = w2i(2, 0) * x + w2i(2, 1) * y + w2i(2, 2) * z + w2i(2, 3);
654 // If point is outside the output domain...
655 if (i2 < 0 || i2 > out._x-1 || j2 < 0 || j2 > out._y-1 || k2 < 0 || k2 > out._z-1) {
656 // Increase output domain by difference along the axes directions
657 if (i2 < 0) {
658 int dv = iceil(abs(i2));
659 out._x += dv;
660 out._xorigin -= out._xaxis[0] * out._dx * dv / 2.0;
661 out._yorigin -= out._xaxis[1] * out._dx * dv / 2.0;
662 out._zorigin -= out._xaxis[2] * out._dx * dv / 2.0;
663 } else if (i2 > out._x - 1) {
664 int dv = iceil(i2 - out._x + 1);
665 out._x += dv;
666 out._xorigin += out._xaxis[0] * out._dx * dv / 2.0;
667 out._yorigin += out._xaxis[1] * out._dx * dv / 2.0;
668 out._zorigin += out._xaxis[2] * out._dx * dv / 2.0;
669 }
670 if (j2 < 0) {
671 int dv = iceil(abs(j2));
672 out._y += dv;
673 out._xorigin -= out._yaxis[0] * out._dy * dv / 2.0;
674 out._yorigin -= out._yaxis[1] * out._dy * dv / 2.0;
675 out._zorigin -= out._yaxis[2] * out._dy * dv / 2.0;
676 } else if (j2 > out._y - 1) {
677 int dv = iceil(j2 - out._y + 1);
678 out._y += dv;
679 out._xorigin += out._yaxis[0] * out._dy * dv / 2.0;
680 out._yorigin += out._yaxis[1] * out._dy * dv / 2.0;
681 out._zorigin += out._yaxis[2] * out._dy * dv / 2.0;
682 }
683 if (k2 < 0) {
684 int dv = iceil(abs(k2));
685 out._z += dv;
686 out._xorigin -= out._zaxis[0] * out._dz * dv / 2.0;
687 out._yorigin -= out._zaxis[1] * out._dz * dv / 2.0;
688 out._zorigin -= out._zaxis[2] * out._dz * dv / 2.0;
689 } else if (k2 > out._z - 1) {
690 int dv = iceil(k2 - out._z + 1);
691 out._z += dv;
692 out._xorigin += out._zaxis[0] * out._dz * dv / 2.0;
693 out._yorigin += out._zaxis[1] * out._dz * dv / 2.0;
694 out._zorigin += out._zaxis[2] * out._dz * dv / 2.0;
695 }
696 // Update transformation matrix
697 w2i = out.GetWorldToImageMatrix();
698 }
699 }
700 }
701
702 // -----------------------------------------------------------------------------
GetCorners(const ImageAttributes & attr,double * corners[3])703 void GetCorners(const ImageAttributes &attr, double *corners[3])
704 {
705 int idx = 0;
706 for (int k = 0; k <= 1; ++k)
707 for (int j = 0; j <= 1; ++j)
708 for (int i = 0; i <= 1; ++i, ++idx) {
709 corners[idx][0] = i * (attr._x - 1);
710 corners[idx][1] = j * (attr._y - 1);
711 corners[idx][2] = k * (attr._z - 1);
712 attr.LatticeToWorld(corners[idx][0], corners[idx][1], corners[idx][2]);
713 }
714 }
715
716 // -----------------------------------------------------------------------------
GetCorners(const ImageAttributes & attr,double corners[][3])717 void GetCorners(const ImageAttributes &attr, double corners[][3])
718 {
719 int idx = 0;
720 for (int k = 0; k <= 1; ++k)
721 for (int j = 0; j <= 1; ++j)
722 for (int i = 0; i <= 1; ++i, ++idx) {
723 corners[idx][0] = i * (attr._x - 1);
724 corners[idx][1] = j * (attr._y - 1);
725 corners[idx][2] = k * (attr._z - 1);
726 attr.LatticeToWorld(corners[idx][0], corners[idx][1], corners[idx][2]);
727 }
728 }
729
730 // -----------------------------------------------------------------------------
OrthogonalFieldOfView(const ImageAttributes & attr)731 ImageAttributes OrthogonalFieldOfView(const ImageAttributes &attr)
732 {
733 ImageAttributes out = attr;
734 // Apply any additional affine transformation to the attributes of the image domain
735 if (!attr._smat.IsIdentity()) {
736 // Reset affine transformation of output attributes
737 out._smat.Ident();
738 // Get corners of transformed image lattice
739 double corners[8][3];
740 GetCorners(attr, corners);
741 // Set output origin
742 out._xorigin = out._yorigin = out._zorigin = .0;
743 for (int i = 0; i < 8; ++i) {
744 out._xorigin += corners[i][0];
745 out._yorigin += corners[i][1];
746 out._zorigin += corners[i][2];
747 }
748 out._xorigin /= 8, out._yorigin /= 8, out._zorigin /= 8;
749 // Set output orientation
750 Matrix3x3 covar(.0);
751 for (int i = 0; i < 8; ++i) {
752 corners[i][0] -= out._xorigin;
753 corners[i][1] -= out._yorigin;
754 corners[i][2] -= out._zorigin;
755 for (int r = 0; r < 3; ++r)
756 for (int c = 0; c < 3; ++c) {
757 covar[r][c] += corners[i][r] * corners[i][c];
758 }
759 }
760 double eigen[3];
761 Vector3 axis [3];
762 covar.EigenSolveSymmetric(eigen, axis);
763 Vector3::GenerateOrthonormalBasis(axis[0], axis[1], axis[2]);
764 for (int d = 0; d < 3; ++d) {
765 out._xaxis[d] = axis[0][d];
766 out._yaxis[d] = axis[1][d];
767 out._zaxis[d] = axis[2][d];
768 }
769 // Set output size
770 out._x = out._y = out._z = 1;
771 AdjustFieldOfView(attr, out);
772 }
773 return out;
774 }
775
776 // -----------------------------------------------------------------------------
OverallFieldOfView(const Array<ImageAttributes> & attr)777 ImageAttributes OverallFieldOfView(const Array<ImageAttributes> &attr)
778 {
779 ImageAttributes out;
780 const int N = static_cast<int>(attr.size());
781 if (N == 0) return out;
782 if (N == 1) return attr[0];
783 // Set output resolution
784 out._dx = out._dy = out._dz = .0;
785 for (int n = 0; n < N; ++n) {
786 out._dx += attr[n]._dx;
787 out._dy += attr[n]._dy;
788 out._dz += attr[n]._dz;
789 }
790 out._dx /= N, out._dy /= N, out._dz /= N;
791 // Set output origin
792 out._xorigin = out._yorigin = out._zorigin = .0;
793 for (int n = 0; n < N; ++n) {
794 out._xorigin += attr[n]._xorigin;
795 out._yorigin += attr[n]._yorigin;
796 out._zorigin += attr[n]._zorigin;
797 }
798 out._xorigin /= N, out._yorigin /= N, out._zorigin /= N;
799 // Set output orientation
800 bool have_same_orientation = true;
801 const Matrix R = attr[0].GetLatticeToWorldOrientation();
802 for (int n = 1; n < N; ++n) {
803 if (attr[n].GetLatticeToWorldOrientation() != R) {
804 have_same_orientation = false;
805 break;
806 }
807 }
808 if (have_same_orientation) {
809 for (int d = 0; d < 3; ++d) {
810 out._xaxis[d] = attr[0]._xaxis[d];
811 out._yaxis[d] = attr[0]._yaxis[d];
812 out._zaxis[d] = attr[0]._zaxis[d];
813 }
814 } else {
815 const int ncorners = 8 * N;
816 double **corners = Allocate<double>(3, ncorners);
817 for (int n = 0; n < N; ++n) {
818 GetCorners(attr[n], corners + 8 * n);
819 }
820 Matrix3x3 covar(.0);
821 for (int i = 0; i < ncorners; ++i) {
822 corners[i][0] -= out._xorigin;
823 corners[i][1] -= out._yorigin;
824 corners[i][2] -= out._zorigin;
825 for (int r = 0; r < 3; ++r)
826 for (int c = 0; c < 3; ++c) {
827 covar[r][c] += corners[i][r] * corners[i][c];
828 }
829 }
830 Deallocate(corners);
831 double eigen[3];
832 Vector3 axis [3];
833 covar.EigenSolveSymmetric(eigen, axis);
834 Vector3::GenerateOrthonormalBasis(axis[0], axis[1], axis[2]);
835 for (int d = 0; d < 3; ++d) {
836 out._xaxis[d] = axis[0][d];
837 out._yaxis[d] = axis[1][d];
838 out._zaxis[d] = axis[2][d];
839 }
840 }
841 // Set output size
842 out._x = out._y = out._z = 1;
843 for (int n = 0; n < N; ++n) {
844 AdjustFieldOfView(attr[n], out);
845 }
846 return out;
847 }
848
849
850 } // namespace mirtk
851