Commit 6384c370 authored by Danny Boxhoorn's avatar Danny Boxhoorn
Browse files

Version from 2005/07/07

parents
void apply_norm(float *imdif, unsigned short *flag0, unsigned short *flag1,
double coeff, int nx, int ny, float bad_value)
{
double norm;
int idx;
norm = 1.0/coeff;
for(idx=0;idx<nx*ny;idx++)
if( !flag0[idx] || !flag1[idx] )
//imdif[idx] = bad_value; modifica fatta il 7/10 per ottenere delle immagini differenza sulle immagini di vincenzo 500x500 !!
// imdif[idx] = bad_value;
imdif[idx] *= norm;
else
imdif[idx] *= norm;
return;
}
#include <stdio.h>
#include <math.h>
void base_func(float sigma, int m, int n, double *ker_x, double *ker_y,
int kerhw, int do_norm)
{
double sig, norm_x, norm_y ;
int i, idx;
sig = (double)(2*sigma*sigma);
norm_x = norm_y = 0.0;
for(i=-kerhw;i<=kerhw;i++) {
idx = kerhw + i;
ker_y[idx] = ker_x[idx] = exp(-i*i/sig);
if( m > 0 )
ker_x[idx] *= pow(i,m);
if( n > 0 )
ker_y[idx] *= pow(i,n);
norm_x += ker_x[idx];
norm_y += ker_y[idx];
}
if( do_norm ) {
if( norm_x != 0.0 && norm_y != 0.0 ) {
for(i=0;i<=2*kerhw;i++) {
ker_x[i] /= norm_x;
ker_y[i] /= norm_y;
}
}
else {
fprintf(stderr,"base_func: warning: normalization to 0.0 ignored!\n");
}
return ;
}
return ;
}
#include "defs.h"
int clean_domains(float *im, float *imref, unsigned short *mask0, unsigned short *mask1,
double **vecs, DOM_TABS *domp, PARAMS par)
{
double **mat0, **mat1, *vec0;
int i, j, m, n, idx, npix, xc, yc, idom;
npix = 0;
for(idom=0;idom<par.ndom;idom++) {
mat0 = domp[idom].mat0;
mat1 = domp[idom].mat1;
vec0 = domp[idom].vec0;
xc = domp[idom].xc;
yc = domp[idom].yc;
for(m=0;m<par.nvecs;m++) {
vec0[m] = 0.0;
for(n=0;n<=m;n++)
mat1[m][n] = mat0[m][n];
}
for(i=xc-par.domhw;i<=xc+par.domhw;i++) {
for(j=yc-par.domhw;j<=yc+par.domhw;j++) {
idx = i + par.nx*j;
if( mask0[idx] ) {
if( mask1[idx] ) {
for(m=0;m<par.nvecs;m++)
vec0[m] += im[idx]*vecs[m][idx]*imref[idx];
npix++;
}
else {
for(m=0;m<par.nvecs;m++)
for(n=0;n<=m;n++)
mat1[m][n] -= vecs[m][idx]*vecs[n][idx]*imref[idx];
}
}
}
}
}
/* mask1 becomes logical product of mask0 and mask1 */
for(idx=0;idx<par.nx*par.ny;idx++)
if( !mask0[idx] )
mask1[idx] = 0;
return (npix);
}
#include <math.h>
#include "defs.h"
#include "funcs.h"
int clip_domains(DOM_TABS *domp, float *sort_buf, int *indx,
int iter, double *chi2_n, PARAMS par)
{
int nkeep, idom, n_sig;
double median, sig, sigma;
n_sig = par.n_sig_dom;
nkeep = 0;
for(idom=0;idom<par.ndom;idom++) {
if( domp[idom].reject != RJCT_NONE )
continue;
sort_buf[nkeep++] = domp[idom].chi2_n;
}
quick_sort(sort_buf,indx,nkeep);
median = sort_buf[indx[nkeep/2]];
sigma = 0.0;
for(idom=0;idom<par.ndom;idom++) {
if( domp[idom].reject != RJCT_NONE )
continue;
sig = domp[idom].chi2_n - median;
sigma += sig*sig;
}
sigma /= nkeep;
nkeep = 0;
for(idom=0;idom<par.ndom;idom++) {
sig = domp[idom].chi2_n - median;
if( domp[idom].reject < RJCT_HIGH ) {
domp[idom].reject = RJCT_NONE;
if( sig > n_sig*sqrt(sigma) )
domp[idom].reject = RJCT_LOW;
else
nkeep++;
}
}
*chi2_n = median;
if( par.verbose >= VERB_MED )
printf("%6d\t%14.4f\t%14.4f\t%7d\n",iter,par.gain*median,par.gain*sigma,nkeep);
if( par.verbose >= VERB_HIGH ) {
for(idom=0;idom<par.ndom;idom++) {
if( domp[idom].reject != RJCT_NONE ) {
printf("rejected domain %4d:\tx=%3d, \ty=%3d,",
idom,domp[idom].xc,domp[idom].yc);
printf("\tchi2_n=%6.3f,\tstatus=%1d\n",
domp[idom].chi2_n*par.gain,domp[idom].reject);
}
}
}
return (nkeep);
}
#include "defs.h"
void clone_domains(DOM_TABS *domp, int ndom, int nvecs)
{
double **mat, **mat1, *vec0, *vec;
int m, n, idom;
for(idom=0;idom<ndom;idom++) {
mat = domp[idom].mat;
vec = domp[idom].vec;
mat1 = domp[idom].mat1;
vec0 = domp[idom].vec0;
for(m=0;m<nvecs;m++) {
vec[m] = vec0[m];
for(n=0;n<=m;n++)
mat[m][n] = mat1[m][n];
}
for(m=0;m<nvecs;m++)
for(n=0;n<m;n++)
mat[n][m] = mat[m][n];
}
return;
}
#define MAX_NCOMP 5
#define MAX_NIM 2000
#define MAX_FNAME 128
#define BIG_FLOAT 1.0e32
#define MAX_HEADLEN 2880*32
#define CHI2_INIT 999.0
#define VERB_NONE 0
#define VERB_LOW 1
#define VERB_MED 2
#define VERB_HIGH 3
#define RJCT_NONE 0
#define RJCT_LOW 1
#define RJCT_HIGH 2
typedef struct {
int xc;
int yc;
int reject;
double chi2_n;
double **mat0;
double **mat1;
double **mat;
double *vec0;
double *vec;
} DOM_TABS;
typedef struct {
int nx0;
int ny0;
int nx;
int ny;
int kerhw;
int ncomp;
int nvecs;
int bdeg;
int *deg;
float *sig;
int deg_inc;
float sig_inc;
float sat_level;
float min_level;
float bad_value;
int bad_growrad1;
int bad_growrad2;
float min_area;
float gain;
int n_iter;
int verbose;
float n_sig;
float max_chi2;
int sdeg;
int wdeg;
int nbkg;
int nwxy;
int nspat;
int ntot;
int ndom;
int domhw;
float n_sig_dom;
float dom_thresh;
float min_area_dom;
int ndom_x;
int ndom_y;
int mohw;
int domain_mode;
int n_iter_dom;
int min_nkeep;
} PARAMS;
#include "defs.h"
void expand_matrix(double **exp_mat, double *exp_vec, double **wxy,
DOM_TABS *domp, PARAMS par)
{
double **mat, *vec;
int m, n, ki, kj, li, lj, xc, yc, idom, idx;
for(m=0;m<par.ntot;m++) {
exp_vec[m] = 0.0;
for(n=0;n<=m;n++)
exp_mat[m+1][n+1] = 0.0;
}
for(idom=0;idom<par.ndom;idom++) {
if( domp[idom].reject != RJCT_NONE )
continue;
vec = domp[idom].vec;
mat = domp[idom].mat;
xc = domp[idom].xc;
yc = domp[idom].yc;
idx = xc + par.nx*yc;
for(m=0;m<par.nvecs;m++) {
exp_vec[m] += vec[m];
for(n=0;n<=m;n++)
exp_mat[m+1][n+1] += mat[m][n];
}
m = par.nvecs;
for(ki=1;ki<par.nwxy;ki++) {
for(li=par.nbkg+1;li<par.nvecs;li++) {
exp_vec[m] += wxy[ki][idx]*vec[li];
for(n=0;n<par.nvecs;n++)
exp_mat[m+1][n+1] += wxy[ki][idx]*mat[li][n];
m++;
}
}
m = par.nvecs;
for(ki=1;ki<par.nwxy;ki++) {
for(li=par.nbkg+1;li<par.nvecs;li++) {
n = par.nvecs;
for(kj=1;kj<par.nwxy && n<=m;kj++) {
for(lj=par.nbkg+1;lj<par.nvecs && n<=m;lj++) {
exp_mat[m+1][n+1] += wxy[ki][idx]*wxy[kj][idx]*mat[li][lj];
n++;
}
}
m++;
}
}
}
for(m=0;m<par.ntot;m++)
for(n=0;n<m;n++)
exp_mat[n+1][m+1] = exp_mat[m+1][n+1];
return ;
}
#include <stdio.h>
/* preparation of arrays for image fitting and clipping */
extern void get_params(char *, PARAMS *par);
extern void expand_matrix(double **, double *, double **, DOM_TABS *, PARAMS);
extern void make_vectors(float *, double **, double **, double **, PARAMS);
extern void make_domains(float *, unsigned short *, double **, DOM_TABS *, PARAMS);
extern void clone_domains(DOM_TABS *, int, int);
extern void mask_badpix(float *, unsigned short *, unsigned short *, int, PARAMS);
extern int clean_domains(float *, float *, unsigned short *, unsigned short *,
double **, DOM_TABS *, PARAMS);
extern void get_domains(float *, unsigned short *, DOM_TABS *, PARAMS, int *);
extern int max(int, int, int, float *, int);
/* sigma clipping and local value of fit */
extern void quick_sort(float *, int *, int);
extern int local_clip(float *, float *, unsigned short *, double **,
double **, double *, DOM_TABS *, PARAMS, double *, int *);
extern int clip_domains(DOM_TABS *, float *, int *, int, double*, PARAMS);
extern void local_solution(int, int, double **, double *, double *, PARAMS);
/* basis function for kernel decomposition */
extern double base_func(float, int, int, double *, double *, int, int);
/* FITS manipulation */
extern void init_difimages(char **, char **, long int *, int, PARAMS par);
extern void wfits_f(char *, char *, long int , float *);
extern void hedit(char *, int, char *, char, char *, char *, ...);
extern FILE *read_header(char *, char *, long int *);
extern int read_sector(char *infn, long int, int, int, int, int,
char *, int, int, int);
extern int write_sector(char *infn, long int, int, int, int, int,
char *, int, int, int, int);
/* convolutions */
extern void xy_convolve(float *, double *, double *, int, int,
double *, double *, int);
extern void spatial_convolve(float *, float *, double **, double **,
double *, PARAMS);
/* solving of linear equation (numerical recipes) */
extern void ludcmp(double **, int, int *, double *);
extern void lubksb(double **, int, int *, double *);
/* output */
extern void write_kernel(char *, double *, int, int, int, PARAMS *, double *);
extern void apply_norm(float *, unsigned short *, unsigned short *,
double, int, int, float);
#include <math.h>
#include "defs.h"
#include "funcs.h"
void get_domains(float *im, unsigned short *mask0, DOM_TABS *domp,
PARAMS par, int *ndom)
{
int i, j, nx_sector, ny_sector, xs, ys;
int i_max=0, j_max=0, idx, idom, idom_x, idom_y, nmarg;
float pixmax, pix;
FILE *inpf;
nmarg = par.kerhw + par.domhw;
idom = 0;
if( par.domain_mode==1 ) {
//printf("--------------->> domain normale");
/* domains around bright stars */
nx_sector = (int)ceil((par.nx-2*nmarg)/(double)par.ndom_x);
ny_sector = (int)ceil((par.ny-2*nmarg)/(double)par.ndom_y);
//inpf = fopen("coordinate.txt","w");
for(idom_x=0;idom_x<par.ndom_x;idom_x++) {
xs = nmarg + idom_x*nx_sector;
for(idom_y=0;idom_y<par.ndom_y;idom_y++) {
ys = nmarg + idom_y*ny_sector;
pixmax = 0.0;
for(i=xs;i<xs+nx_sector && i<par.nx-nmarg;i++) {
for(j=ys;j<ys+ny_sector && j<par.ny-nmarg;j++) {
idx = i + par.nx*j;
if( max(i,j,par.mohw,im,par.nx) && mask0[idx] ) {
pix = im[idx];
if( pix > pixmax ) {
pixmax = pix;
i_max = i;
j_max = j;
}
}
}
}
if( pixmax > par.dom_thresh ) {
domp[idom].xc = i_max;
domp[idom].yc = j_max;
domp[idom].reject = RJCT_NONE; /* #define RJCT_NONE 0 */
domp[idom].chi2_n = CHI2_INIT; /* #define CHI2_INIT 999.0 */
idom++;
//fprintf(inpf,"%d %d \n",i_max,j_max);
}
}
}
*ndom = idom;
//close(inpf);
}
else {
if( par.domain_mode==2 ) {
/* domains take in input */
//printf(" NUOVA FUNZIONALITa' prende coordinate da un file");
nx_sector = (par.nx-2*nmarg-1)/(par.ndom_x-1);
ny_sector = (par.ny-2*nmarg-1)/(par.ndom_y-1);
inpf = fopen("coordinate.txt","r");
for(idom_x=0;idom_x<par.ndom_x;idom_x++) {
for(idom_y=0;idom_y<par.ndom_y;idom_y++) {
if ( fscanf(inpf,"%d %d",&xs,&ys) != EOF){
if( xs<par.nx-nmarg && ys<par.ny-nmarg ) {
domp[idom].xc = xs;
domp[idom].yc = ys;
domp[idom].reject = RJCT_NONE;
domp[idom].chi2_n = CHI2_INIT;
idom++;
}
}
}
}
*ndom = idom;
close(inpf);
}
else
{
/* domains uniformly distributed in the image */
nx_sector = (par.nx-2*nmarg-1)/(par.ndom_x-1);
ny_sector = (par.ny-2*nmarg-1)/(par.ndom_y-1);
for(idom_x=0;idom_x<par.ndom_x;idom_x++) {
xs = nmarg + idom_x*nx_sector;
for(idom_y=0;idom_y<par.ndom_y;idom_y++) {
ys = nmarg + idom_y*ny_sector;
if( xs<par.nx-nmarg && ys<par.ny-nmarg ) {
domp[idom].xc = xs;
domp[idom].yc = ys;
domp[idom].reject = RJCT_NONE;
domp[idom].chi2_n = CHI2_INIT;
idom++;
}
}
}
*ndom = idom;
}
}
if( par.verbose > VERB_HIGH ) {
printf("Domains (%d) :\n",*ndom);
for(idom=0;idom<*ndom;idom++)
printf("%3d\tx=%3d,\ty=%3d\n",idom,domp[idom].xc,domp[idom].yc);
}
return;
}
#include <stdio.h>
#include "defs.h"
#include "funcs.h"
void get_params(char *parfname, PARAMS *par)
{
FILE *inpf;
inpf = fopen(parfname,"r");
fscanf(inpf,"%*s = %d %*s = %d",&par->nx0,&par->ny0);