1 /*
2 pc.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: Stuff related to Principal Component Analysis (PCA).
11 *
12 * Last modify: 27/11/2003
13 *
14 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 */
16
17 #ifdef HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25
26 #include "define.h"
27 #include "globals.h"
28 #include "prefs.h"
29 #include "fits/fitscat.h"
30 #include "check.h"
31 #include "image.h"
32 #include "poly.h"
33 #include "psf.h"
34
35 static obj2struct *obj2 = &outobj2;
36
37 /****** pc_end ***************************************************************
38 PROTO void pc_end(pcstruct *pc)
39 PURPOSE Free a PC structure and everything it contains.
40 INPUT pcstruct pointer.
41 OUTPUT -.
42 NOTES -.
43 AUTHOR E. Bertin (IAP, Leiden observatory & ESO)
44 VERSION 15/07/99
45 ***/
pc_end(pcstruct * pc)46 void pc_end(pcstruct *pc)
47 {
48 int i;
49
50 free(pc->maskcomp);
51 free(pc->omaskcomp);
52 free(pc->omasksize);
53 free(pc->maskcurr);
54 free(pc->masksize);
55 free(pc->mx2);
56 free(pc->my2);
57 free(pc->mxy);
58 free(pc->flux);
59 free(pc->bt);
60 if (pc->code)
61 {
62 free(pc->code->pc);
63 for (i=0; i<pc->code->nparam;i++)
64 free(pc->code->param[i]);
65 free(pc->code->param);
66 free(pc->code);
67 }
68 free(pc);
69
70 return;
71 }
72
73
74 /********************************** pc_load **********************************/
75 /*
76 Load the PC data from a FITS file.
77 */
pc_load(catstruct * cat)78 pcstruct *pc_load(catstruct *cat)
79 {
80 pcstruct *pc;
81 tabstruct *tab;
82 keystruct *key;
83 codestruct *code;
84 char *head, str[80], *ci, *filename;
85 int i, ncode,nparam;
86
87 if (!(tab = name_to_tab(cat, "PC_DATA", 0)))
88 return NULL;
89
90 filename = cat->filename;
91
92 /* OK, we now allocate memory for the PC structure itself */
93 QCALLOC(pc, pcstruct, 1);
94
95 /* Store a short copy of the PC filename */
96 if ((ci=strrchr(filename, '/')))
97 strcpy(pc->name, ci+1);
98 else
99 strcpy(pc->name, filename);
100
101 /* Load important scalars (which are stored as FITS keywords) */
102 head = tab->headbuf;
103
104 /* Dimensionality of the PC mask */
105 if (fitsread(head, "PCNAXIS", &pc->maskdim, H_INT, T_LONG) != RETURN_OK)
106 return NULL;
107 if (pc->maskdim<2 || pc->maskdim>4)
108 error(EXIT_FAILURE, "*Error*: wrong dimensionality for the PC "
109 "mask in ", filename);
110 QMALLOC(pc->masksize, int, pc->maskdim);
111 for (i=0; i<pc->maskdim; i++)
112 pc->masksize[i] = 1;
113 pc->masknpix = 1;
114 for (i=0; i<pc->maskdim; i++)
115 {
116 sprintf(str, "PCAXIS%1d ", i+1);
117 if (fitsread(head, str, &pc->masksize[i], H_INT,T_LONG) != RETURN_OK)
118 goto headerror;
119 pc->masknpix *= pc->masksize[i];
120 }
121
122 pc->npc = pc->masksize[pc->maskdim-1];
123
124 ncode = 0;
125 fitsread(head, "NCODE", &ncode, H_INT, T_LONG);
126 fitsread(head, "NCODEPAR", &nparam, H_INT, T_LONG);
127
128 /* Load the PC mask data */
129 key = read_key(tab, "PC_CONVMASK");
130 pc->maskcomp = key->ptr;
131
132 key = read_key(tab, "PC_MASK");
133 pc->omaskcomp = key->ptr;
134 pc->omaskdim = key->naxis;
135 pc->omasknpix = 1;
136 QMALLOC(pc->omasksize, int, pc->omaskdim);
137 for (i=0; i<pc->omaskdim; i++)
138 pc->omasknpix *= (pc->omasksize[i] = key->naxisn[i]);
139
140 key = read_key(tab, "PC_MX2");
141 pc->mx2 = key->ptr;
142
143 key = read_key(tab, "PC_MY2");
144 pc->my2 = key->ptr;
145
146 key = read_key(tab, "PC_MXY");
147 pc->mxy = key->ptr;
148
149 key = read_key(tab, "PC_FLUX");
150 pc->flux = key->ptr;
151
152 key = read_key(tab, "PC_BRATIO");
153 pc->bt = key->ptr;
154
155 if (ncode)
156 {
157 QMALLOC(pc->code, codestruct, 1);
158 code = pc->code;
159 QMALLOC(code->param, float *, nparam);
160 QMALLOC(code->parammod, int, nparam);
161 code->ncode = ncode;
162 code->nparam = nparam;
163 key = read_key(tab, "CODE_PC");
164 code->pc = (float *)key->ptr;
165 for (i=0; i<nparam; i++)
166 {
167 sprintf(str, "CODE_P%d", i+1);
168 key = read_key(tab, str);
169 code->param[i] = (float *)key->ptr;
170 sprintf(str, "CODE_M%d", i+1);
171 fitsread(head, str, &code->parammod[i], H_INT, T_LONG);
172 }
173 }
174
175 QMALLOC(pc->maskcurr, double, pc->masksize[0]*pc->masksize[1]*pc->npc);
176
177 /* But don't touch my arrays!! */
178 blank_keys(tab);
179
180 return pc;
181
182 headerror:
183 error(EXIT_FAILURE, "*Error*: Incorrect or obsolete PC data in ", filename);
184 return NULL;
185 }
186
187
188 /********************************** pc_fit **********************************/
189 /*
190 Fit the PC data to the current data.
191 */
pc_fit(psfstruct * psf,double * data,double * weight,int width,int height,int ix,int iy,double dx,double dy,int npc,float backrms)192 void pc_fit(psfstruct *psf, double *data, double *weight,
193 int width, int height,int ix, int iy,
194 double dx, double dy, int npc, float backrms)
195 {
196 pcstruct *pc;
197 checkstruct *check;
198 codestruct *code;
199 double *basis,*basis0, *cpix,*cpix0, *pcshift,*wpcshift,
200 *spix,*wspix, *w, *sumopc,*sumopct, *checkbuf,
201 *sol,*solt, *datat,
202 *mx2t, *my2t, *mxyt,
203 val,val2, xm2,ym2,xym,flux, temp,temp2, theta, pmx2,pmy2,
204 wnorm, ellip, norm, snorm;
205 float **param, *ppix, *ospix, *cpc,*cpc2, *fparam,
206 pixstep, fval, fvalmax, fscale, dparam;
207 int *parammod,
208 c,n,p, npix,npix2,nopix, ncoeff, nparam, nmax,nmax2, ncode;
209
210 pc = psf->pc;
211 /* Build the "local PCs", using the basis func. coeffs computed in psf_fit() */
212 if (npc > pc->npc)
213 npc = pc->npc;
214 npix = pc->masksize[0]*pc->masksize[1];
215 npix2 = width*height;
216 ncoeff = psf->poly->ncoeff;
217 pixstep = 1.0/psf->pixstep;
218 dx *= pixstep;
219 dy *= pixstep;
220
221 memset(pc->maskcurr, 0, npix*npc*sizeof(double));
222 basis0 = psf->poly->basis;
223 cpix0 = pc->maskcurr;
224 ppix = pc->maskcomp;
225
226 /* Sum each (PSF-dependent) component */
227 for (c=npc; c--; cpix0 += npix)
228 {
229 basis = basis0;
230 for (n = ncoeff; n--;)
231 {
232 cpix = cpix0;
233 val = *(basis++);
234 for (p=npix; p--;)
235 *(cpix++) += val*(double)*(ppix++);
236 }
237 }
238
239 /* Allocate memory for temporary buffers */
240 QMALLOC(pcshift, double, npix2*npc);
241 QMALLOC(wpcshift, double, npix2*npc);
242 QMALLOC(sol, double, npc);
243
244 /* Now shift and scale to the right position, and weight the PCs */
245 cpix = pc->maskcurr;
246 spix = pcshift;
247 wspix = wpcshift;
248 for (c=npc; c--; cpix += npix)
249 {
250 vignet_resample(cpix, pc->masksize[0], pc->masksize[1],
251 spix, width, height, -dx, -dy, pixstep);
252 w = weight;
253 for (p=npix2; p--;)
254 *(wspix++) = *(spix++)**(w++);
255 }
256
257 /* Compute the weight normalization */
258 wnorm = 0.0;
259 w = weight;
260 for (p=npix2; p--;)
261 {
262 val = *(w++);
263 wnorm += val*val;
264 }
265
266 /* Scalar product of data and (approximately orthogonal) basis functions */
267 wspix = wpcshift;
268 solt = sol;
269 snorm = 0.0;
270 for (c=npc; c--;)
271 {
272 datat = data;
273 val = 0.0;
274 for (p=npix2; p--;)
275 val += *(datat++)**(wspix++);
276 val2 = *(solt++) = val*npix2/wnorm;
277 snorm += val2*val2;
278 }
279
280
281 /* Normalize solution vector */
282 snorm = sqrt(snorm);
283 solt = sol;
284 for (c=npc; c--;)
285 *(solt++) /= snorm;
286
287 if ((code = pc->code))
288 {
289 ncode = code->ncode;
290 /*-- Codebook search */
291 cpc = code->pc;
292 fvalmax = -BIG;
293 nmax = 0;
294 for (n=ncode; n--;)
295 {
296 fval = 0.0;
297 solt = sol;
298 for (p=npc; p--;)
299 fval += *(solt++)**(cpc++);
300 if (fval>fvalmax)
301 {
302 fvalmax = fval;
303 nmax = n;
304 }
305 }
306 nmax = ncode - 1 - nmax;
307
308 /*-- Interpolation */
309 param = code->param;
310 parammod = code->parammod;
311 nparam = code->nparam;
312 QMALLOC(fparam, float, nparam);
313 for (p=0; p<nparam; p++)
314 {
315 dparam = 0.0;
316 if (parammod[p])
317 {
318 val2 = 0.0;
319 if ((nmax2 = nmax+parammod[p]) < ncode)
320 {
321 cpc = code->pc+npc*nmax;
322 cpc2 = code->pc+npc*nmax2;
323 solt = sol;
324 norm = 0.0;
325 for (c=npc; c--;)
326 {
327 val = *(cpc2++)-*cpc;
328 val2 += val*(*(solt++) - *(cpc++));
329 norm += val*val;
330 }
331 if (norm>0.0)
332 dparam = val2/norm*(param[p][nmax2]-param[p][nmax]);
333 else
334 val2 = 0.0;
335 }
336 /*------ If dot product negative of something went wrong, try other side */
337 if (val2<=0.0 && (nmax2 = nmax-parammod[p]) >= 0)
338 {
339 cpc = code->pc+npc*nmax;
340 cpc2 = code->pc+npc*nmax2;
341 solt = sol;
342 norm = val2 = 0.0;
343 for (c=npc; c--;)
344 {
345 val = *(cpc2++)-*cpc;
346 val2 += val*(*(solt++) - *(cpc++));
347 norm += val*val;
348 }
349 if (norm>0.0)
350 dparam = val2/norm*(param[p][nmax2]-param[p][nmax]);
351 }
352 fparam[p] = param[p][nmax] + dparam;
353 }
354 }
355
356 solt = sol;
357 cpc = code->pc+npc*nmax;
358 fscale = fvalmax*code->param[0][nmax]*snorm;
359 for (p=npc; p--;)
360 *(solt++) = fscale**(cpc++);
361
362 /*-- Copy the derived physical quantities to output parameters */
363 /*-- (subject to changes) */
364 obj2->flux_galfit = fscale;
365 obj2->gdposang = fparam[1];
366 if (obj2->gdposang>90.0)
367 obj2->gdposang -= 180.0;
368 else if (obj2->gdposang<-90.0)
369 obj2->gdposang += 180.0;
370 obj2->gdscale = fparam[2];
371 obj2->gdaspect = fparam[3];
372 ellip = (1.0 - obj2->gdaspect)/(1.0 + obj2->gdaspect);
373 obj2->gde1 = (float)(ellip*cos(2*obj2->gdposang*PI/180.0));
374 obj2->gde2 = (float)(ellip*sin(2*obj2->gdposang*PI/180.0));
375 /*---- Copy the best-fitting PCs to the VECTOR_PC output vector */
376 if (FLAG(obj2.vector_pc))
377 {
378 solt = sol;
379 ppix = obj2->vector_pc;
380 for (c=prefs.pc_vectorsize>npc?npc:prefs.pc_vectorsize; c--;)
381 *(ppix++) = *(solt++);
382 }
383
384 free(fparam);
385 }
386
387 xm2 = ym2 = xym = flux = 0.0;
388 solt = sol;
389 mx2t = pc->mx2;
390 my2t = pc->my2;
391 mxyt = pc->mxy;
392 for (c=npc; c--;)
393 {
394 val = *(solt++);
395 xm2 += val**(mx2t++);
396 ym2 += val**(my2t++);
397 xym += val**(mxyt++);
398 }
399
400 obj2->mx2_pc = xm2;
401 obj2->my2_pc = ym2;
402 obj2->mxy_pc = xym;
403
404 if (FLAG(obj2.a_pc))
405 {
406 /* Handle fully correlated x/y (which cause a singularity...) */
407 if ((temp2=xm2*ym2-xym*xym)<0.00694)
408 {
409 xm2 += 0.0833333;
410 ym2 += 0.0833333;
411 temp2 = xm2*ym2-xym*xym;
412 }
413
414 if ((fabs(temp=xm2-ym2)) > 0.0)
415 theta = atan2(2.0 * xym,temp) / 2.0;
416 else
417 theta = PI/4.0;
418
419 temp = sqrt(0.25*temp*temp+xym*xym);
420 pmy2 = pmx2 = 0.5*(xm2+ym2);
421 pmx2 += temp;
422 pmy2 -= temp;
423
424 obj2->a_pc = (float)sqrt(pmx2);
425 obj2->b_pc = (float)sqrt(pmy2);
426 obj2->theta_pc = (float)(theta*180.0/PI);
427 }
428
429 /* CHECK-Images */
430 if (prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS])
431 {
432 spix = pcshift;
433 solt = sol;
434 for (c=npc; c--; solt++)
435 {
436 ppix = checkmask;
437 for (p=npix2; p--;)
438 *(ppix++) = (PIXTYPE)*(spix++);
439 if ((check = prefs.check[CHECK_SUBPCPROTOS]))
440 addcheck(check, checkmask, width,height, ix,iy, -*solt);
441 if ((check = prefs.check[CHECK_PCPROTOS]))
442 addcheck(check, checkmask, width,height, ix,iy, *solt);
443 }
444 }
445 if ((check = prefs.check[CHECK_PCOPROTOS]))
446 {
447 /*- Reconstruct the unconvolved profile */
448 nopix = pc->omasksize[0]*pc->omasksize[1];
449 QCALLOC(sumopc, double, nopix);
450 solt = sol;
451 ospix = pc->omaskcomp;
452 for (c=npc; c--;)
453 {
454 val = *(solt++);
455 sumopct = sumopc;
456 for (p=nopix; p--;)
457 *(sumopct++) += val*(double)*(ospix++);
458 }
459 QMALLOC(checkbuf, double, npix2);
460 vignet_resample(sumopc, pc->omasksize[0], pc->omasksize[1],
461 checkbuf, width, height, -dx, -dy, pixstep);
462 ppix = checkmask;
463 spix = checkbuf;
464 for (p=npix2; p--;)
465 *(ppix++) = (PIXTYPE)*(spix++);
466 addcheck(check, checkmask, width,height, ix,iy, 1.0);
467 free(checkbuf);
468 free(sumopc);
469 }
470
471 /* Free memory */
472 free(pcshift);
473 free(wpcshift);
474 free(sol);
475
476 return;
477 }
478
479