orthanc-stone: Framework/Toolbox/ShearWarpProjectiveTransform.cpp annotate

annotate Framework/Toolbox/ShearWarpProjectiveTransform.cpp @ 199:dabe9982fca3

shear-warp on int16_t

author	Sebastien Jodogne <s.jodogne@gmail.com>
date	Thu, 22 Mar 2018 12:53:47 +0100
parents	4abddd083374
children	e9c7a78a3e77

rev	line source
191 46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	1 /**
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	2 * Stone of Orthanc
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	4 * Department, University Hospital of Liege, Belgium
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	5 * Copyright (C) 2017-2018 Osimis S.A., Belgium
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	6 *
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	7 * This program is free software: you can redistribute it and/or
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	8 * modify it under the terms of the GNU Affero General Public License
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	9 * as published by the Free Software Foundation, either version 3 of
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	10 * the License, or (at your option) any later version.
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	11 *
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	12 * This program is distributed in the hope that it will be useful, but
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	13 * WITHOUT ANY WARRANTY; without even the implied warranty of
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	15 * Affero General Public License for more details.
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	16 *
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	17 * You should have received a copy of the GNU Affero General Public License
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	19 **/
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	20
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	21
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	22 #include "ShearWarpProjectiveTransform.h"
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	23
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	24 #include "ImageGeometry.h"
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	25 #include "Extent2D.h"
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	26 #include "FiniteProjectiveCamera.h"
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	27 #include "GeometryToolbox.h"
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	28
193 4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	29 #include <Core/Images/PixelTraits.h>
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	30 #include <Core/Images/ImageProcessing.h>
191 46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	31 #include <Core/OrthancException.h>
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	32 #include <Core/Logging.h>
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	33
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	34 #include <boost/numeric/ublas/matrix_proxy.hpp>
193 4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	35 #include <boost/math/special_functions/round.hpp>
191 46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	36 #include <cassert>
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	37
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	38
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	39 namespace OrthancStone
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	40 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	41 static bool IsValidShear(const Matrix& M_shear)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	42 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	43 return (LinearAlgebra::IsCloseToZero(M_shear(0, 1)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	44 LinearAlgebra::IsCloseToZero(M_shear(1, 0)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	45 LinearAlgebra::IsCloseToZero(M_shear(2, 0)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	46 LinearAlgebra::IsCloseToZero(M_shear(2, 1)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	47 LinearAlgebra::IsNear(1.0, M_shear(2, 2)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	48 LinearAlgebra::IsCloseToZero(M_shear(2, 3)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	49 LinearAlgebra::IsCloseToZero(M_shear(3, 0)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	50 LinearAlgebra::IsCloseToZero(M_shear(3, 1)) &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	51 LinearAlgebra::IsNear(1.0, M_shear(3, 3)));
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	52 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	53
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	54
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	55 static void ComputeShearParameters(double& scaling,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	56 double& offsetX,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	57 double& offsetY,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	58 const Matrix& shear,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	59 double z)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	60 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	61 // Check out: ../../Resources/Computations/ComputeShearParameters.py
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	62
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	63 if (!LinearAlgebra::IsShearMatrix(shear))
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	64 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	65 LOG(ERROR) << "Not a valid shear matrix";
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	66 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	67 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	68
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	69 scaling = 1.0 / (shear(3,2) * z + 1.0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	70 offsetX = shear(0,2) * z * scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	71 offsetY = shear(1,2) * z * scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	72 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	73
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	74
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	75 ShearWarpProjectiveTransform::
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	76 ShearWarpProjectiveTransform(const Matrix& M_view,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	77 //const Matrix& P, // Permutation applied to the volume
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	78 unsigned int volumeWidth,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	79 unsigned int volumeHeight,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	80 unsigned int volumeDepth,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	81 double pixelSpacingX,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	82 double pixelSpacingY,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	83 unsigned int imageWidth,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	84 unsigned int imageHeight)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	85 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	86 eye_o.resize(4);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	87
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	88 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	89 // Find back the camera center given the "M_view" matrix
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	90 const double m11 = M_view(0, 0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	91 const double m12 = M_view(0, 1);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	92 const double m13 = M_view(0, 2);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	93 const double m14 = M_view(0, 3);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	94 const double m21 = M_view(1, 0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	95 const double m22 = M_view(1, 1);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	96 const double m23 = M_view(1, 2);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	97 const double m24 = M_view(1, 3);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	98 const double m41 = M_view(3, 0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	99 const double m42 = M_view(3, 1);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	100 const double m43 = M_view(3, 2);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	101 const double m44 = M_view(3, 3);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	102
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	103 // Equations (A.8) to (A.11) on page 203. Also check out
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	104 // "Finding the camera center" in "Multiple View Geometry in
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	105 // Computer Vision - 2nd edition", page 163.
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	106 const double vx[9] = { m12, m13, m14, m22, m23, m24, m42, m43, m44 };
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	107 const double vy[9] = { m11, m13, m14, m21, m23, m24, m41, m43, m44 };
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	108 const double vz[9] = { m11, m12, m14, m21, m22, m24, m41, m42, m44 };
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	109 const double vw[9] = { m11, m12, m13, m21, m22, m23, m41, m42, m43 };
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	110
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	111 Matrix m;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	112
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	113 LinearAlgebra::FillMatrix(m, 3, 3, vx);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	114 eye_o[0] = -LinearAlgebra::ComputeDeterminant(m);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	115
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	116 LinearAlgebra::FillMatrix(m, 3, 3, vy);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	117 eye_o[1] = LinearAlgebra::ComputeDeterminant(m);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	118
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	119 LinearAlgebra::FillMatrix(m, 3, 3, vz);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	120 eye_o[2] = -LinearAlgebra::ComputeDeterminant(m);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	121
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	122 LinearAlgebra::FillMatrix(m, 3, 3, vw);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	123 eye_o[3] = LinearAlgebra::ComputeDeterminant(m);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	124
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	125 if (LinearAlgebra::IsCloseToZero(eye_o[3]))
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	126 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	127 LOG(ERROR) << "The shear-warp projective transform is not applicable to affine cameras";
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	128 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	129 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	130 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	131
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	132 #if 0
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	133 // Assume "T_shift = I" (the eye does not lie on plane k = 0)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	134 const Matrix T_shift = LinearAlgebra::IdentityMatrix(4);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	135
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	136 // Equation (A.13) on page 204, given that the inverse of a
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	137 // permutation matrix is its transpose (TODO CHECK). If no T_shift
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	138 // or permutation P is applied, M'_view == M_view
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	139 const Matrix MM_view = LinearAlgebra::Product(
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	140 M_view,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	141 LinearAlgebra::Transpose(P),
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	142 LinearAlgebra::InvertScalingTranslationMatrix(T_shift));
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	143 #else
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	144 // This is a shortcut, as we take "T_shift = I" and "P = I"
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	145 const Matrix MM_view = M_view;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	146 #endif
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	147
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	148 // Equation (A.14) on page 207
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	149 Matrix MM_shear = LinearAlgebra::IdentityMatrix(4);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	150 MM_shear(0, 2) = -eye_o[0] / eye_o[2];
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	151 MM_shear(1, 2) = -eye_o[1] / eye_o[2];
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	152 MM_shear(3, 2) = -eye_o[3] / eye_o[2];
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	153
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	154
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	155 // Compute the extent of the intermediate image
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	156 Extent2D extent;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	157 double maxScaling = 1;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	158
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	159 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	160 // Compute the shearing factors of the two extreme planes of the
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	161 // volume (z=0 and z=volumeDepth)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	162 double scaling, offsetX, offsetY;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	163 ComputeShearParameters(scaling, offsetX, offsetY, MM_shear, 0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	164
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	165 if (scaling > 0)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	166 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	167 extent.AddPoint(offsetX, offsetY);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	168 extent.AddPoint(offsetX + static_cast<double>(volumeWidth) * scaling,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	169 offsetY + static_cast<double>(volumeHeight) * scaling);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	170
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	171 if (scaling > maxScaling)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	172 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	173 maxScaling = scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	174 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	175 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	176
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	177 ComputeShearParameters(scaling, offsetX, offsetY, MM_shear, volumeDepth);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	178
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	179 if (scaling > 0)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	180 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	181 extent.AddPoint(offsetX, offsetY);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	182 extent.AddPoint(offsetX + static_cast<double>(volumeWidth) * scaling,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	183 offsetY + static_cast<double>(volumeHeight) * scaling);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	184
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	185 if (scaling > maxScaling)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	186 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	187 maxScaling = scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	188 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	189 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	190 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	191
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	192 if (LinearAlgebra::IsCloseToZero(extent.GetWidth()) \|\|
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	193 LinearAlgebra::IsCloseToZero(extent.GetHeight()))
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	194 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	195 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	196 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	197
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	198 intermediateWidth_ = std::ceil(extent.GetWidth() / maxScaling);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	199 intermediateHeight_ = std::ceil(extent.GetHeight() / maxScaling);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	200
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	201 // This is the product "T * S" in Equation (A.16) on page 209
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	202 Matrix TS = LinearAlgebra::Product(
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	203 GeometryToolbox::CreateTranslationMatrix(static_cast<double>(intermediateWidth_) / 2.0,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	204 static_cast<double>(intermediateHeight_) / 2.0, 0),
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	205 GeometryToolbox::CreateScalingMatrix(1.0 / maxScaling, 1.0 / maxScaling, 1),
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	206 GeometryToolbox::CreateTranslationMatrix(-extent.GetCenterX(), -extent.GetCenterY(), 0));
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	207
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	208 // This is Equation (A.16) on page 209. WARNING: There is an
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	209 // error in Lacroute's thesis: "inv(MM_shear)" is used instead
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	210 // of "MM_shear".
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	211 M_shear = LinearAlgebra::Product(TS, MM_shear);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	212
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	213 if (!IsValidShear(M_shear))
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	214 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	215 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	216 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	217
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	218 // This is Equation (A.17) on page 209
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	219 Matrix tmp;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	220 LinearAlgebra::InvertMatrix(tmp, M_shear);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	221 M_warp = LinearAlgebra::Product(MM_view, tmp);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	222
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	223 // Intrinsic parameters of the camera
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	224 k_ = LinearAlgebra::ZeroMatrix(3, 4);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	225 k_(0, 0) = 1.0 / pixelSpacingX;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	226 k_(0, 3) = static_cast<double>(imageWidth) / 2.0;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	227 k_(1, 1) = 1.0 / pixelSpacingY;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	228 k_(1, 3) = static_cast<double>(imageHeight) / 2.0;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	229 k_(2, 3) = 1.0;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	230 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	231
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	232
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	233 FiniteProjectiveCamera *ShearWarpProjectiveTransform::CreateCamera() const
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	234 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	235 Matrix p = LinearAlgebra::Product(k_, M_warp, M_shear);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	236 return new FiniteProjectiveCamera(p);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	237 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	238
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	239
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	240 void ShearWarpProjectiveTransform::ComputeShearOnSlice(double& a11,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	241 double& b1,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	242 double& a22,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	243 double& b2,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	244 double& shearedZ,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	245 const double sourceZ)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	246 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	247 // Check out: ../../Resources/Computations/ComputeShearOnSlice.py
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	248 assert(IsValidShear(M_shear));
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	249
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	250 const double s11 = M_shear(0, 0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	251 const double s13 = M_shear(0, 2);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	252 const double s14 = M_shear(0, 3);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	253 const double s22 = M_shear(1, 1);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	254 const double s23 = M_shear(1, 2);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	255 const double s24 = M_shear(1, 3);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	256 const double s43 = M_shear(3, 2);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	257
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	258 double scaling = 1.0 / (s43 * sourceZ + 1.0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	259 shearedZ = sourceZ * scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	260
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	261 a11 = s11 * scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	262 a22 = s22 * scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	263
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	264 b1 = (s13 * sourceZ + s14) * scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	265 b2 = (s23 * sourceZ + s24) * scaling;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	266 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	267
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	268
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	269 Matrix ShearWarpProjectiveTransform::CalibrateView(const Vector& camera,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	270 const Vector& principalPoint,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	271 double angle)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	272 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	273 if (camera.size() != 3 \|\|
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	274 principalPoint.size() != 3)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	275 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	276 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	277 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	278
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	279 const double sid = boost::numeric::ublas::norm_2(camera - principalPoint);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	280
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	281 Matrix a;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	282 GeometryToolbox::AlignVectorsWithRotation(a, camera - principalPoint,
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	283 LinearAlgebra::CreateVector(0, 0, -1));
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	284
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	285 Matrix r = LinearAlgebra::Product(GeometryToolbox::CreateRotationMatrixAlongZ(angle), a);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	286
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	287 a = LinearAlgebra::ZeroMatrix(4, 4);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	288 boost::numeric::ublas::subrange(a, 0, 3, 0, 3) = r;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	289
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	290 const Vector v = LinearAlgebra::Product(r, -camera);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	291 a(0, 3) = v[0];
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	292 a(1, 3) = v[1];
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	293 a(2, 3) = v[2];
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	294 a(3, 3) = 1;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	295
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	296 Matrix perspective = LinearAlgebra::ZeroMatrix(4, 4);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	297 // https://stackoverflow.com/questions/5267866/calculation-of-a-perspective-transformation-matrix
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	298 perspective(0, 0) = sid;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	299 perspective(1, 1) = sid;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	300 perspective(2, 2) = sid;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	301 perspective(3, 2) = 1;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	302
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	303 Matrix M_view = LinearAlgebra::Product(perspective, a);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	304 assert(M_view.size1() == 4 &&
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	305 M_view.size2() == 4);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	306
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	307 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	308 // Sanity checks
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	309 Vector p1 = LinearAlgebra::CreateVector(camera[0], camera[1], camera[2], 1.0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	310 Vector p2 = LinearAlgebra::CreateVector(principalPoint[0], principalPoint[1], principalPoint[2], 1.0);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	311
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	312 Vector v1 = LinearAlgebra::Product(M_view, p1);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	313 Vector v2 = LinearAlgebra::Product(M_view, p2);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	314
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	315 if (!LinearAlgebra::IsCloseToZero(v1[3]) \|\| // Must be mapped to singularity (w=0)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	316 LinearAlgebra::IsCloseToZero(v2[3]))
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	317 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	318 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	319 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	320
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	321 // The principal point must be mapped to (0,0,z,1)
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	322 v2 /= v2[3];
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	323 if (!LinearAlgebra::IsCloseToZero(v2[0]) \|\|
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	324 !LinearAlgebra::IsCloseToZero(v2[1]))
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	325 {
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	326 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	327 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	328 }
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	329
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	330 return M_view;
46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	331 }
193 4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	332
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	333
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	334 template <Orthanc::PixelFormat SourceFormat,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	335 Orthanc::PixelFormat TargetFormat,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	336 bool MIP>
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	337 static void ApplyAxialInternal(Orthanc::ImageAccessor& target,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	338 float& maxValue,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	339 const Matrix& M_view,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	340 const ImageBuffer3D& source,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	341 double pixelSpacing,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	342 unsigned int countSlices,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	343 ImageInterpolation shearInterpolation,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	344 ImageInterpolation warpInterpolation)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	345 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	346 typedef Orthanc::PixelTraits<SourceFormat> SourceTraits;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	347 typedef Orthanc::PixelTraits<TargetFormat> TargetTraits;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	348
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	349 /**
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	350 * Step 1: Precompute some information.
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	351 **/
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	352
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	353 if (target.GetFormat() != TargetFormat \|\|
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	354 source.GetFormat() != SourceFormat \|\|
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	355 !std::numeric_limits<float>::is_iec559 \|\|
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	356 sizeof(float) != 4)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	357 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	358 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	359 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	360
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	361 if (countSlices > source.GetDepth())
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	362 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	363 countSlices = source.GetDepth();
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	364 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	365
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	366 if (countSlices == 0)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	367 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	368 maxValue = 0;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	369 Orthanc::ImageProcessing::Set(target, 0);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	370 return;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	371 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	372
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	373 LOG(INFO) << "Number of rendered slices: " << countSlices;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	374
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	375
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	376 /**
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	377 * Step 2: Extract the shear-warp transform corresponding to
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	378 * M_view.
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	379 **/
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	380
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	381 // Compute the "world" matrix that maps the source volume to the
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	382 // (0,0,0)->(1,1,1) unit cube
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	383 Vector origin = source.GetCoordinates(0, 0, 0);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	384 Vector ps = source.GetVoxelDimensions(VolumeProjection_Axial);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	385 Matrix world = LinearAlgebra::Product(
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	386 GeometryToolbox::CreateScalingMatrix(1.0 / ps[0], 1.0 / ps[1], 1.0 / ps[2]),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	387 GeometryToolbox::CreateTranslationMatrix(-origin[0], -origin[1], -origin[2]));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	388
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	389 Matrix worldInv;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	390 LinearAlgebra::InvertMatrix(worldInv, world);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	391
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	392 ShearWarpProjectiveTransform shearWarp(LinearAlgebra::Product(M_view, worldInv),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	393 /LinearAlgebra::IdentityMatrix(4),/
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	394 source.GetWidth(),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	395 source.GetHeight(),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	396 source.GetDepth(),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	397 pixelSpacing, pixelSpacing,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	398 target.GetWidth(), target.GetHeight());
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	399
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	400 const unsigned int intermediateWidth = shearWarp.GetIntermediateWidth();
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	401 const unsigned int intermediateHeight = shearWarp.GetIntermediateHeight();
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	402
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	403
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	404 /**
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	405 * Step 3: Apply the "shear" part of the transform to form the
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	406 * intermediate image. The sheared images are accumulated into the
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	407 * Float32 image "accumulator". The number of samples available
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	408 * for each pixel is stored in the "counter" image.
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	409 **/
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	410
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	411 std::auto_ptr<Orthanc::ImageAccessor> accumulator, counter, intermediate;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	412
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	413 accumulator.reset(new Orthanc::Image(Orthanc::PixelFormat_Float32,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	414 intermediateWidth, intermediateHeight, false));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	415 counter.reset(new Orthanc::Image(Orthanc::PixelFormat_Grayscale16,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	416 intermediateWidth, intermediateHeight, false));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	417 intermediate.reset(new Orthanc::Image(SourceFormat, intermediateWidth, intermediateHeight, false));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	418
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	419 Orthanc::ImageProcessing::Set(*accumulator, 0);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	420 Orthanc::ImageProcessing::Set(*counter, 0);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	421
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	422 // Loop around the slices of the volume
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	423 for (unsigned int i = 0; i <= countSlices; i++)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	424 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	425 // (3.a) Compute the shear for this specific slice
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	426 unsigned int z = static_cast<unsigned int>(
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	427 boost::math::iround(static_cast<double>(i) /
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	428 static_cast<double>(countSlices) *
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	429 static_cast<double>(source.GetDepth() - 1)));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	430
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	431 double a11, b1, a22, b2, vz;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	432 shearWarp.ComputeShearOnSlice(a11, b1, a22, b2, vz, static_cast<double>(z) + 0.5);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	433
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	434
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	435 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	436 // (3.b) Detect the "useful" portion of the intermediate image
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	437 // for this slice (i.e. the bounding box where the source
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	438 // slice is mapped to by the shear), so as to update "counter"
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	439 Matrix a = LinearAlgebra::ZeroMatrix(3, 3);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	440 a(0,0) = a11;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	441 a(0,2) = b1;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	442 a(1,1) = a22;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	443 a(1,2) = b2;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	444 a(2,2) = 1;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	445
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	446 unsigned int x1, y1, x2, y2;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	447 if (GetProjectiveTransformExtent(x1, y1, x2, y2, a,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	448 source.GetWidth(), source.GetHeight(),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	449 intermediateWidth, intermediateHeight))
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	450 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	451 for (unsigned int y = y1; y <= y2; y++)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	452 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	453 uint16_t* p = reinterpret_cast<uint16_t*>(counter->GetRow(y)) + x1;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	454 for (unsigned int x = x1; x <= x2; x++, p++)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	455 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	456 if (MIP)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	457 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	458 // TODO - In the case of MIP, "counter" could be
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	459 // reduced to "PixelFormat_Grayscale8" to reduce
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	460 // memory usage
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	461 *p = 1;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	462 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	463 else
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	464 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	465 *p += 1;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	466 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	467 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	468 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	469 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	470 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	471
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	472
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	473 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	474 // (3.c) Shear the source slice into a temporary image
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	475 ImageBuffer3D::SliceReader reader(source, VolumeProjection_Axial, z);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	476 ApplyAffineTransform(*intermediate, reader.GetAccessor(),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	477 a11, 0, b1,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	478 0, a22, b2,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	479 shearInterpolation);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	480 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	481
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	482
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	483 for (unsigned int y = 0; y < intermediateHeight; y++)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	484 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	485 // (3.d) Accumulate the pixels of the sheared image into "accumulator"
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	486 const typename SourceTraits::PixelType* p =
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	487 reinterpret_cast<const typename SourceTraits::PixelType*>(intermediate->GetConstRow(y));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	488
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	489 float* q = reinterpret_cast<float*>(accumulator->GetRow(y));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	490
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	491 for (unsigned int x = 0; x < intermediateWidth; x++)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	492 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	493 float pixel = SourceTraits::PixelToFloat(*p);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	494
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	495 if (MIP)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	496 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	497 // Get maximum for MIP
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	498 if (*q < pixel)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	499 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	500 *q = pixel;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	501 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	502 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	503 else
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	504 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	505 *q += pixel;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	506 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	507
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	508 p++;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	509 q++;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	510 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	511 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	512 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	513
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	514
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	515 /**
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	516 * Step 4: The intermediate image (that will be transformed by the
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	517 * "warp") is now available as an accumulator image together with
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	518 * a counter image. "Flatten" these two images into one.
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	519 **/
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	520
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	521 intermediate.reset(new Orthanc::Image
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	522 (TargetFormat, intermediateWidth, intermediateHeight, false));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	523
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	524 maxValue = 0;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	525
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	526 for (unsigned int y = 0; y < intermediateHeight; y++)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	527 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	528 const float qacc = reinterpret_cast<const float>(accumulator->GetConstRow(y));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	529 const uint16_t qcount = reinterpret_cast<const uint16_t>(counter->GetConstRow(y));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	530 typename TargetTraits::PixelType *p =
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	531 reinterpret_cast<typename TargetTraits::PixelType*>(intermediate->GetRow(y));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	532
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	533 for (unsigned int x = 0; x < intermediateWidth; x++)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	534 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	535 if (*qcount == 0)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	536 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	537 TargetTraits::SetZero(*p);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	538 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	539 else
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	540 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	541 p = qacc / static_cast<float>(*qcount);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	542
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	543 if (*p > maxValue)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	544 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	545 maxValue = *p;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	546 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	547 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	548
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	549 p++;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	550 qacc++;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	551 qcount++;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	552 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	553 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	554
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	555 // We don't need the accumulator images anymore
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	556 accumulator.reset(NULL);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	557 counter.reset(NULL);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	558
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	559
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	560 /**
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	561 * Step 6: Apply the "warp" part of the transform to map the
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	562 * intermediate image to the final image.
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	563 **/
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	564
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	565 Matrix warp;
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	566
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	567 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	568 // (5.a) Compute the "warp" matrix by removing the 3rd row and
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	569 // 3rd column from the GetWarp() matrix
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	570 // Check out: ../../Resources/Computations/ComputeWarp.py
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	571
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	572 Matrix fullWarp = LinearAlgebra::Product
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	573 (shearWarp.GetIntrinsicParameters(), shearWarp.GetWarp());
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	574
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	575 const double v[] = {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	576 fullWarp(0,0), fullWarp(0,1), fullWarp(0,3),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	577 fullWarp(1,0), fullWarp(1,1), fullWarp(1,3),
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	578 fullWarp(2,0), fullWarp(2,1), fullWarp(2,3)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	579 };
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	580
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	581 LinearAlgebra::FillMatrix(warp, 3, 3, v);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	582 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	583
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	584 // (5.b) Apply the projective transform to the image
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	585 ApplyProjectiveTransform(target, *intermediate, warp, warpInterpolation);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	586 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	587
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	588
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	589 template <Orthanc::PixelFormat SourceFormat,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	590 Orthanc::PixelFormat TargetFormat>
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	591 static void ApplyAxialInternal2(Orthanc::ImageAccessor& target,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	592 float& maxValue,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	593 const Matrix& M_view,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	594 const ImageBuffer3D& source,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	595 bool mip,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	596 double pixelSpacing,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	597 unsigned int countSlices,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	598 ImageInterpolation shearInterpolation,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	599 ImageInterpolation warpInterpolation)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	600 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	601 if (mip)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	602 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	603 ApplyAxialInternal<SourceFormat, TargetFormat, true>
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	604 (target, maxValue, M_view, source, pixelSpacing,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	605 countSlices, shearInterpolation, warpInterpolation);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	606 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	607 else
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	608 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	609 ApplyAxialInternal<SourceFormat, TargetFormat, false>
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	610 (target, maxValue, M_view, source, pixelSpacing,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	611 countSlices, shearInterpolation, warpInterpolation);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	612 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	613 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	614
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	615
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	616 Orthanc::ImageAccessor*
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	617 ShearWarpProjectiveTransform::ApplyAxial(float& maxValue,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	618 const Matrix& M_view,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	619 const ImageBuffer3D& source,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	620 Orthanc::PixelFormat targetFormat,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	621 unsigned int targetWidth,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	622 unsigned int targetHeight,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	623 bool mip,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	624 double pixelSpacing,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	625 unsigned int countSlices,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	626 ImageInterpolation shearInterpolation,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	627 ImageInterpolation warpInterpolation)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	628 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	629 std::auto_ptr<Orthanc::ImageAccessor> target
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	630 (new Orthanc::Image(targetFormat, targetWidth, targetHeight, false));
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	631
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	632 if (source.GetFormat() == Orthanc::PixelFormat_Grayscale16 &&
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	633 targetFormat == Orthanc::PixelFormat_Grayscale16)
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	634 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	635 ApplyAxialInternal2<Orthanc::PixelFormat_Grayscale16,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	636 Orthanc::PixelFormat_Grayscale16>
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	637 (*target, maxValue, M_view, source, mip, pixelSpacing,
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	638 countSlices, shearInterpolation, warpInterpolation);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	639 }
199 dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	640 else if (source.GetFormat() == Orthanc::PixelFormat_SignedGrayscale16 &&
dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	641 targetFormat == Orthanc::PixelFormat_SignedGrayscale16)
dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	642 {
dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	643 ApplyAxialInternal2<Orthanc::PixelFormat_SignedGrayscale16,
dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	644 Orthanc::PixelFormat_SignedGrayscale16>
dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	645 (*target, maxValue, M_view, source, mip, pixelSpacing,
dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	646 countSlices, shearInterpolation, warpInterpolation);
dabe9982fca3 shear-warp on int16_t Sebastien Jodogne <s.jodogne@gmail.com> parents: 193 diff changeset	647 }
193 4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	648 else
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	649 {
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	650 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	651 }
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	652
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	653 return target.release();
4abddd083374 ShearWarpProjectiveTransform::ApplyAxial() Sebastien Jodogne <s.jodogne@gmail.com> parents: 191 diff changeset	654 }
191 46cb2eedc2e0 ShearWarpProjectiveTransform Sebastien Jodogne <s.jodogne@gmail.com> parents: diff changeset	655 }

Mercurial > hg > orthanc-stone

annotate Framework/Toolbox/ShearWarpProjectiveTransform.cpp @ 199:dabe9982fca3