From b854572eff26a66bc2da5a0e4f456dcc31e31e28 Mon Sep 17 00:00:00 2001 From: Claudio Gheller Date: Tue, 27 Dec 2022 17:29:48 +0100 Subject: [PATCH] Support to CFITSIO added to inverse-imaging --- Build/Makefile.M100 | 5 +- Makefile | 6 ++ inverse-imaging.c | 175 ++++++++++++++++++++++++++------------------ 3 files changed, 114 insertions(+), 72 deletions(-) diff --git a/Build/Makefile.M100 b/Build/Makefile.M100 index a1dde12..d675b68 100644 --- a/Build/Makefile.M100 +++ b/Build/Makefile.M100 @@ -8,6 +8,9 @@ MPIC++ = mpiCC FFTW_INCL= -I/cineca/prod/opt/libraries/fftw/3.3.8/spectrum_mpi--10.4.0--binary/include/ FFTW_LIB= /cineca/prod/opt/libraries/fftw/3.3.8/spectrum_mpi--10.4.0--binary/lib +FITSIO_INCL= -I/m100/home/userexternal/cgheller/cfitsio-3.49 +FITSIO_LIB= /m100/home/userexternal/cgheller/cfitsio-3.49 + NVCC = nvcc NVFLAGS = -arch=sm_70 -Xcompiler -mno-float128 -std=c++11 @@ -15,6 +18,6 @@ NVLIB = -L/cineca/prod/opt/compilers/cuda/10.1/none/lib64/ -lcudart -lcuda OMP= #-fopenmp -CFLAGS += -I. $(FFTW_INCL) $(GSL_INCL) $(MPI_INCL) +CFLAGS += -I. $(FFTW_INCL) $(GSL_INCL) $(MPI_INCL) $(FITSIO_INCL) OPTIMIZE = $(OMP) -O3 -mtune=native diff --git a/Makefile b/Makefile index b825783..93ab885 100644 --- a/Makefile +++ b/Makefile @@ -34,6 +34,12 @@ OPT += -DWRITE_DATA OPT += -DWRITE_IMAGE # perform w-stacking phase correction OPT += -DPHASE_ON +# Support CFITSIO +OPT += -DFITSIO + +ifeq (FITSIO,$(findstring FITSIO,$(OPT))) + LIBS += -L$(FITSIO_LIB) -lcfitsio +endif DEPS = w-stacking.h w-stacking-fftw.c w-stacking.cu phase_correction.cu diff --git a/inverse-imaging.c b/inverse-imaging.c index a3dfbbe..66c1603 100644 --- a/inverse-imaging.c +++ b/inverse-imaging.c @@ -7,6 +7,9 @@ #include #endif #endif +#ifdef FITSIO +#include "fitsio.h" +#endif #include #include #include @@ -80,7 +83,11 @@ int main(int argc, char * argv[]) // Image related files char imagepath[900] = "./"; + #ifdef FITSIO + char imagename[FILENAMELENGTH] = "image_name.fits"; + #else char imagename[FILENAMELENGTH] = "image_name.bin"; + #endif // Visibilities related variables double * uu; @@ -208,10 +215,7 @@ if(rank == 0){ // INPUT FILES (only the first ndatasets entries are used) int ndatasets = 1; - //strcpy(datapath_multi[0],"data/newgauss2noconj_t201806301100_SBL180.binMS/"); - //strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/newgauss4_t201806301100_SBL180.binMS/"); strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.binMS/"); - //strcpy(datapath_multi[1],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_134MHz.pre-cal.binMS/"); strcpy(datapath,datapath_multi[0]); // Read metadata @@ -329,18 +333,66 @@ if(rank == 0){ // define image variable double* image_real = (double*) calloc(xaxis*yaxis,sizeof(double)); + #ifdef WRITE_DATA + double* grid_real = (double*) calloc(xaxis*yaxis,sizeof(double)); + double* grid_img = (double*) calloc(xaxis*yaxis,sizeof(double)); + #endif // reading image strcpy(filename,imagepath); strcat(filename,imagename); printf("Reading Image %s\n",filename); + #ifdef FITSIO + fitsfile *fptr; + int status = 0; + int nkeys; + + fits_open_file(&fptr, filename, READONLY, &status); + + // header stuff + int nocomments; + char *header; + fits_hdr2str(fptr,0,NULL,0, &header,&nkeys,&status); + printf("Number of keywords in the Header = %d\n",nkeys); + int header_aux = (int)(80*nkeys/2880); + int header_size = (header_aux+1)*2880; + printf("Header size (bytes) = %d\n\n",header_size); + + // get image size + int naxis; + long * naxes; + fits_get_img_dim(fptr, &naxis, &status); + naxes = (long *)malloc(sizeof(long)*naxis); + fits_get_img_size(fptr, naxis, naxes, &status); + if (naxes[0] != xaxis){ + printf("Wrong Image Size : %d vs %d\n",naxes[0],xaxis); + exit(2); + } + + long * fpixel = (long *) malloc(sizeof(long)*naxis); + long * lpixel = (long *) malloc(sizeof(long)*naxis); + long * inc = (long *) malloc(sizeof(long)*naxis); + int anynul; + + fpixel[0] = 1; + lpixel[0] = xaxis; + inc[0] = 1; + fpixel[1] = 1; + lpixel[1] = rank*yaxis+1; + inc[1] = 1; + + fits_read_subset(fptr, TFLOAT, fpixel, lpixel, inc, NULL, image_real, &anynul, &status); + fits_close_file(fptr, &status); + + #else // all MPI tasks read together: parallel filesystem required pFilereal = fopen (filename,"rb"); long global_index = rank*(xaxis*yaxis)*sizeof(double); fseek(pFilereal, global_index, SEEK_SET); fread(image_real, xaxis*yaxis, sizeof(double), pFilereal); fclose(pFilereal); + #endif // image read // We read binary images, however we can easily extend to fits files @@ -392,12 +444,18 @@ if(rank == 0){ fftwindex = 2*(fftwindex2D + iw*xaxis*yaxis); grid[fftwindex] = fftwgrid[fftwindex2D][0]; grid[fftwindex+1] = fftwgrid[fftwindex2D][1]; + #ifdef WRITE_DATA + grid_real[fftwindex2D] = fftwgrid[fftwindex2D][0]; + grid_img[fftwindex2D] = fftwgrid[fftwindex2D][1]; + #endif } } } fftw_free(fftwgrid); fftw_free(fftwimage); + +/* // This may be moved inside the WRITE_DATA section: TO BE CHECKED // Create sector grid double * gridss; @@ -413,6 +471,7 @@ if(rank == 0){ #ifndef USE_MPI double * gridtot = (double*) calloc(2*grid_size_x*grid_size_y*num_w_planes,sizeof(double)); #endif +*/ // Open the MPI Memory Window for the slab #ifdef USE_MPI @@ -423,76 +482,44 @@ if(rank == 0){ #ifdef WRITE_DATA // Write results - if (rank == 0) + if(rank == 0) { - printf("WRITING GRIDDED DATA\n"); - pFilereal = fopen (outfile2,"wb"); - pFileimg = fopen (outfile3,"wb"); - #ifdef USE_MPI - for (int isector=0; isector 1) - { - for (int iw=0; iw