Developer documentation
Version 3.0.3-105-gd3941f44
rician.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 __math_rician_h__
18#define __math_rician_h__
19
20#include "math/math.h"
21#include "math/bessel.h"
22#include "math/vector.h"
23
24namespace MR
25{
26 namespace Math
27 {
28 namespace Rician
29 {
30
31 template <typename T> inline T lnP (const T measured, const T actual, const T one_over_noise_squared)
32 {
33 T nm = one_over_noise_squared * measured;
34 T s = abs (actual);
35 return 0.5 * one_over_noise_squared * pow2 (measured - s) - log (nm * Bessel::I0_scaled (nm * s));
36 }
37
38
39 template <typename T> inline T lnP (const T measured, T actual, const T one_over_noise_squared, T& dP_dactual)
40 {
41 assert (measured >= 0.0);
42 assert (one_over_noise_squared > 0.0);
43
44 actual = abs (actual);
45 T nm = one_over_noise_squared * measured;
46 T nms = nm * actual;
47 T F0 = Bessel::I0_scaled (nms);
48 T m_a = measured - actual;
49 T nm_a = one_over_noise_squared * m_a;
50 T F1_F0 = (Bessel::I1_scaled (nms) - F0) / F0;
51 assert (nms >= 0.0); // (nms < 0.0) F1_F0 = -F1_F0;
52 dP_dactual = -nm_a - nm * F1_F0;
53 return 0.5 * nm_a * m_a - log (nm * F0);
54 }
55
56
57
58
59 template <typename T> inline T lnP (const T measured, const T actual, const T one_over_noise_squared, T& dP_dactual, T& dP_dN)
60 {
61 assert (measured >= 0.0);
62 assert (one_over_noise_squared > 0.0);
63
64 bool actual_is_positive = (actual >= 0.0);
65 T actual_pos = abs (actual);
66 T nm = one_over_noise_squared * measured;
67 T nms = nm * actual_pos;
68 T F0 = Bessel::I0_scaled (nms);
69 T m_a = measured - actual_pos;
70 T nm_a = one_over_noise_squared * m_a;
71 T F1_F0 = (Bessel::I1_scaled (nms) - F0) / F0;
72 if (nms < 0.0) F1_F0 = -F1_F0;
73 dP_dactual = -nm_a - nm * F1_F0;
74 actual_pos *= measured * F1_F0;
75 dP_dN = 0.5 * pow2 (m_a) - 1.0/one_over_noise_squared + (actual_is_positive ? -actual_pos : actual_pos);
76 return 0.5 * nm_a * m_a - log (nm * F0);
77 }
78
79
80
81
82
83 template <typename T> inline T lnP (const int N, const T* measured, const T* actual, const T one_over_noise_squared, T* dP_dactual, T& dP_dN)
84 {
85 assert (one_over_noise_squared > 0.0);
86
87 T lnP = 0.0;
88 dP_dN = -T (N) / one_over_noise_squared;
89
90 for (int i = 0; i < N; i++) {
91 assert (measured[i] >= 0.0);
92
93 T actual_pos = abs (actual[i]);
94 T nm = one_over_noise_squared * measured[i];
95 T nms = nm * actual_pos;
96 T F0 = Bessel::I0_scaled (nms);
97 T m_a = measured[i] - actual_pos;
98 T nm_a = one_over_noise_squared * m_a;
99 T F1_F0 = (Bessel::I1_scaled (nms) - F0) / F0;
100 dP_dactual[i] = -nm_a - nm * F1_F0;
101 if (actual[i] < 0.0) dP_dactual[i] = -dP_dactual[i];
102 dP_dN += 0.5 * pow2 (m_a) - measured[i] * actual_pos * F1_F0;
103 lnP += 0.5 * nm_a * m_a - log (nm * F0);
104 }
105
106 return lnP;
107 }
108
109
110
111 template <typename T> inline T lnP (const Vector<T>& measured, const Vector<T>& actual, const T one_over_noise_squared, Vector<T>& dP_dactual)
112 {
113 assert (one_over_noise_squared > 0.0);
114 assert (measured.size() == actual.size());
115 assert (measured.size() == dP_dactual.size());
116
117 T lnP = 0.0;
118
119 for (size_t i = 0; i < measured.size(); i++) {
120 assert (measured[i] > 0.0);
121
122 T actual_pos = abs (actual[i]);
123 T nm = one_over_noise_squared * measured[i];
124 T nms = nm * actual_pos;
125 T F0 = Bessel::I0_scaled (nms);
126 T m_a = measured[i] - actual_pos;
127 T nm_a = one_over_noise_squared * m_a;
128 T F1_F0 = (Bessel::I1_scaled (nms) - F0) / F0;
129 dP_dactual[i] = -nm_a - nm * F1_F0;
130 if (actual[i] < 0.0) dP_dactual[i] = -dP_dactual[i];
131 lnP += 0.5 * nm_a * m_a - log (nm * F0);
132 assert (std::isfinite (lnP));
133 }
134
135 return lnP;
136 }
137
138
139 template <typename T> inline T lnP (const Vector<T>& measured, const Vector<T>& actual, const T one_over_noise_squared, Vector<T>& dP_dactual, T& dP_dN)
140 {
141 assert (one_over_noise_squared > 0.0);
142 assert (measured.size() == actual.size());
143 assert (measured.size() == dP_dactual.size());
144
145 T lnP = 0.0;
146 dP_dN = -T (measured.size()) / one_over_noise_squared;
147
148 for (size_t i = 0; i < measured.size(); i++) {
149 assert (measured[i] > 0.0);
150
151 T actual_pos = abs (actual[i]);
152 T nm = one_over_noise_squared * measured[i];
153 T nms = nm * actual_pos;
154 T F0 = Bessel::I0_scaled (nms);
155 T m_a = measured[i] - actual_pos;
156 T nm_a = one_over_noise_squared * m_a;
157 T F1_F0 = (Bessel::I1_scaled (nms) - F0) / F0;
158 dP_dactual[i] = -nm_a - nm * F1_F0;
159 if (actual[i] < 0.0) dP_dactual[i] = -dP_dactual[i];
160 dP_dN += 0.5 * pow2 (m_a) - measured[i] * actual_pos * F1_F0;
161 lnP += 0.5 * nm_a * m_a - std::log (nm * F0);
162 assert (std::isfinite (dP_dN));
163 assert (std::isfinite (lnP));
164 }
165
166 return lnP;
167 }
168
169 }
170 }
171}
172
173#endif
MR::Math::pow2
constexpr T pow2(const T &v)
Definition: math.h:53
MR::Math::Bessel::I0_scaled
T I0_scaled(const T x)
Definition: bessel.h:42
MR::Math::Bessel::I1_scaled
T I1_scaled(const T x)
Definition: bessel.h:53
MR::Math::Rician::lnP
T lnP(const T measured, const T actual, const T one_over_noise_squared)
Definition: rician.h:31
MR
Definition: base.h:24
MR::abs
constexpr std::enable_if< std::is_arithmetic< X >::value &&std::is_unsigned< X >::value, X >::type abs(X x)
Definition: types.h:297