From fd1e730130d4c6b8471c09b2e6aa3059be746bd9 Mon Sep 17 00:00:00 2001
From: Stefano Zibetti <stefano.zibetti@inaf.it>
Date: Tue, 27 Jul 2021 10:18:21 +0000
Subject: [PATCH] Upload New File

---
 adaptsmooth.c | 821 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 821 insertions(+)
 create mode 100644 adaptsmooth.c

diff --git a/adaptsmooth.c b/adaptsmooth.c
new file mode 100644
index 0000000..8f40977
--- /dev/null
+++ b/adaptsmooth.c
@@ -0,0 +1,821 @@
+#define VERSION "V.2.7.1 15.11.2010 by Stefano Zibetti (MPIA&DARK)"
+
+#include "adaptsmooth.h"
+
+void usage()
+{
+    fprintf (stderr, "Adaptsmooth %s\n",VERSION);
+    fprintf (stderr, "Usage: adaptsmooth -p|b|l|m [-a] [-r RMS_SKY] [-G EFFECTIVE_GAIN] [-s SN] [-L L] [-c SLCUT] input_image output_image mask\n");
+    fprintf (stderr, "       -p : noise mode: poisson+background\n");
+    fprintf (stderr, "       -b : noise mode: background-dominated\n");
+    fprintf (stderr, "       -l : noise mode: local\n");
+    fprintf (stderr, "       -m : switch to input-mask mode\n");
+    fprintf (stderr, "       -a : use arithmetic mean instead of median\n");
+    fprintf (stderr, "       -r RMS_SKY (requested in poisson+bkg and bkg-dominated noise mode)\n");
+    fprintf (stderr, "       -G EFFECTIVE_GAIN : e- per ADU (requested in poisson+bkg mode)\n");
+    fprintf (stderr, "       -s SN : minimum S/N requested (defaults to %4.1f)\n",SNMIN);
+    fprintf (stderr, "       -L L : maximum smoothing radius (defaults to max=%2d)\n",MAXRING);
+    fprintf (stderr, "       -c SLCUT : smoothing level cut\n\
+          0: no exclusions; 1: exclude first level (1pix scale) from all others;\n\
+          N<0: exclude up to the (-N)th level below current. Defaults to %d\n",SLCUT);
+
+
+}
+
+
+int main(int argc, char *argv[])
+{
+    float skyRMS=-1.;
+    char ima1_name[FLEN_FILENAME],ima1out_name[FLEN_FILENAME],mask_name[FLEN_FILENAME];
+    fitsfile *imaptr,*imaptr1,*maskptr;
+    int status = 0,anynul,bitpix;
+    long ii,jj;
+    float **image1, **image1_out;
+    int **mask;
+    float nulval=NULVAL;
+    int c;
+    long naxesm[2];
+    short int mode_flag=-1; /* 0:bkg ; 1:local ; 2: poisson ; -1 indef */
+    short int avg_flag=0; /* 0=median is default ; 1=mean */
+
+
+/*********** INITIALIZATION **********/
+/* parse arguments */ 
+/* set default parameters */
+    opterr = 0;
+    fsmooth=&medsmooth_local;
+
+    while ((c = getopt (argc, argv, ":mblpar:s:c:L:G:")) != -1)
+	switch (c)
+	{
+	    case 'm':
+		maskflag = 1;
+		fsmooth=&medsmooth_mask;
+		break;
+	    case 'b':
+		fsmooth=&medsmooth_bkg;
+		mode_flag=0;
+		break;
+	    case 'l':
+		fsmooth=&medsmooth_local;
+		mode_flag=1;
+		break;
+	    case 'p':
+		fsmooth=&medsmooth_poisson;
+		mode_flag=2;
+		break;
+	    case 'a':
+		avg_flag=1;
+		break;
+	    case 'r':
+		skyRMS = atof(optarg);
+		break;
+	    case 's':
+		sn = atof(optarg);
+		break;
+	    case 'c':
+		slconst = atoi(optarg);
+		break;
+	    case 'L':
+		maxsmooth = atoi(optarg);
+		if (maxsmooth>MAXRING) {
+		  fprintf (stderr, "Max smoothing radius too large (%d). adaptsmooth will be run with the maximum radius %d\n", maxsmooth, MAXRING);
+		  maxsmooth=MAXRING;
+		}
+		if (maxsmooth<1) {
+		  fprintf (stderr,"Max smoothing radius < 1 makes no sense... Returning.\n");
+		  usage();
+		  return 7;
+		}
+
+		break;
+	    case 'G':
+		inv_e_gain = 1./atof(optarg);
+		break;
+	    case '?':
+		if (isprint (optopt))
+		    fprintf (stderr, "Unknown option `-%c'.\n", optopt);
+		else
+		  fprintf (stderr,
+			   "Unknown option character `\\x%x'.\n",
+			   optopt);
+		usage();
+		return 1;
+	    case ':':
+		fprintf (stderr, "Option `-%c' requires an argument!\n", optopt);
+		usage();
+		return 2;
+	    default:
+		abort ();
+	}
+    if (mode_flag==-1 && maskflag==0) {
+      fprintf (stderr, "One of the options p|b|l|m must be selected!\n");
+      usage();
+      return 3;
+    }
+    if (avg_flag) {
+      if (maskflag) {
+	fsmooth=&meansmooth_mask;
+      } else {
+	switch (mode_flag)
+	  {
+	  case 0:
+	    fsmooth=&meansmooth_bkg;
+	    break;
+	  case 1:
+	    fsmooth=&meansmooth_local;
+	    break;
+	  case 2:
+	    fsmooth=&meansmooth_poisson;
+	    break;
+	  }
+      }
+    }
+    if (argc-optind!=3) {
+	fprintf (stderr, "Wrong number of arguments given (%d)!\n",argc-optind);
+	usage();
+	return 4;
+    }
+    strcpy(ima1_name,argv[optind]);
+    strcpy(ima1out_name,argv[optind+1]);
+    strcpy(mask_name,argv[optind+2]);
+
+    if (skyRMS==-1. && !maskflag && (mode_flag != 1)) {
+	fprintf (stderr, "Option `-r RMS' must be provided if no smoothing mask is provided nor\
+ local rms computation is requested!\n");
+	usage();
+	return 5;
+    }
+    if (inv_e_gain==-1. && !maskflag && (mode_flag == 2)) {
+	fprintf (stderr, "Option `-G EFFECTIVE_GAIN' must be provided if poisson-noise mode is requested!\n");
+	usage();
+	return 6;
+    }
+    if (slconst<=0)
+	slcut=&slcut_absolute;
+    else 
+	slcut=&slcut_relative;
+
+/***************/
+
+/* open image file */
+    if(fits_open_file(&imaptr,ima1_name,READONLY,&status)) {
+	fits_report_error(stderr,status);
+	return status;
+    };
+/* check image format (bitpix) */
+    if(fits_get_img_type(imaptr,&bitpix,&status)) {
+	fprintf(stderr,"Unreadable image numeric format\nExiting to system now.\n");
+	fits_close_file(imaptr,&status);
+	return status;
+    };
+ 
+    if (bitpix!=FLOAT_IMG ) {
+	fprintf(stderr,"Warning: image numeric format (%d) will be converted to -32 (FLOAT_IMG).\n",bitpix);
+    };
+    
+/* get dimensions and allocate memory */
+    if(fits_get_img_size(imaptr,2,naxes1,&status)) {
+	fits_report_error(stderr,status);
+	fits_close_file(imaptr,&status);
+	return status;
+    };
+/*  printf("%d %d \n",naxes[0],naxes[1]); */
+/* input image 1 */
+    image1=(float **)malloc(naxes1[1]*sizeof(float*));
+    image1[0]=(float *)malloc(naxes1[0]*naxes1[1]*sizeof(float));
+    for (ii=1;ii<naxes1[1];ii++)
+	image1[ii]=image1[ii-1]+naxes1[0];
+/* output image 1 */
+    image1_out=(float **)malloc(naxes1[1]*sizeof(float*));
+    image1_out[0]=(float *)malloc(naxes1[0]*naxes1[1]*sizeof(float));
+    for (ii=1;ii<naxes1[1];ii++)
+	image1_out[ii]=image1_out[ii-1]+naxes1[0];
+/* mask image containing the smooth level */
+    mask=(int **)malloc(naxes1[1]*sizeof(int*));
+    mask[0]=(int *)malloc(naxes1[0]*naxes1[1]*sizeof(int));
+    for (ii=1;ii<naxes1[1];ii++)
+	mask[ii]=mask[ii-1]+naxes1[0];
+ 
+/* load image */
+    if(fits_read_img(imaptr,TFLOAT,1,naxes1[0]*naxes1[1],&nulval,image1[0],&anynul,&status)) {
+	fits_report_error(stderr,status);
+	fits_close_file(imaptr,&status);
+	goto exit_onima;
+    };
+/* do not close input file yet: header must be copied to the output first */
+
+/* initialize output image */
+    for (ii=0;ii<naxes1[0];ii++)
+	for (jj=0;jj<naxes1[1];jj++) {
+	    image1_out[jj][ii]=nulval;
+	}
+ 
+   if (maskflag) {
+/* read in smooth mask */
+/* open mask file */
+       if(fits_open_file(&maskptr,mask_name,READONLY,&status)) {
+	   fits_report_error(stderr,status);
+	   return status;
+       };
+/* check image format (bitpix) */
+       if(fits_get_img_type(maskptr,&bitpix,&status)) {
+	   fprintf(stderr,"Unreadable image numeric format\nExiting to system now.\n");
+	   fits_close_file(maskptr,&status);
+	   return status;
+       };
+       
+       if (bitpix<0 ) { /* if not an integer type */
+	   fprintf(stderr,"Wrong image numeric format. FLOAT (%d) given, INTEGER requested\n Exiting to\
+system now.\n",bitpix);
+	   fits_close_file(maskptr,&status);
+	   return status;
+       };
+       
+/* check get dimensions  */
+       if(fits_get_img_size(maskptr,2,naxesm,&status)) {
+	   fits_report_error(stderr,status);
+	   fits_close_file(maskptr,&status);
+	   return status;
+       };
+       if (naxesm[0]!=naxes1[0] || naxesm[1]!=naxes1[1]) {
+	   fprintf(stderr,"Image and mask dimensions mismatch. Exiting to system...\n");
+	   fits_close_file(maskptr,&status);
+	   goto exit_onima;
+       }
+/* read in the mask */
+       if(fits_read_img(maskptr,TINT,1,naxesm[0]*naxesm[1],&nulval,mask[0],&anynul,&status)) {
+	   fits_report_error(stderr,status);
+	   fits_close_file(maskptr,&status);
+	   goto exit_onima;
+       };
+   } else {
+/* initialize mask to default MSKMAX */
+    for (ii=0;ii<naxes1[0];ii++)
+	for (jj=0;jj<naxes1[1];jj++)
+	    mask[jj][ii]=MSKMAX;
+   }
+   
+   
+/* create smoothed image */
+    for (ii=0;ii<naxes1[0];ii++)
+	for (jj=0;jj<naxes1[1];jj++)
+	  fsmooth(ii,jj,image1,mask,skyRMS,image1_out);
+
+    if (slconst!=0)
+      /* make a second smooth pass to take all slevel cuts into account properly */
+      for (ii=0;ii<naxes1[0];ii++)
+	for (jj=0;jj<naxes1[1];jj++)
+	  fsmooth(ii,jj,image1,mask,skyRMS,image1_out);
+
+/* write output image */
+    long fpix[2]={1,1};
+    if (fits_create_file(&imaptr1,ima1out_name,&status)) {
+	fits_report_error(stderr,status);
+	goto exit_onima;
+    };
+
+/* copy input header to the output */
+    if (fits_copy_header (imaptr,imaptr1, &status)) {
+	fits_report_error(stderr,status);
+	fits_close_file(imaptr,&status);
+	goto exit_onima;
+    };
+    fits_modify_key_lng(imaptr1,(char *)"BITPIX",-32,(char *)"Bits per pixel",&status);
+    fits_write_pix(imaptr1,TFLOAT,fpix,naxes1[0]*naxes1[1],image1_out[0],&status);
+/* update header */
+    fits_write_history(imaptr1,(char *)"Adaptively smoothed image created by ADAPTSMOOTH",&status);
+    fits_write_history(imaptr1,(char *)VERSION,&status);
+    fits_update_key_str(imaptr1,(char *)"IN_PIX",ima1_name,(char *)"Original unsmoothed image",&status);
+    fits_update_key_str(imaptr1,(char *)"MASK",mask_name,(char *)"Smoothing mask",&status);
+    fits_update_key_log(imaptr1,(char *)"MASKMODE",maskflag,(char *)"T=input mask, F=output mask",&status);
+    if (!maskflag)
+      fits_update_key_flt(imaptr1,(char *)"SN",sn,4,(char *)"Minimum S/N",&status);  
+    if (avg_flag) 
+      fits_update_key_str(imaptr1,(char *)"AVERAGE",(char *)"MEAN",(char *)"Average estimator",&status);  
+    else
+      fits_update_key_str(imaptr1,(char *)"AVERAGE",(char *)"MEDIAN",(char *)"Average estimator",&status);  
+    if (!maskflag) {
+      switch (mode_flag) {
+      case 0:
+	fits_update_key_str(imaptr1,(char *)"NOISE_M",(char *)"BKG",(char *)"Background-dominated noise estimates",&status);
+	fits_update_key_flt(imaptr1,(char *)"BKGRMS",skyRMS,4,(char *)"Background noise RMS",&status);    
+	break;
+      case 1:
+	fits_update_key_str(imaptr1,(char *)"NOISE_M",(char *)"LOCAL",(char *)"Noise estimated locally",&status);
+	break;
+      case 2:
+	fits_update_key_str(imaptr1,(char *)"NOISE_M",(char *)"POISSON",(char *)"Poisson+background noise estimates",&status);
+	fits_update_key_flt(imaptr1,(char *)"BKGRMS",skyRMS,4,(char *)"Background noise RMS",&status);    
+	fits_update_key_flt(imaptr1,(char *)"E_GAIN",1./inv_e_gain,4,(char *)"Effective gain (e-/ADU)",&status);    
+	break;
+      }
+    }
+    fits_update_key_lng(imaptr1,(char *)"SLCONST",slconst,(char *)"Smoothing level cut",&status);
+/* close file */
+    fits_close_file(imaptr1,&status);
+    
+
+    if (maskflag==0) {
+/* write output mask */
+	if (fits_create_file(&maskptr,mask_name,&status)) {
+	    fits_report_error(stderr,status);
+	    goto exit_onima;
+	};
+/* copy input header to the output */
+    if (fits_copy_header (imaptr,maskptr, &status)) {
+	fits_report_error(stderr,status);
+	fits_close_file(imaptr,&status);
+	goto exit_onima;
+    };
+	fits_modify_key_lng(maskptr,(char *)"BITPIX",16,(char *)"Bits per pixel",&status);    
+	/* fits_create_img(maskptr,SHORT_IMG,2,naxes1,&status); */
+	fits_write_pix(maskptr,TINT,fpix,naxes1[0]*naxes1[1],mask[0],&status);
+/* update header */
+	fits_write_history(maskptr,(char *)"Smoothing mask created by ADAPTSMOOTH",&status);
+	fits_write_history(maskptr,(char *)VERSION,&status);
+	fits_update_key_str(maskptr,(char *)"IN_PIX",ima1_name,(char *)"Original unsmoothed image",&status);
+	fits_update_key_str(maskptr,(char *)"OUT_PIX",ima1out_name,(char *)"Output smoothed image",&status);
+	fits_update_key_str(maskptr,(char *)"MASK",mask_name,(char *)"Smoothing mask",&status);
+	fits_update_key_log(maskptr,(char *)"MASKMODE",maskflag,(char *)"T=input mask, F=output mask",&status);
+	fits_update_key_flt(maskptr,(char *)"SN",sn,4,(char *)"Minimum S/N",&status);  
+	if (avg_flag) 
+	  fits_update_key_str(maskptr,(char *)"AVERAGE",(char *)"MEAN",(char *)"Average estimator",&status);  
+	else
+	  fits_update_key_str(maskptr,(char *)"AVERAGE",(char *)"MEDIAN",(char *)"Average estimator",&status);  
+	if (!maskflag) {
+	  switch (mode_flag) {
+	  case 0:
+	    fits_update_key_str(maskptr,(char *)"NOISE_M",(char *)"BKG",(char *)"Background-dominated noise estimates",&status);
+	    fits_update_key_flt(maskptr,(char *)"BKGRMS",skyRMS,4,(char *)"Background noise RMS",&status);    
+	    break;
+	  case 1:
+	    fits_update_key_str(maskptr,(char *)"NOISE_M",(char *)"LOCAL",(char *)"Noise estimated locally",&status);
+	    break;
+	  case 2:
+	    fits_update_key_str(maskptr,(char *)"NOISE_M",(char *)"POISSON",(char *)"Poisson+background noise estimates",&status);
+	    fits_update_key_flt(maskptr,(char *)"BKGRMS",skyRMS,4,(char *)"Background noise RMS",&status);    
+	    fits_update_key_flt(maskptr,(char *)"E_GAIN",1./inv_e_gain,4,(char *)"Effective gain (e-/ADU)",&status);    
+	    break;
+	  }
+	}
+	fits_update_key_lng(maskptr,(char *)"SLCONST",slconst,(char *)"Smoothing level cut",&status);
+/* close file */
+	fits_close_file(maskptr,&status);
+/* close input file */
+    fits_close_file(imaptr,&status);
+    }
+ exit_onima:
+    free(image1[0]);
+    free(image1);
+    free(image1_out[0]);
+    free(image1_out);
+    free(mask[0]);
+    free(mask);
+    return status;
+    
+}
+
+
+int slcut_relative(int sl)
+{
+    extern int slconst;
+    return sl+slconst;
+}
+
+int slcut_absolute(int sl)
+{
+    extern int slconst;
+    return slconst;
+}
+
+int medsmooth_mask (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN];
+    int sl,i,ibad,xp,yp;
+    unsigned int n_el;
+
+    ibad=0;
+    if (mask[y][x]<0 || mask[y][x]>=MSKMAX) /* bad pixel or SN too low */
+      return MSKMAX; /* nothing to do */
+    sl=mask[y][x];
+    if (sl==1) {
+      outimage[y][x]=image[y][x];
+      return sl;
+    }
+    for (i=0;i<kernel_limits[sl-1];i++) {
+      yp=y-kernel_pix[i][1];
+      if (yp<0) yp=0;
+      if (yp>=naxes1[1]) yp=naxes1[1]-1;
+      xp=x-kernel_pix[i][0];
+      if (xp<0) xp=0;
+      if (xp>=naxes1[0]) xp=naxes1[0]-1;
+      if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the median subsample */
+	pix_values[i-ibad]=image[yp][xp];
+      else ibad++;
+    }
+    n_el=kernel_limits[sl-1]-ibad;
+    outimage[y][x]=select_FR(pix_values,0,n_el-1,n_el >> 1);
+    return sl;
+}
+
+int medsmooth_local (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN], pix_values_sigma[MAXKERN];
+    float median,sigma;
+    int sl,sl_sigma,i,ibad,ibad_sigma,xp,yp;
+    unsigned int n_el,n_el_sigma;
+
+    ibad=0;
+    ibad_sigma=0;
+    /* determine smoothing level from image */
+    sl=1;
+    sl_sigma=sl+DELTA_SL_RMS;
+    if (sl_sigma>=maxsmooth) sl_sigma=maxsmooth-1;
+    /* copy pixel values into the arrays to be sorted */
+    for (i=0;i<kernel_limits[sl_sigma-1];i++) {
+      yp=y-kernel_pix[i][1];
+      if (yp<0) yp=0;
+      if (yp>=naxes1[1]) yp=naxes1[1]-1;
+      xp=x-kernel_pix[i][0];
+      if (xp<0) xp=0;
+      if (xp>=naxes1[0]) xp=naxes1[0]-1;
+      if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the sigma subsample */
+	pix_values_sigma[i-ibad_sigma]=image[yp][xp];
+      else ibad_sigma++;
+    }
+    n_el_sigma=kernel_limits[sl_sigma-1]-ibad_sigma;
+    sigma=(select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.84)) - \
+	   select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.16)))/2.;
+    if (image[y][x]>sn*sigma) {
+      outimage[y][x]=image[y][x];
+      mask[y][x]=1;
+      return sl; /* everything fine, no need to smooth */
+    }
+    sl++;
+    pix_values[0]=image[y][x];
+    while (sl<=maxsmooth) {
+      /* copy pixel values into the arrays to be sorted */
+      for (i=kernel_limits[sl-2];i<kernel_limits[sl-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the median subsample */
+	  pix_values[i-ibad]=image[yp][xp];
+	else ibad++;
+      }
+      sl_sigma=sl+DELTA_SL_RMS;
+      if (sl_sigma>=maxsmooth) sl_sigma=maxsmooth-1;
+      for (i=kernel_limits[sl_sigma-2];i<kernel_limits[sl_sigma-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the sigma subsample */
+	  pix_values_sigma[i-ibad_sigma]=image[yp][xp];
+	else ibad_sigma++;
+      }
+      n_el=kernel_limits[sl-1]-ibad;
+      median=select_FR(pix_values,0,n_el-1,n_el >> 1);
+      n_el_sigma=kernel_limits[sl_sigma-1]-ibad_sigma;
+      sigma=(select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.84)) - \
+	     select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.16)))/2.;
+      if (median > sn * sigma / sqrt(kernel_limits[sl-1]-ibad) ) {
+	outimage[y][x]=median;
+	mask[y][x]=sl;
+	return sl;
+      }
+      sl++;
+    }
+    /* minimum S/N cannot be reached */
+    outimage[y][x]=NULVAL;
+    mask[y][x]=MSKMAX;
+    return MSKMAX;
+}
+
+int medsmooth_bkg (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN];
+    float median;
+    int sl,i,ibad,xp,yp;
+    unsigned int n_el;
+
+    ibad=0;
+    /* determine smoothing level from image */
+    sl=1;
+    if (image[y][x]>sn*rms) {
+      outimage[y][x]=image[y][x];
+      mask[y][x]=1;
+      return sl; /* everything fine, no need to smooth */
+    }
+    sl++;
+    pix_values[0]=image[y][x];
+    while (sl<=maxsmooth) {
+      /* copy pixel values into the array to be sorted */
+      for (i=kernel_limits[sl-2];i<kernel_limits[sl-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the median subsample */
+	  pix_values[i-ibad]=image[yp][xp];
+	else ibad++;
+      }
+      n_el=kernel_limits[sl-1]-ibad;
+      median=select_FR(pix_values,0,n_el-1,n_el >> 1);
+      if (median > sn * rms / sqrt(kernel_limits[sl-1]-ibad) ) {
+	outimage[y][x]=median;
+	mask[y][x]=sl;
+	return sl;
+      }
+      sl++;
+    }
+    
+    /* minimum S/N cannot be reached */
+    outimage[y][x]=NULVAL;
+    mask[y][x]=MSKMAX;    
+    return MSKMAX; 
+}
+
+int medsmooth_poisson (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN];
+    float median;
+    int sl,i,ibad,xp,yp;
+    unsigned int n_el;
+
+    ibad=0;
+    /* determine smoothing level from image */
+    sl=1;
+    /* 1/sn = sqrt ((rms/i)^2+1/(e_gain*i)) */
+    median=image[y][x];
+    if (median > 0)
+      if ((rms*rms/median+inv_e_gain) < median/(sn*sn)) {
+	outimage[y][x]=median;
+	mask[y][x]=1;
+	return sl; /* everything fine, no need to smooth */
+      }
+    sl++;
+    pix_values[0]=image[y][x];
+    while (sl<=maxsmooth) {
+      /* copy pixel values into the array to be sorted */
+      for (i=kernel_limits[sl-2];i<kernel_limits[sl-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the median subsample */
+	  pix_values[i-ibad]=image[yp][xp];
+	else ibad++;
+      }
+      n_el=kernel_limits[sl-1]-ibad;
+      median=select_FR(pix_values,0,n_el-1,n_el >> 1);
+      if (median > 0)
+	if ((rms*rms/median+inv_e_gain) < median/(sn*sn/(kernel_limits[sl-1]-ibad))) {
+	  outimage[y][x]=median;
+	  mask[y][x]=sl;
+	  return sl;
+	}
+      sl++;
+    }
+    
+    /* minimum S/N cannot be reached */
+    outimage[y][x]=NULVAL;
+    mask[y][x]=MSKMAX;    
+    return MSKMAX; 
+}
+
+int meansmooth_mask (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN];
+    int sl,i,ibad,xp,yp;
+    float mean;
+    unsigned int n_el;
+    
+    ibad=0;
+    if (mask[y][x]<0 || mask[y][x]>=MSKMAX) /* bad pixel or SN too low */
+      return MSKMAX; /* nothing to do */
+    sl=mask[y][x];
+    if (sl==1) {
+      outimage[y][x]=image[y][x];
+      return sl;
+    }
+    for (i=0;i<kernel_limits[sl-1];i++) {
+      yp=y-kernel_pix[i][1];
+      if (yp<0) yp=0;
+      if (yp>=naxes1[1]) yp=naxes1[1]-1;
+      xp=x-kernel_pix[i][0];
+      if (xp<0) xp=0;
+      if (xp>=naxes1[0]) xp=naxes1[0]-1;
+      if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the mean subsample */
+	pix_values[i-ibad]=image[yp][xp];
+      else ibad++;
+    }
+    n_el=kernel_limits[sl-1]-ibad;
+    for (i=0,mean=0.0;i<n_el;i++) mean+=pix_values[i];
+    outimage[y][x]=mean/n_el;
+    return sl;
+}
+
+int meansmooth_local (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN], pix_values_sigma[MAXKERN];
+    float mean,sigma;
+    int sl,sl_sigma,i,ibad,ibad_sigma,xp,yp;
+    unsigned int n_el,n_el_sigma;
+
+    ibad=0;
+    ibad_sigma=0;
+    /* determine smoothing level from image */
+    sl=1;
+    sl_sigma=sl+DELTA_SL_RMS;
+    if (sl_sigma>=maxsmooth) sl_sigma=maxsmooth-1;
+    /* copy pixel values into the arrays to be sorted */
+    for (i=0;i<kernel_limits[sl_sigma-1];i++) {
+      yp=y-kernel_pix[i][1];
+      if (yp<0) yp=0;
+      if (yp>=naxes1[1]) yp=naxes1[1]-1;
+      xp=x-kernel_pix[i][0];
+      if (xp<0) xp=0;
+      if (xp>=naxes1[0]) xp=naxes1[0]-1;
+      if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the sigma subsample */
+	pix_values_sigma[i-ibad_sigma]=image[yp][xp];
+      else ibad_sigma++;
+    }
+    n_el_sigma=kernel_limits[sl_sigma-1]-ibad_sigma;
+    sigma=(select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.84)) - \
+	   select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.16)))/2.;
+    if (image[y][x]>sn*sigma) {
+      outimage[y][x]=image[y][x];
+      mask[y][x]=1;
+      return sl; /* everything fine, no need to smooth */
+    }
+    sl++;
+    pix_values[0]=image[y][x];
+    while (sl<=maxsmooth) {
+      /* copy pixel values into the arrays to be sorted */
+      for (i=kernel_limits[sl-2];i<kernel_limits[sl-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the mean subsample */
+	  pix_values[i-ibad]=image[yp][xp];
+	else ibad++;
+      }
+      sl_sigma=sl+DELTA_SL_RMS;
+      if (sl_sigma>=maxsmooth) sl_sigma=maxsmooth-1;
+      for (i=kernel_limits[sl_sigma-2];i<kernel_limits[sl_sigma-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the sigma subsample */
+	  pix_values_sigma[i-ibad_sigma]=image[yp][xp];
+	else ibad_sigma++;
+      }
+      n_el=kernel_limits[sl-1]-ibad;
+      for (i=0,mean=0.0;i<n_el;i++) mean+=pix_values[i];
+      mean/=n_el;
+      n_el_sigma=kernel_limits[sl_sigma-1]-ibad_sigma;
+      sigma=(select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.84)) - \
+	     select_FR(pix_values_sigma,0,n_el_sigma-1, (int)(n_el_sigma * 0.16)))/2.;
+      if (mean > sn * sigma / sqrt(kernel_limits[sl-1]-ibad) ) {
+	outimage[y][x]=mean;
+	mask[y][x]=sl;
+	return sl;
+      }
+      sl++;
+    }
+    /* minimum S/N cannot be reached */
+    outimage[y][x]=NULVAL;
+    mask[y][x]=MSKMAX;
+    return MSKMAX;
+}
+
+int meansmooth_bkg (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN];
+    float mean;
+    int sl,i,ibad,xp,yp;
+    unsigned int n_el;
+
+    ibad=0;
+    /* determine smoothing level from image */
+    sl=1;
+    if (image[y][x]>sn*rms) {
+      outimage[y][x]=image[y][x];
+      mask[y][x]=1;
+      return sl; /* everything fine, no need to smooth */
+    }
+    sl++;
+    pix_values[0]=image[y][x];
+    while (sl<=maxsmooth) {
+      /* copy pixel values into the array to be sorted */
+      for (i=kernel_limits[sl-2];i<kernel_limits[sl-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the mean subsample */
+	  pix_values[i-ibad]=image[yp][xp];
+	else ibad++;
+      }
+      n_el=kernel_limits[sl-1]-ibad;
+      for (i=0,mean=0.0;i<n_el;i++) mean+=pix_values[i];
+      mean/=n_el;
+      if (mean > sn * rms / sqrt(kernel_limits[sl-1]-ibad) ) {
+	outimage[y][x]=mean;
+	mask[y][x]=sl;
+	return sl;
+      }
+      sl++;
+    }
+    
+    /* minimum S/N cannot be reached */
+    outimage[y][x]=NULVAL;
+    mask[y][x]=MSKMAX;    
+    return MSKMAX; 
+}
+
+int meansmooth_poisson (long int x, long int y, float **image, int **mask, float rms, float **outimage )
+{
+    extern int kernel_pix[][2];
+    extern int kernel_limits[];
+    static float pix_values[MAXKERN];
+    float mean;
+    int sl,i,ibad,xp,yp;
+    unsigned int n_el;
+
+    ibad=0;
+    /* determine smoothing level from image */
+    sl=1;
+    /* 1/sn = sqrt ((rms/i)^2+1/(e_gain*i)) */
+    mean=image[y][x];
+    if (mean > 0)
+      if ((rms*rms/mean+inv_e_gain) < mean/(sn*sn)) {
+	outimage[y][x]=mean;
+	mask[y][x]=1;
+	return sl; /* everything fine, no need to smooth */
+      }
+    sl++;
+    pix_values[0]=image[y][x];
+    while (sl<=maxsmooth) {
+      /* copy pixel values into the array to be sorted */
+      for (i=kernel_limits[sl-2];i<kernel_limits[sl-1];i++) {
+	yp=y-kernel_pix[i][1];
+	if (yp<0) yp=0;
+	if (yp>=naxes1[1]) yp=naxes1[1]-1;
+	xp=x-kernel_pix[i][0];
+	if (xp<0) xp=0;
+	if (xp>=naxes1[0]) xp=naxes1[0]-1;
+	if (mask[yp][xp]>slcut(sl)) /* exclude higher S/N pixels from the mean subsample */
+	  pix_values[i-ibad]=image[yp][xp];
+	else ibad++;
+      }
+      n_el=kernel_limits[sl-1]-ibad;
+      for (i=0,mean=0.0;i<n_el;i++) mean+=pix_values[i];
+      mean/=n_el;
+      if (mean > 0)
+	if ((rms*rms/mean+inv_e_gain) < mean/(sn*sn/(kernel_limits[sl-1]-ibad))) {
+	  outimage[y][x]=mean;
+	  mask[y][x]=sl;
+	  return sl;
+	}
+      sl++;
+    }
+    
+    /* minimum S/N cannot be reached */
+    outimage[y][x]=NULVAL;
+    mask[y][x]=MSKMAX;    
+    return MSKMAX; 
+}
-- 
GitLab