Developer documentation
Version 3.0.3-105-gd3941f44
gaussian.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_gaussian_h__
18#define __math_gaussian_h__
19
20#include "math/math.h"
21#include "math/vector.h"
22
23namespace MR {
24 namespace Math {
25 namespace Gaussian {
26
27 template <typename T> inline T lnP (const T measured, const T actual, const T one_over_noise_squared)
28 {
29 return (0.5 * (one_over_noise_squared * pow2 (actual - measured) - log (one_over_noise_squared)));
30 }
31
32
33
34
35 template <typename T> inline T lnP (const T measured, const T actual, const T one_over_noise_squared, T& dP_dactual, T& dP_dN)
36 {
37 assert (one_over_noise_squared > 0.0);
38 dP_dN = pow2 (actual - measured);
39 dP_dactual = one_over_noise_squared * dP_dN;
40 dP_dN = 0.5 * (dP_dN - 1.0/one_over_noise_squared);
41 return (0.5 * (dP_dactual - log(one_over_noise_squared)));
42 }
43
44
45
46
47
48 template <typename T> inline T lnP (const int N, const T* measured, const T* actual, const T one_over_noise_squared)
49 {
50 assert (one_over_noise_squared > 0.0);
51
52 T lnP = 0.0;
53 for (int i = 0; i < N; i++)
54 lnP += pow2 (actual[i] - measured[i]);
55 lnP *= one_over_noise_squared;
56 lnP -= N * log (one_over_noise_squared);
57 return (0.5*lnP);
58 }
59
60
61
62 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)
63 {
64 assert (one_over_noise_squared > 0.0);
65
66 T lnP = 0.0;
67 dP_dN = 0.0;
68 for (int i = 0; i < N; i++) {
69 T diff = actual[i] - measured[i];
70 dP_dactual[i] = one_over_noise_squared * diff;
71 lnP += pow2 (diff);
72 }
73 dP_dN = 0.5 * (lnP - T(N)/one_over_noise_squared);
74 lnP *= one_over_noise_squared;
75 lnP -= N * log (one_over_noise_squared);
76
77 return (0.5*lnP);
78 }
79
80
81
82 template <typename T> inline T lnP (const Vector<T>& measured, const Vector<T>& actual, const T one_over_noise_squared)
83 {
84 assert (one_over_noise_squared > 0.0);
85 assert (measured.size() == actual.size());
86
87 T lnP = 0.0;
88 for (size_t i = 0; i < measured.size(); i++)
89 lnP += pow2 (actual[i] - measured[i]);
90 lnP *= one_over_noise_squared;
91 lnP -= measured.size() * log (one_over_noise_squared);
92 return (0.5*lnP);
93 }
94
95
96
97
98 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)
99 {
100 assert (one_over_noise_squared > 0.0);
101 assert (measured.size() == actual.size());
102 assert (measured.size() == dP_dactual.size());
103
104 T lnP = 0.0;
105 dP_dN = 0.0;
106 for (size_t i = 0; i < measured.size(); i++) {
107 T diff = actual[i] - measured[i];
108 dP_dactual[i] = one_over_noise_squared * diff;
109 lnP += pow2 (diff);
110 }
111 dP_dN = 0.5 * (lnP - T(measured.size())/one_over_noise_squared);
112 lnP *= one_over_noise_squared;
113 lnP -= measured.size() * log (one_over_noise_squared);
114
115 return (0.5*lnP);
116 }
117
118 }
119 }
120}
121
122#endif
constexpr T pow2(const T &v)
Definition: math.h:53
T lnP(const T measured, const T actual, const T one_over_noise_squared)
Definition: gaussian.h:27
Definition: base.h:24