From 1db6071087434db4b73479a2eac9794c29d58b31 Mon Sep 17 00:00:00 2001
From: Claudio Gheller <cgheller@login01.m100.cineca.it>
Date: Mon, 20 Mar 2023 17:52:55 +0100
Subject: [PATCH] Kaiser-Bessel and tabulated Gaussian available on CPU. Bug
fixed in uvw normalization
---
Makefile | 2 +-
w-stacking-fftw.c | 8 ++++----
w-stacking.cu | 38 +++++++++++++++++++++++++++++---------
3 files changed, 34 insertions(+), 14 deletions(-)
diff --git a/Makefile b/Makefile
index ce52b64..dc58e4e 100644
--- a/Makefile
+++ b/Makefile
@@ -40,7 +40,7 @@ OPT += -DWRITE_IMAGE
#OPT += -DPARALLELIO
# Normalize uvw in case it is not done in the binMS
OPT += -DNORMALIZE_UVW
-# Gridding kernel: GAUSS, GAUSS_HI_PRECISION
+# Gridding kernel: GAUSS, GAUSS_HI_PRECISION, KAISERBESSEL
OPT += -DGAUSS_HI_PRECISION
ifeq (FITSIO,$(findstring FITSIO,$(OPT)))
diff --git a/w-stacking-fftw.c b/w-stacking-fftw.c
index ef56aa4..93a0c97 100644
--- a/w-stacking-fftw.c
+++ b/w-stacking-fftw.c
@@ -354,12 +354,12 @@ if(rank == 0){
maxv_all = maxv;
maxw_all = maxw;
#endif
- double offset = 0.0;
- double ming = MIN(minu_all,minv_all);
- double maxg = MAX(maxu_all,maxv_all);
+ double offset = 0.001;
+ double ming = MAX(abs(minu_all),abs(minv_all));
+ double maxg = MAX(abs(maxu_all),abs(maxv_all));
+ maxg = MAX(maxg,ming);
minw = minw_all;
maxw = maxw_all;
- ming = ming-offset*ming;
maxg = maxg+offset*maxg;
for (long inorm=0; inorm<Nmeasures; inorm++)
{
diff --git a/w-stacking.cu b/w-stacking.cu
index 049dfb5..cbad05c 100644
--- a/w-stacking.cu
+++ b/w-stacking.cu
@@ -31,11 +31,12 @@ void makeGaussKernel(double * kernel,
double norm = std22/PI;
int n = increaseprecision*KernelLen, mid = n / 2;
for (int i = 0; i != mid + 1; i++) {
- double term = (double)i / mid;
+ double term = (double)i/(double)increaseprecision;
kernel[mid + i] = sqrt(norm) * exp(-(term*term)*std22);
}
for (int i = 0; i != mid; i++) kernel[i] = kernel[n - 1 - i];
+// for (int i = 0; i < n; i++) printf("%f\n",kernel[i]);
}
@@ -53,10 +54,11 @@ double bessel0(double x, double precision) {
}
void makeKaiserBesselKernel(double * kernel,
int KernelLen,
+ int increaseprecision,
double alpha,
double overSamplingFactor,
int withSinc) {
- int n = KernelLen, mid = n / 2;
+ int n = increaseprecision*KernelLen, mid = n / 2;
double * sincKernel = malloc((mid + 1)*sizeof(*sincKernel));
const double filterRatio = 1.0 / overSamplingFactor;
sincKernel[0] = filterRatio;
@@ -69,10 +71,11 @@ void makeKaiserBesselKernel(double * kernel,
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;
+ bessel0(alpha * sqrt(1.0 - (term * term)), 1e-8) *
+ normFactor;
}
for (int i = 0; i != mid; i++) kernel[i] = kernel[n - 1 - i];
+ //for (int i = 0; i < n; i++) printf("%f\n",kernel[i]);
}
@@ -202,12 +205,21 @@ void wstack(
// initialize the convolution kernel
// gaussian:
int KernelLen = (w_support-1)/2;
- int increaseprecision = 21; // this number must be odd: increaseprecison*w_support must be odd (w_support must be odd)
+ int increaseprecision = 5; // this number must be odd: increaseprecison*w_support must be odd (w_support must be odd)
double std = 1.0;
double std22 = 1.0/(2.0*std*std);
double norm = std22/PI;
double * convkernel = malloc(increaseprecision*w_support*sizeof(*convkernel));
+ double overSamplingFactor = 1.0;
+ int withSinc = 0;
+ double alpha = 8.6;
+ #ifdef GAUSS
makeGaussKernel(convkernel,w_support,increaseprecision,std22);
+ #endif
+ #ifdef KAISERBESSEL
+ makeKaiserBesselKernel(convkernel, w_support, increaseprecision, alpha, overSamplingFactor, withSinc);
+ #endif
+
// Loop over visibilities.
// Switch between CUDA and GPU versions
@@ -325,21 +337,29 @@ void wstack(
for (k = kmin; k <= kmax; k++)
{
- double v_dist = (double)k+0.5 - pos_v;
+ //double v_dist = (double)k+0.5 - pos_v;
+ double v_dist = (double)k - pos_v;
for (j = jmin; j <= jmax; j++)
{
- double u_dist = (double)j+0.5 - pos_u;
+ //double u_dist = (double)j+0.5 - pos_u;
+ double u_dist = (double)j - pos_u;
long iKer = 2 * (j + k*grid_size_x + grid_w*grid_size_x*grid_size_y);
- int jKer = (int)(increaseprecision*u_dist) + KernelLen;
- int kKer = (int)(increaseprecision*v_dist) + KernelLen;
+ int jKer = (int)(increaseprecision * (fabs(u_dist+(double)KernelLen)));
+ int kKer = (int)(increaseprecision * (fabs(v_dist+(double)KernelLen)));
#ifdef GAUSS_HI_PRECISION
double conv_weight = gauss_kernel_norm(norm,std22,u_dist,v_dist);
#endif
#ifdef GAUSS
double conv_weight = convkernel[jKer]*convkernel[kKer];
+ //if(jKer < 0 || jKer >= 35 || kKer < 0 || kKer >= 35)
+ // printf("%f %d %f %d\n",fabs(u_dist+(double)KernelLen),jKer,fabs(v_dist+(double)KernelLen),kKer);
+ //printf("%d %d %d %d %f %f %f %f %f\n",jKer, j, kKer, k, pos_u, pos_v, u_dist,v_dist,conv_weight);
#endif
+ #ifdef KAISERBESSEL
+ double conv_weight = convkernel[jKer]*convkernel[kKer];
+ #endif
// Loops over frequencies and polarizations
double add_term_real = 0.0;
double add_term_img = 0.0;
--
GitLab