Developer documentation
Version 3.0.3-105-gd3941f44
mean_squared.h
Go to the documentation of this file.
1/* Copyright (c) 2008-2022 the MRtrix3 contributors.
2 *
3 * This Source Code Form is subject to the terms of the Mozilla Public
4 * License, v. 2.0. If a copy of the MPL was not distributed with this
5 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
6 *
7 * Covered Software is provided under this License on an "as is"
8 * basis, without warranty of any kind, either expressed, implied, or
9 * statutory, including, without limitation, warranties that the
10 * Covered Software is free of defects, merchantable, fit for a
11 * particular purpose or non-infringing.
12 * See the Mozilla Public License v. 2.0 for more details.
13 *
14 * For more details, see http://www.mrtrix.org/.
15 */
16
17#ifndef __registration_metric_mean_squared_h__
18#define __registration_metric_mean_squared_h__
19
20#include "math/math.h"
21#include "registration/metric/linear_base.h"
22#include "transform.h"
23
24namespace MR
25{
26 namespace Registration
27 {
28 namespace Metric
29 {
30
31 class MeanSquared : public LinearBase { MEMALIGN(MeanSquared)
32 public:
33 template <class Params>
34 default_type operator() (Params& params,
35 const Eigen::Vector3d& im1_point,
36 const Eigen::Vector3d& im2_point,
37 const Eigen::Vector3d& midway_point,
38 Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
39
40 assert (!this->weighted && "FIXME: set_weights not implemented for 3D metric");
41
42 typename Params::Im1ValueType im1_value;
43 typename Params::Im2ValueType im2_value;
44 Eigen::Matrix<typename Params::Im1ValueType, 1, 3> im1_grad;
45 Eigen::Matrix<typename Params::Im2ValueType, 1, 3> im2_grad;
46
47 params.im1_image_interp->value_and_gradient_wrt_scanner (im1_value, im1_grad);
48 if (std::isnan (default_type (im1_value)))
49 return 0.0;
50 params.im2_image_interp->value_and_gradient_wrt_scanner (im2_value, im2_grad);
51 if (std::isnan (default_type (im2_value)))
52 return 0.0;
53
54 default_type diff = (default_type) im1_value - (default_type) im2_value;
55
56 const auto jacobian_vec = params.transformation.get_jacobian_vector_wrt_params (midway_point);
57 const Eigen::Vector3d g = diff * (im1_grad + im2_grad);
58 gradient.segment<4>(0) += g(0) * jacobian_vec;
59 gradient.segment<4>(4) += g(1) * jacobian_vec;
60 gradient.segment<4>(8) += g(2) * jacobian_vec;
61
62 return diff * diff;
63 }
64 };
65
66 class MeanSquaredNoGradient : public LinearBase { MEMALIGN(MeanSquaredNoGradient)
67 public:
68 template <class Params>
69 default_type operator() (Params& params,
70 const Eigen::Vector3d& im1_point,
71 const Eigen::Vector3d& im2_point,
72 const Eigen::Vector3d& midway_point,
73 Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
74
75 assert (!this->weighted && "FIXME: set_weights not implemented for 3D metric");
76
77 typename Params::Im1ValueType im1_value;
78 typename Params::Im2ValueType im2_value;
79
80 im1_value = params.im1_image_interp->value ();
81 if (std::isnan (default_type (im1_value)))
82 return 0.0;
83 im2_value = params.im2_image_interp->value ();
84 if (std::isnan (default_type (im2_value)))
85 return 0.0;
86
87 default_type diff = (default_type) im1_value - (default_type) im2_value;
88
89 return diff * diff;
90 }
91 };
92
93 class MeanSquaredNonSymmetric : public LinearBase { MEMALIGN(MeanSquaredNonSymmetric)
94 public:
98 using is_asymmetric_type = int;
99
100 template <class Params>
101 default_type operator() (Params& params,
102 const Eigen::Vector3d& im1_point,
103 const Eigen::Vector3d& im2_point,
104 const Eigen::Vector3d& midway_point,
105 Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
106
107 assert (!this->weighted && "FIXME: set_weights not implemented for 3D metric");
108
109 typename Params::Im1ValueType im1_value;
110 default_type im2_value;
111 Eigen::Matrix<typename Params::Im1ValueType, 1, 3> im1_grad;
112
113 params.im1_image_interp->value_and_gradient_wrt_scanner (im1_value, im1_grad);
114 if (std::isnan (default_type (im1_value)))
115 return 0.0;
116
117 assert(params.im2_image.index(0) == std::round((MR::Transform(params.midway_image).voxel2scanner.inverse() * midway_point)[0]));
118 assert(params.im2_image.index(1) == std::round((MR::Transform(params.midway_image).voxel2scanner.inverse() * midway_point)[1]));
119 assert(params.im2_image.index(2) == std::round((MR::Transform(params.midway_image).voxel2scanner.inverse() * midway_point)[2]));
120 im2_value = params.im2_image.value();
121
122 if (std::isnan (im2_value))
123 return 0.0;
124
125 default_type diff = (default_type) im1_value - im2_value;
126
127 const auto jacobian_vec = params.transformation.get_jacobian_vector_wrt_params (im2_point);
128 const Eigen::Vector3d g = 2.0 * diff * im1_grad;
129 gradient.segment<4>(0) += g(0) * jacobian_vec;
130 gradient.segment<4>(4) += g(1) * jacobian_vec;
131 gradient.segment<4>(8) += g(2) * jacobian_vec;
132
133 return diff * diff;
134 }
135 };
136
137 template <class Im1Type, class Im2Type>
138 class MeanSquared4D : public LinearBase { MEMALIGN(MeanSquared4D<Im1Type,Im2Type>)
139 public:
140 template <class Params>
141 default_type operator() (Params& params,
142 const Eigen::Vector3d& im1_point,
143 const Eigen::Vector3d& im2_point,
144 const Eigen::Vector3d& midway_point,
145 Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
146
147 const ssize_t volumes = params.im1_image_interp->size(3);
148
149 if (im1_values.rows() != volumes) {
150 im1_values.resize (volumes);
151 im1_grad.resize (volumes, 3);
152 diff_values.resize (volumes);
153 im2_values.resize (volumes);
154 im2_grad.resize (volumes, 3);
155 }
156
157 params.im1_image_interp->value_and_gradient_row_wrt_scanner (im1_values, im1_grad);
158 if (im1_values.hasNaN())
159 return 0.0;
160
161 params.im2_image_interp->value_and_gradient_row_wrt_scanner (im2_values, im2_grad);
162 if (im2_values.hasNaN())
163 return 0.0;
164
165 const Eigen::Matrix<default_type, 4, 1> jacobian_vec = params.transformation.get_jacobian_vector_wrt_params (midway_point);
166 diff_values = im1_values - im2_values;
167
168 if (this->weighted)
169 diff_values.array() *= this->mc_weights.array();
170
171 for (ssize_t i = 0; i < volumes; ++i) {
172 const Eigen::Vector3d g = diff_values[i] * (im1_grad.row(i) + im2_grad.row(i));
173 gradient.segment<4>(0) += g(0) * jacobian_vec;
174 gradient.segment<4>(4) += g(1) * jacobian_vec;
175 gradient.segment<4>(8) += g(2) * jacobian_vec;
176 }
177
178 return diff_values.squaredNorm() / volumes;
179 }
180 private:
181 Eigen::Matrix<typename Im1Type::value_type, Eigen::Dynamic, 1> im1_values;
182 Eigen::Matrix<typename Im2Type::value_type, Eigen::Dynamic, 1> im2_values;
183 Eigen::Matrix<typename Im1Type::value_type, Eigen::Dynamic, 3> im1_grad;
184 Eigen::Matrix<typename Im2Type::value_type, Eigen::Dynamic, 3> im2_grad;
185 Eigen::VectorXd diff_values;
186 };
187
188 template <class Im1Type, class Im2Type>
189 class MeanSquaredNoGradient4D : public LinearBase { MEMALIGN(MeanSquaredNoGradient4D<Im1Type,Im2Type>)
190 public:
191 MeanSquaredNoGradient4D ( ) {}
192
193 template <class Params>
194 default_type operator() (Params& params,
195 const Eigen::Vector3d& im1_point,
196 const Eigen::Vector3d& im2_point,
197 const Eigen::Vector3d& midway_point,
198 Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
199 const ssize_t volumes = params.im1_image_interp->size(3);
200
201 if (im1_values.rows() != volumes) {
202 im1_values.resize (volumes);
203 im2_values.resize (volumes);
204 diff_values.resize (volumes);
205 }
206
207 params.im1_image_interp->row (im1_values);
208 if (im1_values.hasNaN())
209 return 0.0;
210
211 params.im2_image_interp->row (im2_values);
212 if (im2_values.hasNaN())
213 return 0.0;
214
215 diff_values = im1_values - im2_values;
216 if (this->weighted)
217 diff_values.array() *= this->mc_weights.array();
218
219 return diff_values.squaredNorm() / volumes;
220 }
221 private:
222 Eigen::Matrix<typename Im1Type::value_type, Eigen::Dynamic, 1> im1_values;
223 Eigen::Matrix<typename Im2Type::value_type, Eigen::Dynamic, 1> im2_values;
224 Eigen::VectorXd diff_values;
225 };
226
227 template <class Im1Type, class Im2Type>
228 class MeanSquared4DNonSymmetric : public LinearBase { MEMALIGN(MeanSquared4DNonSymmetric<Im1Type,Im2Type>)
229 public:
230 template <class Params>
231 default_type operator() (Params& params,
232 const Eigen::Vector3d& im1_point,
233 const Eigen::Vector3d& im2_point,
234 const Eigen::Vector3d& midway_point,
235 Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
236
237 const ssize_t volumes = params.im1_image_interp->size(3);
238
239 if (im1_values.rows() != volumes) {
240 im1_values.resize (volumes);
241 im1_grad.resize (volumes, 3);
242 diff_values.resize (volumes);
243 im2_values.resize (volumes);
244 im2_grad.resize (volumes, 3);
245 }
246
247 params.im1_image_interp->value_and_gradient_row_wrt_scanner (im1_values, im1_grad);
248 if (im1_values.hasNaN())
249 return 0.0;
250
251 params.im2_image_interp->value_and_gradient_row_wrt_scanner (im2_values, im2_grad);
252 if (im2_values.hasNaN())
253 return 0.0;
254
255 assert(params.im2_image.index(0) == std::round((MR::Transform(params.midway_image).voxel2scanner.inverse() * midway_point)[0]));
256 assert(params.im2_image.index(1) == std::round((MR::Transform(params.midway_image).voxel2scanner.inverse() * midway_point)[1]));
257 assert(params.im2_image.index(2) == std::round((MR::Transform(params.midway_image).voxel2scanner.inverse() * midway_point)[2]));
258 im2_values = params.im2_image.row(3);
259
260 if (im2_values.hasNaN())
261 return 0.0;
262
263 const Eigen::Matrix<default_type, 4, 1> jacobian_vec = params.transformation.get_jacobian_vector_wrt_params (im2_point);
264 diff_values = im1_values - im2_values;
265
266 if (this->weighted)
267 diff_values.array() *= this->mc_weights.array();
268
269 for (ssize_t i = 0; i < volumes; ++i) {
270 const Eigen::Vector3d g = 2.0 * diff_values[i] * im1_grad.row(i);
271 gradient.segment<4>(0) += g(0) * jacobian_vec;
272 gradient.segment<4>(4) += g(1) * jacobian_vec;
273 gradient.segment<4>(8) += g(2) * jacobian_vec;
274 }
275
276 return diff_values.squaredNorm() / volumes;
277 }
278 private:
279 Eigen::Matrix<typename Im1Type::value_type, Eigen::Dynamic, 1> im1_values;
280 Eigen::Matrix<typename Im2Type::value_type, Eigen::Dynamic, 1> im2_values;
281 Eigen::Matrix<typename Im1Type::value_type, Eigen::Dynamic, 3> im1_grad;
282 Eigen::Matrix<typename Im2Type::value_type, Eigen::Dynamic, 3> im2_grad;
283 Eigen::VectorXd diff_values;
284 };
285
286 template <class Im1Type, class Im2Type>
287 class MeanSquaredVectorNoGradient4D : public LinearBase { MEMALIGN(MeanSquaredVectorNoGradient4D<Im1Type,Im2Type>)
288 public:
289 MeanSquaredVectorNoGradient4D () = delete;
290 MeanSquaredVectorNoGradient4D ( const Im1Type im1, const Im2Type im2 ):
291 volumes ( (im1.ndim() == 3) ? 1 : im1.size(3) ) {
292 assert (im1.ndim() == 4);
293 assert (im1.ndim() == im2.ndim());
294 assert (im1.size(3) == im2.size(3));
295 diff_values.resize (volumes);
296 im1_values.resize (volumes);
297 im2_values.resize (volumes);
298 }
299
300 //type_trait to indicate return type of operator()
301 using is_vector_type = int;
302
303 template <class Params>
304 Eigen::Matrix<default_type, Eigen::Dynamic, 1> operator() (Params& params,
305 const Eigen::Vector3d& im1_point,
306 const Eigen::Vector3d& im2_point,
307 const Eigen::Vector3d& midway_point,
308 Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
309
310 assert (volumes == params.im1_image_interp->size(3));
311
312 im1_values = params.im1_image_interp->row (3);
313 if (im1_values.hasNaN()) {
314 diff_values.setZero();
315 return diff_values;
316 }
317
318 im2_values = params.im2_image_interp->row (3);
319 if (im2_values.hasNaN()) {
320 diff_values.setZero();
321 return diff_values;
322 }
323
324 diff_values = im1_values - im2_values;
325 if (this->weighted)
326 diff_values.array() *= this->mc_weights.array();
327
328 return diff_values.array().square();
329 }
330 private:
331 ssize_t volumes;
332 Eigen::Matrix<typename Im1Type::value_type, Eigen::Dynamic, 1> im1_values;
333 Eigen::Matrix<typename Im2Type::value_type, Eigen::Dynamic, 1> im2_values;
334 Eigen::VectorXd diff_values;
335 };
336 }
337 }
338}
339#endif
MR::Transform::voxel2scanner
const transform_type voxel2scanner
Definition: transform.h:43
MR::Math::round
constexpr I round(const T x)
Definition: math.h:64
linear_base.h
MR
Definition: base.h:24
MR::default_type
double default_type
the default type used throughout MRtrix
Definition: types.h:228
MEMALIGN
#define MEMALIGN(...)
Definition: types.h:185