diff --git a/adaptsmooth.c b/adaptsmooth.c
new file mode 100644
index 0000000000000000000000000000000000000000..8f40977813c1548bc4af7fe0b8bcb782966b9e98
--- /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;
+}