Skip to content
Snippets Groups Projects
Commit fd1e7301 authored by Stefano Zibetti's avatar Stefano Zibetti
Browse files

Upload New File

parent 172b8ad6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 821 additions and 0 deletions
adaptsmooth.c 0 → 100644
+ 821
− 0
View file @ fd1e7301
#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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment