2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/blob_functions.cpp $
5 $LastChangedRevision: 45 $
6 $Date: 2014-06-12 19:31:23 -0400 (Thu, 12 Jun 2014) $
8 ***********************************************************
10 There are three subroutines:
12 1. dc computes the DC term of the blob defined by the function kaiser
13 (i.e., Robert Lewitt's original definition of a blob) blob function at
16 2. radial computes the normalized blob function at distance "dist", and
18 3. lintegral computes the line integral of a blob at distance "dist".
29 REAL kaiser_proj(REAL s);
30 INTEGER kaiser(REAL* temp, REAL dist);
31 static REAL bessi3(REAL shape);
32 static REAL bessi2(REAL shape);
37 REAL r = Blob.support;
38 REAL shape = Blob.shape;
39 //REAL fshape = fabs(shape); shape is supposted to be non-negative
41 ans = 2 * Consts.pi * r * r * bessi3(shape) / (shape * bessi2(shape));
43 //fprintf(output,"\ndc = %f ", ans);
49 REAL Blob_class::radial(REAL dist)
55 // now comes the normalization
56 ans = ans * Consts.sqrt3 / 2.0 * Blob.delta * Blob.delta / dc;
58 //fprintf(output,"\nblob factor = %f ", sqrt3 / 2.0 * Blob.delta * Blob.delta / dc);
63 REAL Blob_class::lintegral(REAL dist)
68 ans = kaiser_proj(dist);
69 ans = ans * Consts.sqrt3 / 2.0 * Blob.delta * Blob.delta / dc; //normalisation
75 // Function to compute the Bessel 0
77 static REAL bessi0(REAL x)
83 y = (x * x) / (3.75 * 3.75);
84 ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
85 + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
89 ans = (exp(ax) / sqrt(ax)) * (0.39894228 + y * (0.1328592e-1
90 + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
91 + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
92 + y * 0.392377e-2))))))));
97 // Function to compute the Bessel 1
98 static REAL bessi1(REAL x)
106 ans = ax * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934
107 + y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3))))));
111 ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1
113 ans = 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2
114 + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
115 ans *= (exp(ax) / sqrt(ax));
118 return (x < 0.0)? -ans : ans;
122 // Function to compute the Bessel 2
124 static REAL bessi2(REAL x)
128 if(x <= Consts.zero) {
132 ans = bessi0(x) - (2.0 / x) * bessi1(x);
138 // Function to compute the Bessel 3
139 static REAL bessi3(REAL x)
143 if(x <= Consts.zero) {
147 ans = bessi1(x) - (4.0 / x) * bessi2(x);
151 /* Function to compute the value of the Spatial Kaiser-Bessel Function */
153 * dist stands for the radial distance from the blob center (variable)
155 * Blob.support stands for the radius
156 * Blob.shape for the decay (shape of the blob)
160 INTEGER kaiser(REAL* ans, REAL dist)
164 rda = dist / Blob.support;
167 arg = Blob.shape * sqrt(1.0 - rdas);
169 if(fabs(Blob.shape) < Consts.zero) {
170 // if(Blob.shape==0.0)
171 w = (1.0 - rdas) * (1.0 - rdas);
174 w = (1.0 - rdas) * bessi2(arg) / bessi2(Blob.shape);
175 //fprintf(output,"\ny=%f",*y);
188 // Function to compute the value of the Spatial Kaiser-Bessel Function
190 // Value of line integral through Kaiser-Bessel radial function
191 // (n >=2 dimensions) at distance s from center of function.
192 // Blob.support stands for the radius
193 // Blob.shape for the decay (shape of the blob)
195 REAL kaiser_proj(REAL s)
197 REAL sda, sdas, w, arg, p;
199 sda = s / Blob.support;
203 arg = Blob.shape * sqrt(w);
205 if(Blob.shape == 0.0) {
206 p = 2.0 * Blob.support * w * w * sqrt(w) * (8.0 / 15.0);
209 p = (2.0 * Blob.support/Blob.shape) * w *\
210 ((3.0/(arg * arg) + 1.0) * sinh(arg) - (3.0/arg) * cosh(arg) ) / bessi2(Blob.shape);