1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2008-2015 Imperial College London
5 * Copyright 2008-2015 Daniel Rueckert, Julia Schnabel
6 * Copyright 2013-2015 Andreas Schuh
7 *
8 * Licensed under the Apache License, Version 2.0 (the "License");
9 * you may not use this file except in compliance with the License.
10 * You may obtain a copy of the License at
11 *
12 * http://www.apache.org/licenses/LICENSE-2.0
13 *
14 * Unless required by applicable law or agreed to in writing, software
15 * distributed under the License is distributed on an "AS IS" BASIS,
16 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 * See the License for the specific language governing permissions and
18 * limitations under the License.
19 */
20
21 #include "mirtk/ScalarGaussian.h"
22 #include "mirtk/Math.h"
23
24
25 namespace mirtk {
26
27
28 // =============================================================================
29 // Construction/Destruction
30 // =============================================================================
31
32 // -----------------------------------------------------------------------------
ComputeNorm()33 void ScalarGaussian::ComputeNorm()
34 {
35 _Norm[0] = (_Variance._x > .0 ? (1.0 / sqrt(two_pi * _Variance._x)) : 1.0);
36 _Norm[1] = (_Variance._y > .0 ? (_Norm[0] / sqrt(two_pi * _Variance._y)) : 1.0);
37 _Norm[2] = (_Variance._z > .0 ? (_Norm[1] / sqrt(two_pi * _Variance._z)) : 1.0);
38 _Norm[3] = (_Variance._t > .0 ? (_Norm[2] / sqrt(two_pi * _Variance._t)) : 1.0);
39 if (_Variance._x <= .0) _Variance._x = 1.0;
40 if (_Variance._y <= .0) _Variance._y = 1.0;
41 if (_Variance._z <= .0) _Variance._z = 1.0;
42 if (_Variance._t <= .0) _Variance._t = 1.0;
43 }
44
45 // -----------------------------------------------------------------------------
ScalarGaussian()46 ScalarGaussian::ScalarGaussian()
47 :
48 _Variance(1.0),
49 _Center(.0)
50 {
51 ComputeNorm();
52 }
53
54 // -----------------------------------------------------------------------------
ScalarGaussian(double sigma)55 ScalarGaussian::ScalarGaussian(double sigma)
56 :
57 _Variance(sigma * sigma),
58 _Center(.0)
59 {
60 ComputeNorm();
61 }
62
63 // -----------------------------------------------------------------------------
ScalarGaussian(double sigma,double x_0,double y_0,double z_0,double t_0)64 ScalarGaussian::ScalarGaussian(double sigma, double x_0, double y_0, double z_0, double t_0)
65 :
66 _Variance(sigma * sigma),
67 _Center(x_0, y_0, z_0, t_0)
68 {
69 ComputeNorm();
70 }
71
72 // -----------------------------------------------------------------------------
ScalarGaussian(double sigma_x,double sigma_y,double sigma_z,double x_0,double y_0,double z_0)73 ScalarGaussian::ScalarGaussian(double sigma_x, double sigma_y, double sigma_z,
74 double x_0, double y_0, double z_0)
75 :
76 _Variance(sigma_x*sigma_x, sigma_y*sigma_y, sigma_z*sigma_z, .0),
77 _Center(x_0, y_0, z_0, .0)
78 {
79 ComputeNorm();
80 }
81
82 // -----------------------------------------------------------------------------
ScalarGaussian(double sigma_x,double sigma_y,double sigma_z,double sigma_t,double x_0,double y_0,double z_0,double t_0)83 ScalarGaussian::ScalarGaussian(double sigma_x, double sigma_y, double sigma_z, double sigma_t,
84 double x_0, double y_0, double z_0, double t_0)
85 :
86 _Variance(sigma_x*sigma_x, sigma_y*sigma_y, sigma_z*sigma_z, sigma_t*sigma_t),
87 _Center(x_0, y_0, z_0, t_0)
88 {
89 ComputeNorm();
90 }
91
92 // -----------------------------------------------------------------------------
~ScalarGaussian()93 ScalarGaussian::~ScalarGaussian()
94 {
95 }
96
97 // =============================================================================
98 // Evaluation
99 // =============================================================================
100
101 // -----------------------------------------------------------------------------
Evaluate(double x)102 double ScalarGaussian::Evaluate(double x)
103 {
104 x -= _Center._x;
105 return _Norm[0] * exp(- 0.5 * (x * x) / _Variance._x);
106 }
107
108 // -----------------------------------------------------------------------------
Evaluate(double x,double y)109 double ScalarGaussian::Evaluate(double x, double y)
110 {
111 x -= _Center._x, y -= _Center._y;
112 return _Norm[1] * exp(- 0.5 * ( (x * x) / _Variance._x
113 + (y * y) / _Variance._y));
114 }
115
116 // -----------------------------------------------------------------------------
Evaluate(double x,double y,double z)117 double ScalarGaussian::Evaluate(double x, double y, double z)
118 {
119 x -= _Center._x, y -= _Center._y, z -= _Center._z;
120 return _Norm[2] * exp(- 0.5 * ( (x * x) / _Variance._x
121 + (y * y) / _Variance._y
122 + (z * z) / _Variance._z));
123 }
124
125 // -----------------------------------------------------------------------------
Evaluate(double x,double y,double z,double t)126 double ScalarGaussian::Evaluate(double x, double y, double z, double t)
127 {
128 x -= _Center._x, y -= _Center._y, z -= _Center._z, t -= _Center._t;
129 return _Norm[3] * exp(- 0.5 * ( (x * x) / _Variance._x
130 + (y * y) / _Variance._y
131 + (z * z) / _Variance._z
132 + (t * t) / _Variance._t));
133 }
134
135
136 } // namespace mirtk
137