CR2RE Pipeline Reference Manual 1.6.11
hdrl_cat_table.c
1/*
2 * This file is part of the HDRL
3 * Copyright (C) 2017 European Southern Observatory
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 */
19
20#include "hdrl_cat_table.h"
21
22#include "hdrl_cat_areals.h"
23#include "hdrl_cat_radii.h"
24#include "hdrl_cat_phopt.h"
25#include "hdrl_cat_seeing.h"
26#include "hdrl_cat_overlp.h"
27#include "hdrl_cat_moments.h"
28#include "hdrl_cat_extend.h"
29#include "hdrl_cat_background.h"
30
31
32#define INITROWS 2048 /* Allocation size for rows in the output table */
33
34/* Number assigned of columns in the table */
35#define COL_NUMBER 1
36#define COL_FLUXISO 2
37#define COL_X 3
38#define COL_XERR 4
39#define COL_Y 5
40#define COL_YERR 6
41#define COL_SIGMA 7
42#define COL_ELLIPT 8
43#define COL_PA 9
44#define COL_AREAL1 10
45#define COL_AREAL2 11
46#define COL_AREAL3 12
47#define COL_AREAL4 13
48#define COL_AREAL5 14
49#define COL_AREAL6 15
50#define COL_AREAL7 16
51#define COL_AREAL8 17
52#define COL_PEAKHEIGHT 18
53#define COL_PKHTERR 19
54#define COL_APFLUX1 20
55#define COL_APFLUX1ERR 21
56#define COL_APFLUX2 22
57#define COL_APFLUX2ERR 23
58#define COL_APFLUX3 24
59#define COL_APFLUX3ERR 25
60#define COL_APFLUX4 26
61#define COL_APFLUX4ERR 27
62#define COL_APFLUX5 28
63#define COL_APFLUX5ERR 29
64#define COL_APFLUX6 30
65#define COL_APFLUX6ERR 31
66#define COL_APFLUX7 32
67#define COL_APFLUX7ERR 33
68#define COL_APFLUX8 34
69#define COL_APFLUX8ERR 35
70#define COL_APFLUX9 36
71#define COL_APFLUX9ERR 37
72#define COL_APFLUX10 38
73#define COL_APFLUX10ERR 39
74#define COL_APFLUX11 40
75#define COL_APFLUX11ERR 41
76#define COL_APFLUX12 42
77#define COL_APFLUX12ERR 43
78#define COL_APFLUX13 44
79#define COL_APFLUX13ERR 45
80#define COL_PETRAD 46
81#define COL_KRONRAD 47
82#define COL_HALFRAD 48
83#define COL_PETFLUX 49
84#define COL_PETFLUXERR 50
85#define COL_KRONFLUX 51
86#define COL_KRONFLUXERR 52
87#define COL_HALFFLUX 53
88#define COL_HALFFLUXERR 54
89#define COL_ERRFLAG 55
90#define COL_SKYLEVEL 56
91#define COL_SKYSIGMA 57
92#define COL_AVCONF 58
93#define COL_RA 59
94#define COL_DEC 60
95#define COL_CLASS 61
96#define COL_STAT 62
97#define COL_FWHM 63
98
99/* Name assigned of columns in the table */
100static const char *ttype[NCOLS]={"Sequence_number","Isophotal_flux",
101 "X_coordinate","X_coordinate_err",
102 "Y_coordinate","Y_coordinate_err",
103 "Gaussian_sigma","Ellipticity","Position_angle",
104 "Areal_1_profile","Areal_2_profile","Areal_3_profile",
105 "Areal_4_profile","Areal_5_profile","Areal_6_profile",
106 "Areal_7_profile","Areal_8_profile",
107 "Peak_height","Peak_height_err",
108 "Aper_flux_1","Aper_flux_1_err",
109 "Aper_flux_2","Aper_flux_2_err",
110 "Aper_flux_3","Aper_flux_3_err",
111 "Aper_flux_4","Aper_flux_4_err",
112 "Aper_flux_5","Aper_flux_5_err",
113 "Aper_flux_6","Aper_flux_6_err",
114 "Aper_flux_7","Aper_flux_7_err",
115 "Aper_flux_8","Aper_flux_8_err",
116 "Aper_flux_9","Aper_flux_9_err",
117 "Aper_flux_10","Aper_flux_10_err",
118 "Aper_flux_11","Aper_flux_11_err",
119 "Aper_flux_12","Aper_flux_12_err",
120 "Aper_flux_13","Aper_flux_13_err",
121 "Petr_radius","Kron_radius","Half_radius",
122 "Petr_flux","Petr_flux_err",
123 "Kron_flux","Kron_flux_err","Half_flux","Half_flux_err",
124 "Error_bit_flag","Sky_level","Sky_rms",
125 "Av_conf",
126 "RA","DEC","Classification","Statistic",
127 "FWHM"};
128
129static const char *tunit[NCOLS]={"Number","ADU",
130 "Pixels","Pixels",
131 "Pixels","Pixels",
132 "Pixels","Number","Degrees",
133 "Pixels","Pixels","Pixels",
134 "Pixels","Pixels","Pixels",
135 "Pixels","Pixels",
136 "ADU","ADU",
137 "ADU","ADU",
138 "ADU","ADU",
139 "ADU","ADU",
140 "ADU","ADU",
141 "ADU","ADU",
142 "ADU","ADU",
143 "ADU","ADU",
144 "ADU","ADU",
145 "ADU","ADU",
146 "ADU","ADU",
147 "ADU","ADU",
148 "ADU","ADU",
149 "ADU","ADU",
150 "Pixels","Pixels","Pixels",
151 "ADU","ADU",
152 "ADU","ADU","ADU","ADU",
153 "Number","ADU","ADU","Number",
154 "Degrees","Degrees","Flag","N-sigma",
155 "Pixels"};
156
157static cpl_type g_tform[NCOLS]={CPL_TYPE_INT, CPL_TYPE_DOUBLE,
158 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
159 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
160 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
161 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
162 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
163 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
164 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
165 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
166 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
167 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
168 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
169 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
170 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
171 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
172 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
173 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
174 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
175 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
176 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
177 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
178 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
179 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
180 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
181 CPL_TYPE_DOUBLE,
182 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
183 CPL_TYPE_DOUBLE,
184 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE,
185 CPL_TYPE_DOUBLE, CPL_TYPE_DOUBLE};
186
187/* Apertures and areals */
188static double g_apertures[NRADS];
189
190static double g_rmults[] = {0.5,
191 1. / CPL_MATH_SQRT2,
192 1.,
193 CPL_MATH_SQRT2,
194 2.,
195 2. * CPL_MATH_SQRT2,
196 4., 5., 6., 7., 8., 10., 12.0};
197
198static cpl_size g_nrcore = 2;
199static cpl_size g_n2rcore = 4;
200
201static cpl_size g_areal_cols[NAREAL] = {COL_AREAL1, COL_AREAL2, COL_AREAL3,
202 COL_AREAL4, COL_AREAL5, COL_AREAL6,
203 COL_AREAL7, COL_AREAL8};
204
205
206/*** Prototypes ***/
207
208static cpl_error_code hdrl_do_seeing_gen(
209 ap_t *ap, const char *col_ellipt, const char *col_pkht,
210 const char *col_areals[NAREAL], cpl_size nobjects, cpl_table *tab);
211
212static cpl_error_code hdrl_tabinit_gen(
213 cpl_size ncols, cpl_type tform[], cpl_table **tab);
214
215
216/*----------------------------------------------------------------------------*/
224/*----------------------------------------------------------------------------*/
225
227
228/* ---------------------------------------------------------------------------*/
245/* ---------------------------------------------------------------------------*/
246cpl_error_code hdrl_tabinit(
247 ap_t *ap, cpl_size *xcol, cpl_size *ycol, hdrl_catalogue_options cattype,
248 cpl_table **tab, hdrl_casu_result *res)
249{
250 /* Define RA and Dec columns */
251 *xcol = COL_X;
252 *ycol = COL_Y;
253
254 /* Call the generic routine to open a new output table */
255 cpl_error_code e = hdrl_tabinit_gen(NCOLS, g_tform, tab);
256 if (e == CPL_ERROR_NONE){
257
258 if (cattype & HDRL_CATALOGUE_SEGMAP) {
259 res->segmentation_map = cpl_image_new(ap->lsiz, ap->csiz, CPL_TYPE_INT);
260 } else {
261 res->segmentation_map = NULL;
262 }
263
264 if (cattype & HDRL_CATALOGUE_BKG) {
265 res->background = cpl_image_new(ap->lsiz, ap->csiz, CPL_TYPE_DOUBLE);
266 } else {
267 res->background = NULL;
268 }
269 }
270
271 return e;
272}
273
274/* ---------------------------------------------------------------------------*/
289/* ---------------------------------------------------------------------------*/
290cpl_error_code hdrl_do_seeing(
291 ap_t *ap, cpl_size nobjects, cpl_table *tab)
292{
293 /* Sort out the areal profile column names */
294 const char *areal_colnames[NAREAL];
295 for (cpl_size i = 0; i < NAREAL; i++) {
296 areal_colnames[i] = (const char *)ttype[g_areal_cols[i] - 1];
297 }
298
299 /* Just call the generic seeing routine */
300 return hdrl_do_seeing_gen( ap,
301 ttype[COL_ELLIPT - 1],
302 ttype[COL_PEAKHEIGHT - 1],
303 areal_colnames, nobjects, tab);
304}
305
306/* ---------------------------------------------------------------------------*/
323/* ---------------------------------------------------------------------------*/
324cpl_error_code hdrl_process_results(
325 ap_t *ap, double gain, cpl_size *nobjects, cpl_table *tab, hdrl_casu_result *res)
326{
327 /* Do a basic moments analysis and work out the areal profiles*/
328 double momresults[8];
329 hdrl_moments(ap, momresults);
330 if (momresults[0] < 0) {
331 return CPL_ERROR_ILLEGAL_INPUT;
332 }
333
334 cpl_size iareal[NAREAL];
335 hdrl_areals(ap, iareal);
336
337 /* See if this object makes the cut in terms of its size. If not, then just return with good status */
338 if (iareal[0] < ap->ipnop || momresults[3] < ap->xintmin) {
339 return CPL_ERROR_NONE;
340 }
341
342 /* Work out the total flux */
343 double ttotal;
344 hdrl_extend( ap, momresults[3], momresults[1], momresults[2],
345 momresults[4], momresults[5], momresults[6], (double)(iareal[0]),
346 momresults[7], &ttotal);
347
348 /* Try and deblend the images if it is requested and justified */
349 cpl_size nbit;
350 double parmall[IMNUM][NPAR];
351 if (iareal[0] >= ap->mulpix && ap->icrowd) {
352 hdrl_overlp( ap,
353 parmall,
354 &nbit,
355 momresults[1],
356 momresults[2],
357 momresults[3],
358 iareal[0],
359 momresults[7]);
360 } else {
361 nbit = 1;
362 }
363
364 if (nbit == 1) {
365
366 parmall[0][0] = momresults[3];
367 parmall[0][1] = momresults[1];
368 parmall[0][2] = momresults[2];
369 parmall[0][3] = ap->thresh;
370
371 for (cpl_size i = 4; i < 8; i++) {
372 parmall[0][i] = momresults[i];
373 }
374
375 for (cpl_size i = 0; i < NAREAL; i++) {
376 parmall[0][i + 8] = (double)(iareal[i]);
377 }
378
379 } else {
380
381 cpl_size mbit = 0;
382 for (cpl_size i = 0; i < nbit; i++) {
383
384 if ( parmall[i][1] > 1. && parmall[i][1] < ap->lsiz
385 && parmall[i][2] > 1. && parmall[i][2] < ap->csiz) {
386
387 for (cpl_size j = 0; j < NPAR; j++) {
388 parmall[mbit][j] = parmall[i][j];
389 }
390
391 mbit++;
392 }
393 }
394 nbit = mbit;
395
396 if (nbit == 0) {
397 return CPL_ERROR_NONE;
398 }
399 }
400
401 /* Create a list of g_apertures */
402 double badpix[IMNUM];
403 double skyvar[IMNUM];
404 double avconf[IMNUM];
405 for (cpl_size i = 0; i < NRADS; i++) {
406 g_apertures[i] = g_rmults[i] * (ap->rcore);
407 skyvar[i] = CPL_MATH_PI * g_apertures[i] * g_apertures[i];
408 }
409
410 double rcore_area = CPL_MATH_PI * pow(ap->rcore, 2.);
411
412 /* Initialise the badpix accumulator */
413 for (cpl_size i = 0; i < nbit; i++) {
414 badpix[i] = 0.;
415 avconf[i] = 0.;
416 }
417
418 /* Get the core fluxes in all g_apertures */
419 double cflux[NRADS * IMNUM];
420 hdrl_phopt(ap, parmall, nbit, NRADS, g_apertures, cflux, badpix, g_nrcore, avconf);
421 for (cpl_size i = 0; i < nbit; i++) {
422 avconf[i] /= rcore_area;
423 }
424
425 /* Get half-light radius for all images */
426 double half_flux[IMNUM];
427 double half_rad[ IMNUM];
428 for (cpl_size k = 0; k < nbit; k++) {
429
430 half_flux[k] = 0.5 * (CPL_MAX(parmall[k][0], cflux[k * NRADS + g_n2rcore]));
431
432 half_rad[k] = hdrl_halflight( g_apertures,
433 cflux+k*NRADS,
434 half_flux[k],
435 parmall[k][7],
436 NRADS);
437 }
438
439 /* Get Kron radius for all images and get the flux */
440 double areal;
441 double kron_flux[IMNUM];
442 double kron_rad[ IMNUM];
443 for (cpl_size k = 0; k < nbit; k++) {
444 areal = parmall[k][8];
445 kron_rad[k] = hdrl_kronrad(areal, g_apertures, cflux + k * NRADS, NRADS);
446 }
447 hdrl_flux(ap, parmall, nbit, kron_rad, kron_flux, NRADS, g_apertures, cflux);
448
449 /* Get Petrosian radius for all images and get the flux */
450 double petr_flux[IMNUM];
451 double petr_rad[ IMNUM];
452 for (cpl_size k = 0; k < nbit; k++) {
453 areal = parmall[k][8];
454 petr_rad[k] = hdrl_petrad(areal, g_apertures, cflux + k * NRADS, NRADS);
455 }
456 hdrl_flux(ap, parmall, nbit, petr_rad, petr_flux, NRADS, g_apertures, cflux);
457
458 /* Massage the results and write them to the fits table */
459 double sigsq = pow(ap->sigma, 2.);
460 double radeg = 180. / CPL_MATH_PI;
461
462 for (cpl_size k = 0; k < nbit; k++) {
463
464 double sxx = parmall[k][4];
465 double sxy = parmall[k][5];
466 double syy = parmall[k][6];
467
468 double srr = CPL_MAX(0.5, sxx + syy);
469
470 if (sxy > 0.) {
471 sxy = CPL_MAX( 1.e-4, CPL_MIN(sxy, sqrt(sxx * syy)));
472 } else {
473 sxy = CPL_MIN(-1.e-4, CPL_MAX(sxy, -sqrt(sxx * syy)));
474 }
475
476 double ecc = sqrt((syy - sxx) * (syy - sxx) + 4. * sxy * sxy) / srr;
477
478 double temp = CPL_MAX((1. - ecc) / (1. + ecc), 0.);
479 double ell = CPL_MIN(0.99, CPL_MAX(0., 1. - sqrt(temp)));
480
481 double xx = 0.5 * (1. + ecc) * srr - sxx;
482 double yy;
483
484 double theta;
485 if (xx == 0.) {
486 theta = 0.;
487 } else {
488 theta = 90. - radeg * atan(sxy / xx);
489 }
490 double theta_ra = theta/radeg;
491
492 double cc = (1. + ecc) * pow(cos(theta_ra), 2.) + (1. - ecc) * pow(sin(theta_ra), 2.);
493 double dd = (1. + ecc) * pow(sin(theta_ra), 2.) + (1. - ecc) * pow(cos(theta_ra), 2.);
494
495 /* Create a list of values */
496 cpl_size nrows = cpl_table_get_nrow(tab);
497
498 (*nobjects)++;
499 if (*nobjects > nrows) {
500 cpl_table_set_size(tab, nrows + INITROWS);
501 }
502 cpl_size nr = *nobjects - 1;
503
504 double iso_flux = parmall[k][0];
505
506 double apflux1 = cflux[k * NRADS + 0];
507 double apflux2 = cflux[k * NRADS + 1];
508 double apflux3 = cflux[k * NRADS + 2];
509 double apflux4 = cflux[k * NRADS + 3];
510 double apflux5 = cflux[k * NRADS + 4];
511 double apflux6 = cflux[k * NRADS + 5];
512 double apflux7 = cflux[k * NRADS + 6];
513 double apflux8 = cflux[k * NRADS + 7];
514 double apflux9 = cflux[k * NRADS + 8];
515 double apflux10 = cflux[k * NRADS + 9];
516 double apflux11 = cflux[k * NRADS + 10];
517 double apflux12 = cflux[k * NRADS + 11];
518 double apflux13 = cflux[k * NRADS + 12];
519
520 double peak = parmall[k][7];
521
522 xx = parmall[k][1];
523 yy = parmall[k][2];
524
525 double xxe = sqrt( (2. * sigsq / (CPL_MATH_PI * peak * peak))
526 + cc / (2. * CPL_MATH_PI * gain * peak)
527 + 0.0001);
528
529 double yye = sqrt( (2. * sigsq / (CPL_MATH_PI * peak * peak))
530 + dd / (2. * CPL_MATH_PI * gain * peak)
531 + 0.0001);
532
533 double sigma = sqrt(srr);
534 double fwhm = sqrt(sigma * sigma / 2.) * CPL_MATH_FWHM_SIG;
535 /* heuristic correction of moment based fwhm obtained via simulated 2d gaussians,
536 * with gaussians 4.3 corrects slighltly better but it is not very
537 * significant and 4.0 is the same factor sextractor uses */
538 fwhm -= 1. / (4. * fwhm);
539
540 double areal1 = parmall[k][8];
541 double areal2 = parmall[k][9];
542 double areal3 = parmall[k][10];
543 double areal4 = parmall[k][11];
544 double areal5 = parmall[k][12];
545 double areal6 = parmall[k][13];
546 double areal7 = parmall[k][14];
547
548 double areal8;
549 if (nbit > 1 && k == 0) {
550 areal8 = 0.;
551 } else {
552 areal8 = parmall[k][15];
553 }
554
555 double skylev;
556 double skyrms;
557 hdrl_backest(ap, xx, yy, &skylev, &skyrms);
558
559 double kron_fluxe = sqrt(kron_flux[k] / gain + (sigsq + skyrms * skyrms) * CPL_MATH_PI * pow(kron_rad[k], 2.));
560 double petr_fluxe = sqrt(petr_flux[k] / gain + (sigsq + skyrms * skyrms) * CPL_MATH_PI * pow(petr_rad[k], 2.));
561 double half_fluxe = sqrt(CPL_MAX(0., half_flux[k]) / gain + (sigsq + skyrms * skyrms) * CPL_MATH_PI * pow(half_rad[k], 2.));
562
563 double apflux1e = sqrt(CPL_MAX(0., apflux1 / gain) + skyvar[0 ] * (sigsq + skyrms * skyrms));
564 double apflux2e = sqrt(CPL_MAX(0., apflux2 / gain) + skyvar[1 ] * (sigsq + skyrms * skyrms));
565 double apflux3e = sqrt(CPL_MAX(0., apflux3 / gain) + skyvar[2 ] * (sigsq + skyrms * skyrms));
566 double apflux4e = sqrt(CPL_MAX(0., apflux4 / gain) + skyvar[3 ] * (sigsq + skyrms * skyrms));
567 double apflux5e = sqrt(CPL_MAX(0., apflux5 / gain) + skyvar[4 ] * (sigsq + skyrms * skyrms));
568 double apflux6e = sqrt(CPL_MAX(0., apflux6 / gain) + skyvar[5 ] * (sigsq + skyrms * skyrms));
569 double apflux7e = sqrt(CPL_MAX(0., apflux7 / gain) + skyvar[6 ] * (sigsq + skyrms * skyrms));
570 double apflux8e = sqrt(CPL_MAX(0., apflux8 / gain) + skyvar[7 ] * (sigsq + skyrms * skyrms));
571 double apflux9e = sqrt(CPL_MAX(0., apflux9 / gain) + skyvar[8 ] * (sigsq + skyrms * skyrms));
572 double apflux10e = sqrt(CPL_MAX(0., apflux10 / gain) + skyvar[9 ] * (sigsq + skyrms * skyrms));
573 double apflux11e = sqrt(CPL_MAX(0., apflux11 / gain) + skyvar[10] * (sigsq + skyrms * skyrms));
574 double apflux12e = sqrt(CPL_MAX(0., apflux12 / gain) + skyvar[11] * (sigsq + skyrms * skyrms));
575 double apflux13e = sqrt(CPL_MAX(0., apflux13 / gain) + skyvar[12] * (sigsq + skyrms * skyrms));
576
577 double peake = sqrt(peak / gain + sigsq + skyrms * skyrms) ;
578
579 /* Store away the results for this object */
580 cpl_table_set_int( tab, ttype[COL_NUMBER - 1], nr, *nobjects);
581 cpl_table_set_double(tab, ttype[COL_FLUXISO - 1], nr, iso_flux);
582 cpl_table_set_double(tab, ttype[COL_X - 1], nr, xx);
583 cpl_table_set_double(tab, ttype[COL_XERR - 1], nr, xxe);
584 cpl_table_set_double(tab, ttype[COL_Y - 1], nr, yy);
585 cpl_table_set_double(tab, ttype[COL_YERR - 1], nr, yye);
586 cpl_table_set_double(tab, ttype[COL_SIGMA - 1], nr, sigma);
587 cpl_table_set_double(tab, ttype[COL_ELLIPT - 1], nr, ell);
588 cpl_table_set_double(tab, ttype[COL_PA - 1], nr, theta);
589 cpl_table_set_double(tab, ttype[COL_AREAL1 - 1], nr, areal1);
590 cpl_table_set_double(tab, ttype[COL_AREAL2 - 1], nr, areal2);
591 cpl_table_set_double(tab, ttype[COL_AREAL3 - 1], nr, areal3);
592 cpl_table_set_double(tab, ttype[COL_AREAL4 - 1], nr, areal4);
593 cpl_table_set_double(tab, ttype[COL_AREAL5 - 1], nr, areal5);
594 cpl_table_set_double(tab, ttype[COL_AREAL6 - 1], nr, areal6);
595 cpl_table_set_double(tab, ttype[COL_AREAL7 - 1], nr, areal7);
596 cpl_table_set_double(tab, ttype[COL_AREAL8 - 1], nr, areal8);
597 cpl_table_set_double(tab, ttype[COL_PEAKHEIGHT - 1], nr, peak);
598 cpl_table_set_double(tab, ttype[COL_PKHTERR - 1], nr, peake);
599 cpl_table_set_double(tab, ttype[COL_APFLUX1 - 1], nr, apflux1);
600 cpl_table_set_double(tab, ttype[COL_APFLUX1ERR - 1], nr, apflux1e);
601 cpl_table_set_double(tab, ttype[COL_APFLUX2 - 1], nr, apflux2);
602 cpl_table_set_double(tab, ttype[COL_APFLUX2ERR - 1], nr, apflux2e);
603 cpl_table_set_double(tab, ttype[COL_APFLUX3 - 1], nr, apflux3);
604 cpl_table_set_double(tab, ttype[COL_APFLUX3ERR - 1], nr, apflux3e);
605 cpl_table_set_double(tab, ttype[COL_APFLUX4 - 1], nr, apflux4);
606 cpl_table_set_double(tab, ttype[COL_APFLUX4ERR - 1], nr, apflux4e);
607 cpl_table_set_double(tab, ttype[COL_APFLUX5 - 1], nr, apflux5);
608 cpl_table_set_double(tab, ttype[COL_APFLUX5ERR - 1], nr, apflux5e);
609 cpl_table_set_double(tab, ttype[COL_APFLUX6 - 1], nr, apflux6);
610 cpl_table_set_double(tab, ttype[COL_APFLUX6ERR - 1], nr, apflux6e);
611 cpl_table_set_double(tab, ttype[COL_APFLUX7 - 1], nr, apflux7);
612 cpl_table_set_double(tab, ttype[COL_APFLUX7ERR - 1], nr, apflux7e);
613 cpl_table_set_double(tab, ttype[COL_APFLUX8 - 1], nr, apflux8);
614 cpl_table_set_double(tab, ttype[COL_APFLUX8ERR - 1], nr, apflux8e);
615 cpl_table_set_double(tab, ttype[COL_APFLUX9 - 1], nr, apflux9);
616 cpl_table_set_double(tab, ttype[COL_APFLUX9ERR - 1], nr, apflux9e);
617 cpl_table_set_double(tab, ttype[COL_APFLUX10 - 1], nr, apflux10);
618 cpl_table_set_double(tab, ttype[COL_APFLUX10ERR - 1], nr, apflux10e);
619 cpl_table_set_double(tab, ttype[COL_APFLUX11 - 1], nr, apflux11);
620 cpl_table_set_double(tab, ttype[COL_APFLUX11ERR - 1], nr, apflux11e);
621 cpl_table_set_double(tab, ttype[COL_APFLUX12 - 1], nr, apflux12);
622 cpl_table_set_double(tab, ttype[COL_APFLUX12ERR - 1], nr, apflux12e);
623 cpl_table_set_double(tab, ttype[COL_APFLUX13 - 1], nr, apflux13);
624 cpl_table_set_double(tab, ttype[COL_APFLUX13ERR - 1], nr, apflux13e);
625 cpl_table_set_double(tab, ttype[COL_PETRAD - 1], nr, 0.5 * petr_rad[k]);
626 cpl_table_set_double(tab, ttype[COL_KRONRAD - 1], nr, 0.5 * kron_rad[k]);
627 cpl_table_set_double(tab, ttype[COL_HALFRAD - 1], nr, half_rad[k]);
628 cpl_table_set_double(tab, ttype[COL_PETFLUX - 1], nr, petr_flux[k]);
629 cpl_table_set_double(tab, ttype[COL_PETFLUXERR - 1], nr, petr_fluxe);
630 cpl_table_set_double(tab, ttype[COL_KRONFLUX - 1], nr, kron_flux[k]);
631 cpl_table_set_double(tab, ttype[COL_KRONFLUXERR - 1], nr, kron_fluxe);
632 cpl_table_set_double(tab, ttype[COL_HALFFLUX - 1], nr, half_flux[k]);
633 cpl_table_set_double(tab, ttype[COL_HALFFLUXERR - 1], nr, half_fluxe);
634 cpl_table_set_double(tab, ttype[COL_ERRFLAG - 1], nr, badpix[k]);
635 cpl_table_set_double(tab, ttype[COL_SKYLEVEL - 1], nr, skylev);
636 cpl_table_set_double(tab, ttype[COL_SKYSIGMA - 1], nr, skyrms);
637 cpl_table_set_double(tab, ttype[COL_AVCONF - 1], nr, avconf[k]);
638 cpl_table_set_double(tab, ttype[COL_FWHM - 1], nr, fwhm);
639
640 /* Store away some dummy values to avoid problems later on */
641 double zero = 0.;
642 cpl_table_set_double(tab, ttype[COL_RA - 1], nr, zero);
643 cpl_table_set_double(tab, ttype[COL_DEC - 1], nr, zero);
644 cpl_table_set_double(tab, ttype[COL_CLASS - 1], nr, 100.);
645 cpl_table_set_double(tab, ttype[COL_STAT - 1], nr, zero);
646 }
647
648 /* Now that everything is okay - fill in the segmenataion map */
649 cpl_msg_info(cpl_func, "Num objects found in catalogue: %lld", *nobjects);
650 if (res->segmentation_map) {
651 for (cpl_size index = 0; index < ap->npl_pix; index++) {
652 cpl_image_set( res->segmentation_map,
653 ap->plarray[index].x,
654 ap->plarray[index].y,
655 *nobjects);
656 }
657 }
658
659 return CPL_ERROR_NONE;
660}
661
663
664/* ---------------------------------------------------------------------------*/
679/* ---------------------------------------------------------------------------*/
680static cpl_error_code hdrl_tabinit_gen(
681 cpl_size ncols, cpl_type tform[], cpl_table **tab)
682{
683 /* First, create the table with a default number of rows. */
684 if ((*tab = cpl_table_new(0)) == NULL) {
685
686 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
687 "hdrl_cat_tabinit_gen - Unable to open cpl table!");
688 return CPL_ERROR_ILLEGAL_INPUT;
689
690 } else {
691
692 /* Now define all of the columns */
693 for (cpl_size i = 0; i < ncols; i++) {
694
695 cpl_table_new_column( *tab, ttype[i], tform[i]);
696 cpl_table_set_column_unit(*tab, ttype[i], tunit[i]);
697 }
698
699 return CPL_ERROR_NONE;
700 }
701}
702
703/* ---------------------------------------------------------------------------*/
719/* ---------------------------------------------------------------------------*/
720static cpl_error_code hdrl_do_seeing_gen(
721 ap_t *ap, const char *col_ellipt, const char *col_pkht,
722 const char *col_areals[NAREAL], cpl_size nobjects, cpl_table *tab)
723{
724 /* Get some space and read the relevant columns */
725 double fwhm;
726 if (nobjects >= 3) {
727
728 double *areal[NAREAL];
729 for (cpl_size i = 0; i < NAREAL; i++) {
730 areal[i] = cpl_table_get_data_double(tab, col_areals[i]);
731 }
732
733 double *ellipt = cpl_table_get_data_double(tab, col_ellipt);
734 double *pkht = cpl_table_get_data_double(tab, col_pkht);
735
736 /* Do the seeing calculation */
737 double *work = cpl_malloc(nobjects * sizeof(double));
738 hdrl_seeing(ap, nobjects, ellipt, pkht, areal, work, &fwhm);
739 cpl_free(work);
740
741 } else {
742 fwhm = 0.;
743 }
744
745 ap->fwhm = fwhm;
746
747 return CPL_ERROR_NONE;
748}
hdrl_areals
void hdrl_areals(ap_t *ap, cpl_size iareal[NAREAL])
Work out the areal profiles for an object.
Definition hdrl_cat_areals.c:52
hdrl_backest
void hdrl_backest(ap_t *ap, double x, double y, double *skylev, double *skyrms)
Work out estimated sky for a pixel position.
Definition hdrl_cat_background.c:315
hdrl_process_results
cpl_error_code hdrl_process_results(ap_t *ap, double gain, cpl_size *nobjects, cpl_table *tab, hdrl_casu_result *res)
Process the results for each object and store them in the table.
Definition hdrl_cat_table.c:324
hdrl_do_seeing
cpl_error_code hdrl_do_seeing(ap_t *ap, cpl_size nobjects, cpl_table *tab)
Do seeing estimate.
Definition hdrl_cat_table.c:290
hdrl_tabinit
cpl_error_code hdrl_tabinit(ap_t *ap, cpl_size *xcol, cpl_size *ycol, hdrl_catalogue_options cattype, cpl_table **tab, hdrl_casu_result *res)
Initialize catalogues.
Definition hdrl_cat_table.c:246
hdrl_extend
cpl_error_code hdrl_extend(ap_t *ap, double xniso, double xbar, double ybar, double sxx, double sxy, double syy, double areal0, double tmax, double *ttotal)
Do aperture integration.
Definition hdrl_cat_extend.c:67
hdrl_moments
void hdrl_moments(ap_t *ap, double results[])
Do moments analysis on an object in a Plessey array.
Definition hdrl_cat_moments.c:47
hdrl_overlp
cpl_error_code hdrl_overlp(ap_t *ap, double parm[IMNUM][NPAR], cpl_size *nbit, double xbar, double ybar, double total, cpl_size npix, double tmax)
Deblend overlapping images.
Definition hdrl_cat_overlp.c:88
hdrl_phopt
cpl_error_code hdrl_phopt(ap_t *ap, double parm[IMNUM][NPAR], cpl_size nbit, cpl_size naper, double apertures[], double cflux[], double badpix[], cpl_size nrcore, double avconf[])
Does multiple profile fitting to determine intensities.
Definition hdrl_cat_phopt.c:59
hdrl_flux
void hdrl_flux(ap_t *ap, double parm[IMNUM][NPAR], cpl_size nbit, double apers[], double fluxes[], cpl_size nr, double rcores[], double rfluxes[])
Work out the fluxes for special radii.
Definition hdrl_cat_radii.c:225
hdrl_petrad
double hdrl_petrad(double areal0, double rcores[], double cflux[], cpl_size naper)
Work out the Petrosian radius for an object.
Definition hdrl_cat_radii.c:169
hdrl_kronrad
double hdrl_kronrad(double areal0, double rcores[], double cflux[], cpl_size naper)
Work out the Kron radius for an object.
Definition hdrl_cat_radii.c:132
hdrl_halflight
double hdrl_halflight(double rcores[], double cflux[], double halflight, double peak, cpl_size naper)
Work out the half-light radius for an object.
Definition hdrl_cat_radii.c:54
hdrl_seeing
cpl_error_code hdrl_seeing(ap_t *ap, cpl_size nrows, double *ellipt, double *pkht, double **areal, double *work, double *fwhm)
Work out the median seeing.
Definition hdrl_cat_seeing.c:57
Generated on Tue Jun 16 2026 00:11:35 for CR2RE Pipeline Reference Manual by doxygen 1.13.2