src/snark/blob_functions.cpp

   1 /*
   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) $
   7 $Author: olangthaler $
   8 ***********************************************************
   9
  10  There are three subroutines:
  11
  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
  14     distance "dist",
  15
  16  2. radial computes the normalized blob function at distance "dist", and
  17
  18  3. lintegral computes the line integral of a blob at distance "dist".
  19
  20 */
  21
  22 #include <cstdio>
  23 #include <cmath>
  24
  25 #include "consts.h"
  26 #include "blkdta.h"
  27 #include "blob.h"
  28
  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);
  33
  34 REAL Blob_class::dc()
  35 {
  36   REAL ans;
  37   REAL r = Blob.support;
  38   REAL shape = Blob.shape;
  39   //REAL fshape = fabs(shape); shape is supposted to be non-negative
  40
  41   ans = 2 * Consts.pi * r * r * bessi3(shape) / (shape * bessi2(shape));
  42   //tracer
  43   //fprintf(output,"\ndc = %f  ", ans);
  44   //tracer
  45   return ans;
  46 }
  47
  48
  49 REAL Blob_class::radial(REAL dist)
  50 {
  51   REAL ans;
  52   REAL dc = Blob.dc();
  53
  54   kaiser(&ans, dist);
  55   // now comes the normalization
  56   ans = ans * Consts.sqrt3 / 2.0 * Blob.delta * Blob.delta / dc;
  57  //tracer
  58   //fprintf(output,"\nblob factor = %f  ", sqrt3 / 2.0 * Blob.delta * Blob.delta / dc);
  59   //tracer
  60    return ans;
  61 }
  62
  63 REAL Blob_class::lintegral(REAL dist)
  64 {
  65   REAL ans;
  66   REAL dc = Blob.dc();
  67
  68   ans = kaiser_proj(dist);
  69   ans = ans * Consts.sqrt3 / 2.0 * Blob.delta * Blob.delta / dc;  //normalisation
  70
  71   return ans;
  72 }
  73
  74
  75 // Function to compute the Bessel 0
  76
  77 static REAL bessi0(REAL x)
  78 {
  79   REAL y,ax,ans;
  80
  81   ax = fabs(x);
  82   if(ax < 3.75) {
  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)))));
  86   }
  87   else {
  88     y = 3.75 / ax;
  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))))))));
  93   }
  94   return ans;
  95 }
  96
  97 // Function to compute the Bessel 1
  98 static REAL bessi1(REAL x)
  99 {
 100   REAL ax, ans, y;
 101
 102   ax = fabs(x);
 103   if(ax < 3.75) {
 104     y = x / 3.75;
 105     y *= y;
 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))))));
 108   }
 109   else {
 110     y = 3.75 / ax;
 111     ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1
 112         - y * 0.420059e-2));
 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));
 116   }
 117
 118   return (x < 0.0)? -ans : ans;
 119 }
 120
 121
 122 // Function to compute the Bessel 2
 123
 124 static REAL bessi2(REAL x)
 125 {
 126   REAL ans;
 127
 128   if(x <= Consts.zero) {
 129     ans = 0.0;
 130   }
 131   else {
 132     ans = bessi0(x) - (2.0 / x) * bessi1(x);
 133   }
 134
 135   return ans;
 136 }
 137
 138 // Function to compute the Bessel 3
 139 static REAL bessi3(REAL x)
 140 {
 141   REAL ans;
 142
 143   if(x <= Consts.zero) {
 144     ans = 0.0;
 145   }
 146   else {
 147     ans = bessi1(x) - (4.0 / x) * bessi2(x);
 148   }
 149   return ans;
 150 }
 151 /* Function to compute the value of the Spatial Kaiser-Bessel Function */
 152 /*
 153  * dist stands for the radial distance from the blob center (variable)
 154  * CONSTANTS:
 155  * Blob.support stands for the radius
 156  * Blob.shape for the decay (shape of the blob)
 157  *
 158  */
 159
 160 INTEGER kaiser(REAL* ans, REAL dist)
 161 {
 162   REAL rda,rdas,arg,w;
 163
 164   rda = dist / Blob.support;
 165   rdas = rda * rda;
 166   if(rdas <= 1.0) {
 167     arg = Blob.shape * sqrt(1.0 - rdas);
 168
 169     if(fabs(Blob.shape) < Consts.zero) {
 170       // if(Blob.shape==0.0)
 171       w = (1.0 - rdas) * (1.0 - rdas);
 172     }
 173     else {
 174       w = (1.0 - rdas) * bessi2(arg) / bessi2(Blob.shape);
 175       //fprintf(output,"\ny=%f",*y);
 176     }
 177   }
 178   else {
 179     w = 0.0;
 180   }
 181
 182   *ans = w;
 183
 184   return 1;
 185 }
 186
 187
 188 // Function to compute the value of the Spatial Kaiser-Bessel Function
 189 //
 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)
 194
 195 REAL kaiser_proj(REAL s)
 196 {
 197   REAL sda, sdas, w, arg, p;
 198
 199   sda = s / Blob.support;
 200   sdas = sda * sda;
 201   w = 1.0 - sdas;
 202   if(w > 1.0e-10) {
 203     arg = Blob.shape * sqrt(w);
 204
 205     if(Blob.shape == 0.0) {
 206       p = 2.0 * Blob.support * w * w * sqrt(w) * (8.0 / 15.0);
 207     }
 208     else {
 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);
 211     }
 212   }
 213   else {
 214     p = 0.0;
 215   }
 216
 217   return p;
 218 }
 219
 220