From 02194d4beb5f4c0492210a9bf44a13ca892dce9e Mon Sep 17 00:00:00 2001
From: Claudio Gheller <cgheller@login01.m100.cineca.it>
Date: Sat, 18 Mar 2023 23:15:08 +0100
Subject: [PATCH] Tabulated gaussian kernel with increased precision added
(only CPU version)
---
Makefile | 2 ++
w-stacking.cu | 19 ++++++++++++-------
2 files changed, 14 insertions(+), 7 deletions(-)
diff --git a/Makefile b/Makefile
index 9bdba20..ce52b64 100644
--- a/Makefile
+++ b/Makefile
@@ -40,6 +40,8 @@ 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
+OPT += -DGAUSS_HI_PRECISION
ifeq (FITSIO,$(findstring FITSIO,$(OPT)))
LIBS += -L$(FITSIO_LIB) -lcfitsio
diff --git a/w-stacking.cu b/w-stacking.cu
index 8728fda..049dfb5 100644
--- a/w-stacking.cu
+++ b/w-stacking.cu
@@ -24,11 +24,12 @@ gauss_kernel_norm(double norm, double std22, double u_dist, double v_dist)
void makeGaussKernel(double * kernel,
int KernelLen,
+ int increaseprecision,
double std22)
{
double norm = std22/PI;
- int n = KernelLen, mid = n / 2;
+ int n = increaseprecision*KernelLen, mid = n / 2;
for (int i = 0; i != mid + 1; i++) {
double term = (double)i / mid;
kernel[mid + i] = sqrt(norm) * exp(-(term*term)*std22);
@@ -201,11 +202,12 @@ 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)
double std = 1.0;
double std22 = 1.0/(2.0*std*std);
double norm = std22/PI;
- double * convkernel = malloc(w_support*sizeof(*convkernel));
- makeGaussKernel(convkernel,w_support,std22);
+ double * convkernel = malloc(increaseprecision*w_support*sizeof(*convkernel));
+ makeGaussKernel(convkernel,w_support,increaseprecision,std22);
// Loop over visibilities.
// Switch between CUDA and GPU versions
@@ -327,14 +329,17 @@ void wstack(
for (j = jmin; j <= jmax; j++)
{
- int jKer = (int)j-pos_u+KernelLen;
- int kKer = (int)k-pos_v+KernelLen;
double u_dist = (double)j+0.5 - 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;
-
- //double conv_weight = gauss_kernel_norm(norm,std22,u_dist,v_dist);
+ #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];
+ #endif
// Loops over frequencies and polarizations
double add_term_real = 0.0;
double add_term_img = 0.0;
--
GitLab