| Line 7... |
Line 7... |
*
|
*
|
* Author: E.BERTIN (IAP)
|
* Author: E.BERTIN (IAP)
|
*
|
*
|
* Contents: Stuff related to Principal Component Analysis (PCA).
|
* Contents: Stuff related to Principal Component Analysis (PCA).
|
*
|
*
|
* Last modify: 11/10/2007
|
* Last modify: 13/09/2009
|
*
|
*
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
*/
|
*/
|
|
|
#ifdef HAVE_CONFIG_H
|
#ifdef HAVE_CONFIG_H
|
| Line 170... |
Line 170... |
sprintf(str, "CODE_M%d", i+1);
|
sprintf(str, "CODE_M%d", i+1);
|
fitsread(head, str, &code->parammod[i], H_INT, T_LONG);
|
fitsread(head, str, &code->parammod[i], H_INT, T_LONG);
|
}
|
}
|
}
|
}
|
|
|
QMALLOC(pc->maskcurr, double, pc->masksize[0]*pc->masksize[1]*pc->npc);
|
QMALLOC(pc->maskcurr, float, pc->masksize[0]*pc->masksize[1]*pc->npc);
|
|
|
/* But don't touch my arrays!! */
|
/* But don't touch my arrays!! */
|
blank_keys(tab);
|
blank_keys(tab);
|
|
|
return pc;
|
return pc;
|
| Line 187... |
Line 187... |
|
|
/********************************** pc_fit **********************************/
|
/********************************** pc_fit **********************************/
|
/*
|
/*
|
Fit the PC data to the current data.
|
Fit the PC data to the current data.
|
*/
|
*/
|
void pc_fit(psfstruct *psf, double *data, double *weight,
|
void pc_fit(psfstruct *psf, float *data, float *weight,
|
int width, int height,int ix, int iy,
|
int width, int height,int ix, int iy,
|
double dx, double dy, int npc, float backrms)
|
float dx, float dy, int npc, float backrms)
|
{
|
{
|
pcstruct *pc;
|
pcstruct *pc;
|
checkstruct *check;
|
checkstruct *check;
|
codestruct *code;
|
codestruct *code;
|
double *basis,*basis0, *cpix,*cpix0, *pcshift,*wpcshift,
|
double *basis,*basis0;
|
|
float *cpix,*cpix0, *pcshift,*wpcshift,
|
*spix,*wspix, *w, *sumopc,*sumopct, *checkbuf,
|
*spix,*wspix, *w, *sumopc,*sumopct, *checkbuf,
|
*sol,*solt, *datat,
|
*sol,*solt, *datat,
|
*mx2t, *my2t, *mxyt,
|
*mx2t, *my2t, *mxyt,
|
val,val2, xm2,ym2,xym,flux, temp,temp2, theta, pmx2,pmy2,
|
val,val2, xm2,ym2,xym,flux, temp,temp2, theta, pmx2,pmy2,
|
wnorm, ellip, norm, snorm;
|
wnorm, ellip, norm, snorm;
|
| Line 216... |
Line 217... |
ncoeff = psf->poly->ncoeff;
|
ncoeff = psf->poly->ncoeff;
|
pixstep = 1.0/psf->pixstep;
|
pixstep = 1.0/psf->pixstep;
|
dx *= pixstep;
|
dx *= pixstep;
|
dy *= pixstep;
|
dy *= pixstep;
|
|
|
memset(pc->maskcurr, 0, npix*npc*sizeof(double));
|
memset(pc->maskcurr, 0, npix*npc*sizeof(float));
|
basis0 = psf->poly->basis;
|
basis0 = psf->poly->basis;
|
cpix0 = pc->maskcurr;
|
cpix0 = pc->maskcurr;
|
ppix = pc->maskcomp;
|
ppix = pc->maskcomp;
|
|
|
/* Sum each (PSF-dependent) component */
|
/* Sum each (PSF-dependent) component */
|
| Line 230... |
Line 231... |
for (n = ncoeff; n--;)
|
for (n = ncoeff; n--;)
|
{
|
{
|
cpix = cpix0;
|
cpix = cpix0;
|
val = *(basis++);
|
val = *(basis++);
|
for (p=npix; p--;)
|
for (p=npix; p--;)
|
*(cpix++) += val*(double)*(ppix++);
|
*(cpix++) += val**(ppix++);
|
}
|
}
|
}
|
}
|
|
|
/* Allocate memory for temporary buffers */
|
/* Allocate memory for temporary buffers */
|
QMALLOC(pcshift, double, npix2*npc);
|
QMALLOC(pcshift, float, npix2*npc);
|
QMALLOC(wpcshift, double, npix2*npc);
|
QMALLOC(wpcshift, float, npix2*npc);
|
QMALLOC(sol, double, npc);
|
QMALLOC(sol, float, npc);
|
|
|
/* Now shift and scale to the right position, and weight the PCs */
|
/* Now shift and scale to the right position, and weight the PCs */
|
cpix = pc->maskcurr;
|
cpix = pc->maskcurr;
|
spix = pcshift;
|
spix = pcshift;
|
wspix = wpcshift;
|
wspix = wpcshift;
|
| Line 444... |
Line 445... |
}
|
}
|
if ((check = prefs.check[CHECK_PCOPROTOS]))
|
if ((check = prefs.check[CHECK_PCOPROTOS]))
|
{
|
{
|
/*- Reconstruct the unconvolved profile */
|
/*- Reconstruct the unconvolved profile */
|
nopix = pc->omasksize[0]*pc->omasksize[1];
|
nopix = pc->omasksize[0]*pc->omasksize[1];
|
QCALLOC(sumopc, double, nopix);
|
QCALLOC(sumopc, float, nopix);
|
solt = sol;
|
solt = sol;
|
ospix = pc->omaskcomp;
|
ospix = pc->omaskcomp;
|
for (c=npc; c--;)
|
for (c=npc; c--;)
|
{
|
{
|
val = *(solt++);
|
val = *(solt++);
|
sumopct = sumopc;
|
sumopct = sumopc;
|
for (p=nopix; p--;)
|
for (p=nopix; p--;)
|
*(sumopct++) += val*(double)*(ospix++);
|
*(sumopct++) += val**(ospix++);
|
}
|
}
|
QMALLOC(checkbuf, double, npix2);
|
QMALLOC(checkbuf, float, npix2);
|
vignet_resample(sumopc, pc->omasksize[0], pc->omasksize[1],
|
vignet_resample(sumopc, pc->omasksize[0], pc->omasksize[1],
|
checkbuf, width, height, -dx, -dy, pixstep);
|
checkbuf, width, height, -dx, -dy, pixstep);
|
ppix = checkmask;
|
ppix = checkmask;
|
spix = checkbuf;
|
spix = checkbuf;
|
for (p=npix2; p--;)
|
for (p=npix2; p--;)
|
