From 2d16e4a010abd8e0b8ed26cdc45fa082e10b8cd5 Mon Sep 17 00:00:00 2001
From: Claudio Gheller <cgheller@login01.m100.cineca.it>
Date: Thu, 16 Mar 2023 09:42:41 +0100
Subject: [PATCH] Kaiser-Bessel kernel added (adapted from WSClean). Not yet
used
---
w-stacking.cu | 39 +++++++++++++++++++++++++++++++++++++++
1 file changed, 39 insertions(+)
diff --git a/w-stacking.cu b/w-stacking.cu
index 4fe3625..d5e8e8c 100644
--- a/w-stacking.cu
+++ b/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
--
GitLab