Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkFresnel.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 /****************************************************************************
14  * fresnel.c -
15  * Calculation of Fresnel integrals by expansion to Chebyshev series
16  * Expansions are taken from the book
17  * Y.L. Luke. Mathematical functions and their approximations.
18  * Moscow, "Mir", 1980. PP. 145-149 (Russian edition)
19  ****************************************************************************
20  */
21 
22 #include <math.h>
23 
24 static const double sqrt_2_pi = 0.7978845608028653558798921199; /* sqrt(2/pi) */
25 static const double _1_sqrt_2pi = 0.3989422804014326779399460599; /* 1/sqrt(2*pi) */
26 static const double pi_2 = 1.5707963267948966192313216916; /* pi/2 */
27 
28 static double f_data_a[18] =
29 {
30  0.76435138664186000189,
31  -0.43135547547660179313,
32  0.43288199979726653054,
33  -0.26973310338387111029,
34  0.08416045320876935378,
35  -0.01546524484461381958,
36  0.00187855423439822018,
37  -0.00016264977618887547,
38  0.00001057397656383260,
39  -0.00000053609339889243,
40  0.00000002181658454933,
41  -0.00000000072901621186,
42  0.00000000002037332546,
43  -0.00000000000048344033,
44  0.00000000000000986533,
45  -0.00000000000000017502,
46  0.00000000000000000272,
47  -0.00000000000000000004
48 };
49 
50 static double f_data_b[17] =
51 {
52  0.63041404314570539241,
53  -0.42344511405705333544,
54  0.37617172643343656625,
55  -0.16249489154509567415,
56  0.03822255778633008694,
57  -0.00564563477132190899,
58  0.00057454951976897367,
59  -0.00004287071532102004,
60  0.00000245120749923299,
61  -0.00000011098841840868,
62  0.00000000408249731696,
63  -0.00000000012449830219,
64  0.00000000000320048425,
65  -0.00000000000007032416,
66  0.00000000000000133638,
67  -0.00000000000000002219,
68  0.00000000000000000032
69 };
70 
71 static double fresnel_cos_0_8(double x)
72 {
73  double x_8 = x / 8.0;
74  double xx = 2.0*x_8*x_8 - 1.0;
75 
76  double t0 = 1.0;
77  double t1 = xx;
78  double sumC = f_data_a[0] + f_data_a[1] * t1;
79  double t2;
80  int n;
81  for (n = 2; n < 18; n++)
82  {
83  t2 = 2.0*xx*t1 - t0;
84  sumC += f_data_a[n] * t2;
85  t0 = t1; t1 = t2;
86  }
87  return _1_sqrt_2pi * sqrt(x)*sumC;
88 }
89 
90 static double fresnel_sin_0_8(double x)
91 {
92  double x_8 = x / 8.0;
93  double xx = 2.0*x_8*x_8 - 1.0;
94  double t0 = 1.;
95  double t1 = xx;
96  double ot1 = x_8;
97  double ot2 = 2.0*x_8*t1 - ot1;
98  double sumS = f_data_b[0] * ot1 + f_data_b[1] * ot2;
99  int n;
100  double t2;
101  for (n = 2; n < 17; n++)
102  {
103  t2 = 2.0*xx*t1 - t0;
104  ot1 = ot2;
105  ot2 = 2.0*x_8*t2 - ot1;
106  sumS += f_data_b[n] * ot2;
107  t0 = t1; t1 = t2;
108  }
109  return _1_sqrt_2pi * sqrt(x)*sumS;
110 }
111 
112 static double f_data_e[41] =
113 {
114  0.97462779093296822410,
115  -0.02424701873969321371,
116  0.00103400906842977317,
117  -0.00008052450246908016,
118  0.00000905962481966582,
119  -0.00000131016996757743,
120  0.00000022770820391497,
121  -0.00000004558623552026,
122  0.00000001021567537083,
123  -0.00000000251114508133,
124  0.00000000066704761275,
125  -0.00000000018931512852,
126  0.00000000005689898935,
127  -0.00000000001798219359,
128  0.00000000000594162963,
129  -0.00000000000204285065,
130  0.00000000000072797580,
131  -0.00000000000026797428,
132  0.00000000000010160694,
133  -0.00000000000003958559,
134  0.00000000000001581262,
135  -0.00000000000000646411,
136  0.00000000000000269981,
137  -0.00000000000000115038,
138  0.00000000000000049942,
139  -0.00000000000000022064,
140  0.00000000000000009910,
141  -0.00000000000000004520,
142  0.00000000000000002092,
143  -0.00000000000000000982,
144  0.00000000000000000467,
145  -0.00000000000000000225,
146  0.00000000000000000110,
147  -0.00000000000000000054,
148  0.00000000000000000027,
149  -0.00000000000000000014,
150  0.00000000000000000007,
151  -0.00000000000000000004,
152  0.00000000000000000002,
153  -0.00000000000000000001,
154  0.00000000000000000001
155 };
156 
157 static double f_data_f[35] =
158 {
159  0.99461545179407928910,
160  -0.00524276766084297210,
161  0.00013325864229883909,
162  -0.00000770856452642713,
163  0.00000070848077032045,
164  -0.00000008812517411602,
165  0.00000001359784717148,
166  -0.00000000246858295747,
167  0.00000000050925789921,
168  -0.00000000011653400634,
169  0.00000000002906578309,
170  -0.00000000000779847361,
171  0.00000000000222802542,
172  -0.00000000000067239338,
173  0.00000000000021296411,
174  -0.00000000000007041482,
175  0.00000000000002419805,
176  -0.00000000000000861080,
177  0.00000000000000316287,
178  -0.00000000000000119596,
179  0.00000000000000046444,
180  -0.00000000000000018485,
181  0.00000000000000007527,
182  -0.00000000000000003131,
183  0.00000000000000001328,
184  -0.00000000000000000574,
185  0.00000000000000000252,
186  -0.00000000000000000113,
187  0.00000000000000000051,
188  -0.00000000000000000024,
189  0.00000000000000000011,
190  -0.00000000000000000005,
191  0.00000000000000000002,
192  -0.00000000000000000001,
193  0.00000000000000000001
194 };
195 
196 static double fresnel_cos_8_inf(double x)
197 {
198  double xx = 128.0 / (x*x) - 1.0; /* 2.0*(8/x)^2 - 1 */
199  double t0 = 1.0;
200  double t1 = xx;
201  double sumP = f_data_e[0] + f_data_e[1] * t1;
202  double sumQ = f_data_f[0] + f_data_f[1] * t1;
203  double t2;
204  int n;
205  for (n = 2; n < 35; n++)
206  {
207  t2 = 2.0*xx*t1 - t0;
208  sumP += f_data_e[n] * t2; /* sumP += f_data_e[n]*ChebyshevT(n,xx) */
209  sumQ += f_data_f[n] * t2; /* sumQ += f_data_f[n]*ChebyshevT(n,xx) */
210  t0 = t1; t1 = t2;
211  }
212  for (n = 35; n < 41; n++)
213  {
214  t2 = 2.0*xx*t1 - t0;
215  sumP += f_data_e[n] * t2; /* sumP += f_data_e[n]*ChebyshevT(n,xx) */
216  t0 = t1; t1 = t2;
217  }
218  return 0.5 - _1_sqrt_2pi * (0.5*sumP*cos(x) / x - sumQ * sin(x)) / sqrt(x);
219 }
220 
221 static double fresnel_sin_8_inf(double x)
222 {
223  double xx = 128.0 / (x*x) - 1.0; /* 2.0*(8/x)^2 - 1 */
224  double t0 = 1.0;
225  double t1 = xx;
226  double sumP = f_data_e[0] + f_data_e[1] * t1;
227  double sumQ = f_data_f[0] + f_data_f[1] * t1;
228  double t2;
229  int n;
230  for (n = 2; n < 35; n++)
231  {
232  t2 = 2.0*xx*t1 - t0;
233  sumP += f_data_e[n] * t2; /* sumP += f_data_e[n]*ChebyshevT(n,xx) */
234  sumQ += f_data_f[n] * t2; /* sumQ += f_data_f[n]*ChebyshevT(n,xx) */
235  t0 = t1; t1 = t2;
236  }
237  for (n = 35; n < 41; n++)
238  {
239  t2 = 2.0*xx*t1 - t0;
240  sumP += f_data_e[n] * t2; /* sumQ += f_data_f[n]*ChebyshevT(n,xx) */
241  t0 = t1; t1 = t2;
242  }
243  return 0.5 - _1_sqrt_2pi * (0.5*sumP*sin(x) / x + sumQ * cos(x)) / sqrt(x);
244 }
245 
246 
247 namespace mitk
248 {
249 
250  double fresnel_c(double x)
251  {
252  double xx = x * x*pi_2;
253  double ret_val;
254  if (xx <= 8.0)
255  ret_val = fresnel_cos_0_8(xx);
256  else
257  ret_val = fresnel_cos_8_inf(xx);
258  return (x < 0.0) ? -ret_val : ret_val;
259  }
260 
261  double fresnel_s(double x)
262  {
263  double xx = x * x*pi_2;
264  double ret_val;
265  if (xx <= 8.0)
266  ret_val = fresnel_sin_0_8(xx);
267  else
268  ret_val = fresnel_sin_8_inf(xx);
269  return (x < 0.0) ? -ret_val : ret_val;
270  }
271 
272  double fresnel_c2(double x)
273  {
274  return fresnel_c(x*sqrt_2_pi);
275  }
276 
277  double fresnel_s2(double x)
278  {
279  return fresnel_s(x*sqrt_2_pi);
280  }
281 
282 }
static double fresnel_sin_8_inf(double x)
DataCollection - Class to facilitate loading/accessing structured data.
double fresnel_c(double x)
static double f_data_a[18]
Definition: mitkFresnel.cpp:28
static double f_data_f[35]
static double f_data_e[41]
static double fresnel_sin_0_8(double x)
Definition: mitkFresnel.cpp:90
bool t2(false)
double fresnel_s(double x)
double fresnel_c2(double x)
static double fresnel_cos_0_8(double x)
Definition: mitkFresnel.cpp:71
static const double sqrt_2_pi
Definition: mitkFresnel.cpp:24
static const double _1_sqrt_2pi
Definition: mitkFresnel.cpp:25
static const double pi_2
Definition: mitkFresnel.cpp:26
static double f_data_b[17]
Definition: mitkFresnel.cpp:50
double fresnel_s2(double x)
static double fresnel_cos_8_inf(double x)