Kaiser-Bessel kernel added (adapted from WSClean). Not yet used

w-stacking.cu

@@ -14,12 +14,51 @@ double __device__
 #else
 double
 #endif
+// Gaussian Kernel
 gauss_kernel_norm(double norm, double std22, double u_dist, double v_dist)
 {
      double conv_weight;
      conv_weight = norm * exp(-((u_dist*u_dist)+(v_dist*v_dist))*std22);
      return conv_weight;
 }
+
+// Kaiser-Bessel Kernel: it is adapted from WSClean
+double bessel0(double x, double precision) {
+  // Calculate I_0 = SUM of m 0 -> inf [ (x/2)^(2m) ]
+  // This is the unnormalized bessel function of order 0.
+  double d = 0.0, ds = 1.0, sum = 1.0;
+  do {
+    d += 2.0;
+    ds *= x * x / (d * d);
+    sum += ds;
+  } while (ds > sum * precision);
+  return sum;
+}
+void makeKaiserBesselKernel(double * kernel,
+		            int KernelLen,
+                            double alpha,
+                            double overSamplingFactor,
+                            int withSinc) {
+  int n = KernelLen, mid = n / 2;
+  double * sincKernel = malloc((mid + 1)*sizeof(*sincKernel));
+  const double filterRatio = 1.0 / overSamplingFactor;
+  sincKernel[0] = filterRatio;
+  for (int i = 1; i != mid + 1; i++) {
+    double x = i;
+    sincKernel[i] =
+        withSinc ? (sin(PI * filterRatio * x) / (PI * x)) : filterRatio;
+  }
+  const double normFactor = overSamplingFactor / bessel0(alpha, 1e-8);
+  for (int i = 0; i != mid + 1; i++) {
+    double term = (double)i / mid;
+    kernel[mid + i] = sincKernel[i] *
+                      bessel0(alpha * sqrt(1.0 - (term * term)), 1e-8) *
+                      normFactor;
+  }
+  for (int i = 0; i != mid; i++) kernel[i] = kernel[n - 1 - i];
+}
+
+
 #ifdef ACCOMP
 #pragma omp end declare target
 #endif