"Fossies" - the Fresh Open Source Software Archive

Member "gnuastro-0.9/bin/mkcatalog/columns.c" (11 Apr 2019, 73477 Bytes) of package /linux/privat/gnuastro-0.9.tar.lz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "columns.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 0.8_vs_0.9.

    1 /*********************************************************************
    2 MakeCatalog - Make a catalog from an input and labeled image.
    3 MakeCatalog is part of GNU Astronomy Utilities (Gnuastro) package.
    4 
    5 Original author:
    6      Mohammad Akhlaghi <mohammad@akhlaghi.org>
    7 Contributing author(s):
    8 Copyright (C) 2016-2019, Free Software Foundation, Inc.
    9 
   10 Gnuastro is free software: you can redistribute it and/or modify it
   11 under the terms of the GNU General Public License as published by the
   12 Free Software Foundation, either version 3 of the License, or (at your
   13 option) any later version.
   14 
   15 Gnuastro is distributed in the hope that it will be useful, but
   16 WITHOUT ANY WARRANTY; without even the implied warranty of
   17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   18 General Public License for more details.
   19 
   20 You should have received a copy of the GNU General Public License
   21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
   22 **********************************************************************/
   23 #include <config.h>
   24 
   25 #include <math.h>
   26 #include <stdio.h>
   27 #include <errno.h>
   28 #include <error.h>
   29 #include <stdlib.h>
   30 #include <string.h>
   31 #include <pthread.h>
   32 
   33 #include <gnuastro/pointer.h>
   34 
   35 #include <gnuastro-internal/checkset.h>
   36 
   37 #include "main.h"
   38 #include "mkcatalog.h"
   39 
   40 #include "ui.h"
   41 #include "columns.h"
   42 
   43 
   44 
   45 
   46 /******************************************************************/
   47 /*******************     Intermediate arrays     ******************/
   48 /******************************************************************/
   49 /* Allocate RA-DEC internal arrays. These arrays are defined to keep all
   50    the positions in one place and do the RA-DEC conversion once in the
   51    end. They are all allocated together, but we don't know if RA is
   52    requested first or Dec or if they are requested multiple times. So
   53    before the allocation, we'll check the first one.
   54 
   55    The space that is allocated in `columns_define_alloc' is for the final
   56    values that are written in the output file. */
   57 static void
   58 columns_alloc_radec(struct mkcatalogparams *p)
   59 {
   60   size_t i;
   61 
   62   /* For objects. */
   63   if(p->wcs_vo==NULL)
   64     for(i=0;i<p->objects->ndim;++i)
   65       gal_list_data_add_alloc(&p->wcs_vo, NULL, GAL_TYPE_FLOAT64, 1,
   66                               &p->numobjects, NULL, 0, p->cp.minmapsize,
   67                               NULL, NULL, NULL);
   68 
   69   /* For clumps */
   70   if(p->clumps && p->wcs_vc==NULL)
   71     for(i=0;i<p->objects->ndim;++i)
   72       gal_list_data_add_alloc(&p->wcs_vc, NULL, GAL_TYPE_FLOAT64, 1,
   73                               &p->numclumps, NULL, 0, p->cp.minmapsize,
   74                               NULL, NULL, NULL);
   75 }
   76 
   77 
   78 
   79 
   80 
   81 /* Similar to `columns_alloc_radec'. */
   82 static void
   83 columns_alloc_georadec(struct mkcatalogparams *p)
   84 {
   85   size_t i;
   86 
   87   /* For objects. */
   88   if(p->wcs_go==NULL)
   89     for(i=0;i<p->objects->ndim;++i)
   90       gal_list_data_add_alloc(&p->wcs_go, NULL, GAL_TYPE_FLOAT64, 1,
   91                               &p->numobjects, NULL, 0, p->cp.minmapsize,
   92                               NULL, NULL, NULL);
   93 
   94   /* For clumps */
   95   if(p->clumps && p->wcs_gc==NULL)
   96     for(i=0;i<p->objects->ndim;++i)
   97       gal_list_data_add_alloc(&p->wcs_gc, NULL, GAL_TYPE_FLOAT64, 1,
   98                               &p->numclumps, NULL, 0, p->cp.minmapsize,
   99                               NULL, NULL, NULL);
  100 }
  101 
  102 
  103 
  104 
  105 
  106 /* Similar to `columns_alloc_radec'. */
  107 static void
  108 columns_alloc_clumpsradec(struct mkcatalogparams *p)
  109 {
  110   size_t i;
  111 
  112   if(p->wcs_vcc==NULL)
  113     for(i=0;i<p->objects->ndim;++i)
  114       gal_list_data_add_alloc(&p->wcs_vcc, NULL, GAL_TYPE_FLOAT64, 1,
  115                               &p->numobjects, NULL, 0, p->cp.minmapsize,
  116                               NULL, NULL, NULL);
  117 }
  118 
  119 
  120 
  121 
  122 
  123 /* Similar to `columns_alloc_radec'. */
  124 static void
  125 columns_alloc_clumpsgeoradec(struct mkcatalogparams *p)
  126 {
  127   size_t i;
  128 
  129   if(p->wcs_gcc==NULL)
  130     for(i=0;i<p->objects->ndim;++i)
  131       gal_list_data_add_alloc(&p->wcs_gcc, NULL, GAL_TYPE_FLOAT64, 1,
  132                               &p->numobjects, NULL, 0, p->cp.minmapsize,
  133                               NULL, NULL, NULL);
  134 }
  135 
  136 
  137 
  138 
  139 
  140 /* Set pointers to fascilitate filling in the values. */
  141 #define SET_WCS_PREPARE(ARR, LIST, ARRNAME) {                           \
  142     d=0;                                                                \
  143     errno=0;                                                            \
  144     (ARR)=malloc(p->objects->ndim * sizeof (ARR) );                     \
  145     if( (ARR)==NULL )                                                   \
  146       error(EXIT_FAILURE, 0, "%s: %zu bytes for %s", __func__,          \
  147             p->objects->ndim * sizeof (ARR), ARRNAME);                  \
  148     for(tmp=(LIST);tmp!=NULL;tmp=tmp->next) (ARR)[d++]=tmp->array;      \
  149   }
  150 
  151 static void
  152 columns_set_wcs_pointers(struct mkcatalogparams *p, double ***vo,
  153                          double ***vc, double ***go, double ***gc,
  154                          double ***vcc, double ***gcc)
  155 {
  156   size_t d;
  157   gal_data_t *tmp;
  158 
  159   if(p->wcs_vo)  SET_WCS_PREPARE(*vo,  p->wcs_vo,  "vo" );
  160   if(p->wcs_vc)  SET_WCS_PREPARE(*vc,  p->wcs_vc,  "vc" );
  161   if(p->wcs_go)  SET_WCS_PREPARE(*go,  p->wcs_go,  "go" );
  162   if(p->wcs_gc)  SET_WCS_PREPARE(*gc,  p->wcs_gc,  "gc" );
  163   if(p->wcs_vcc) SET_WCS_PREPARE(*vcc, p->wcs_vcc, "vcc");
  164   if(p->wcs_gcc) SET_WCS_PREPARE(*gcc, p->wcs_gcc, "gcc");
  165 }
  166 
  167 
  168 
  169 
  170 
  171 
  172 
  173 
  174 
  175 
  176 
  177 
  178 
  179 
  180 
  181 
  182 
  183 
  184 
  185 
  186 /******************************************************************/
  187 /**********       Column definition/allocation      ***************/
  188 /******************************************************************/
  189 static void
  190 columns_wcs_preparation(struct mkcatalogparams *p)
  191 {
  192   size_t i;
  193   gal_list_i32_t *colcode;
  194   int continue_wcs_check=1;
  195 
  196   /* Make sure a WCS structure is present if we need it. */
  197   for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
  198     {
  199       if(continue_wcs_check)
  200         {
  201           switch(colcode->v)
  202             {
  203             /* High-level. */
  204             case UI_KEY_RA:
  205             case UI_KEY_DEC:
  206 
  207             /* Low-level. */
  208             case UI_KEY_W1:
  209             case UI_KEY_W2:
  210             case UI_KEY_GEOW1:
  211             case UI_KEY_GEOW2:
  212             case UI_KEY_CLUMPSW1:
  213             case UI_KEY_CLUMPSW2:
  214             case UI_KEY_CLUMPSGEOW1:
  215             case UI_KEY_CLUMPSGEOW2:
  216               if(p->objects->wcs)
  217                 continue_wcs_check=0;
  218               else
  219                 error(EXIT_FAILURE, 0, "%s (hdu: %s): no WCS meta-data "
  220                       "found by WCSLIB. Atleast one of the requested columns "
  221                       "requires world coordinate system meta-data",
  222                       p->objectsfile, p->cp.hdu);
  223               break;
  224             }
  225         }
  226       else
  227         break;
  228     }
  229 
  230   /* Convert the high-level WCS columns to low-level ones. */
  231   for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
  232     switch(colcode->v)
  233       {
  234       case UI_KEY_RA:
  235       case UI_KEY_DEC:
  236         /* Check all the CTYPES. */
  237         for(i=0;i<p->objects->ndim;++i)
  238           if( !strcmp(p->ctype[i], colcode->v==UI_KEY_RA ? "RA" : "DEC") )
  239             {
  240               colcode->v = i==0 ? UI_KEY_W1 : UI_KEY_W2;
  241               break;
  242             }
  243 
  244         /* Make sure it actually existed. */
  245         if(i==p->objects->ndim)
  246           error(EXIT_FAILURE, 0, "%s (hdu: %s): %s not present in any of "
  247                 "the WCS axis types (CTYPE)", p->objectsfile, p->cp.hdu,
  248                 colcode->v==UI_KEY_RA ? "RA" : "DEC");
  249         break;
  250       }
  251 }
  252 
  253 
  254 
  255 
  256 
  257 /* Set the necessary parameters for each output column and allocate the
  258    space necessary to keep the values. */
  259 void
  260 columns_define_alloc(struct mkcatalogparams *p)
  261 {
  262   gal_list_i32_t *colcode;
  263   gal_list_str_t *strtmp, *noclumpimg=NULL;
  264   int disp_fmt=0, disp_width=0, disp_precision=0;
  265   char *name=NULL, *unit=NULL, *ocomment=NULL, *ccomment=NULL;
  266   uint8_t otype=GAL_TYPE_INVALID, ctype=GAL_TYPE_INVALID, *oiflag, *ciflag;
  267 
  268   /* If there is any columns that need WCS, the input image needs to have a
  269      WCS in its headers. So before anything, we need to check if a WCS is
  270      present or not. This can't be done after the initial setting of column
  271      properties because the WCS-related columns use information that is
  272      based on it (for units and names). */
  273   columns_wcs_preparation(p);
  274 
  275   /* Allocate the array for which intermediate parameters are
  276      necessary. The basic issue is that higher-level calculations require a
  277      smaller domain of raw measurements. So to avoid having to calculate
  278      something multiple times, each parameter will flag the intermediate
  279      parameters it requires in these arrays. */
  280   oiflag = p->oiflag = gal_pointer_allocate(GAL_TYPE_UINT8, OCOL_NUMCOLS, 1,
  281                                             __func__, "oiflag");
  282   ciflag = p->ciflag = gal_pointer_allocate(GAL_TYPE_UINT8, CCOL_NUMCOLS, 1,
  283                                             __func__, "ciflag");
  284 
  285   /* Allocate the columns. */
  286   for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
  287     {
  288       /* Set the column-specific parameters, please follow the same order
  289          as `args.h'. IMPORTANT: we want the names to be the same as the
  290          option names. Note that zero `disp_' variables will be
  291          automatically determined.*/
  292       switch(colcode->v)
  293         {
  294         case UI_KEY_OBJID:
  295           name           = "OBJ_ID";
  296           unit           = "counter";
  297           ocomment       = "Object identifier.";
  298           ccomment       = NULL;
  299           otype          = GAL_TYPE_INT32;
  300           ctype          = GAL_TYPE_INVALID;
  301           disp_fmt       = 0;
  302           disp_width     = 6;
  303           disp_precision = 0;
  304           /* Is an internal parameter. */
  305           break;
  306 
  307         case UI_KEY_HOSTOBJID:
  308           name           = "HOST_OBJ_ID";
  309           unit           = "counter";
  310           ocomment       = NULL;
  311           ccomment       = "Object identifier hosting this clump.";
  312           otype          = GAL_TYPE_INVALID;
  313           ctype          = GAL_TYPE_INT32;
  314           disp_fmt       = 0;
  315           disp_width     = 6;
  316           disp_precision = 0;
  317           /* Is an internal parameter. */
  318           break;
  319 
  320         case UI_KEY_IDINHOSTOBJ:
  321           name           = "ID_IN_HOST_OBJ";
  322           unit           = "counter";
  323           ocomment       = NULL;
  324           ccomment       = "ID of clump in its host object.";
  325           otype          = GAL_TYPE_INVALID;
  326           ctype          = GAL_TYPE_INT32;
  327           disp_fmt       = 0;
  328           disp_width     = 6;
  329           disp_precision = 0;
  330           /* Is an internal parameter. */
  331           break;
  332 
  333         case UI_KEY_NUMCLUMPS:
  334           name           = "NUM_CLUMPS";
  335           unit           = "counter";
  336           ocomment       = "Number of clumps in this object.";
  337           ccomment       = NULL;
  338           otype          = GAL_TYPE_INT32;
  339           ctype          = GAL_TYPE_INVALID;
  340           disp_fmt       = 0;
  341           disp_width     = 5;
  342           disp_precision = 0;
  343           /* Is an internal parameter. */
  344           break;
  345 
  346         case UI_KEY_AREA:
  347           name           = "AREA";
  348           unit           = "counter";
  349           ocomment       = "Number of non-blank pixels.";
  350           ccomment       = ocomment;
  351           otype          = GAL_TYPE_INT32;
  352           ctype          = GAL_TYPE_INT32;
  353           disp_fmt       = 0;
  354           disp_width     = 6;
  355           disp_precision = 0;
  356           oiflag[ OCOL_NUM ] = ciflag[ CCOL_NUM ] = 1;
  357           break;
  358 
  359         case UI_KEY_CLUMPSAREA:
  360           name           = "AREA_CLUMPS";
  361           unit           = "counter";
  362           ocomment       = "Total number of clump pixels in object.";
  363           ccomment       = NULL;
  364           otype          = GAL_TYPE_INT32;
  365           ctype          = GAL_TYPE_INVALID;
  366           disp_fmt       = 0;
  367           disp_width     = 6;
  368           disp_precision = 0;
  369           oiflag[ OCOL_C_NUM ] = 1;
  370           break;
  371 
  372         case UI_KEY_WEIGHTAREA:
  373           name           = "AREA_WEIGHT";
  374           unit           = "counter";
  375           ocomment       = "Area used for flux-weighted positions.";
  376           ccomment       = ocomment;
  377           otype          = GAL_TYPE_INT32;
  378           ctype          = GAL_TYPE_INT32;
  379           disp_fmt       = 0;
  380           disp_width     = 6;
  381           disp_precision = 0;
  382           oiflag[ OCOL_NUMWHT ] = ciflag[ CCOL_NUMWHT ] = 1;
  383           break;
  384 
  385         case UI_KEY_GEOAREA:
  386           name           = "AREA_FULL";
  387           unit           = "counter";
  388           ocomment       = "Full area of label (irrespective of values).";
  389           ccomment       = ocomment;
  390           otype          = GAL_TYPE_INT32;
  391           ctype          = GAL_TYPE_INT32;
  392           disp_fmt       = 0;
  393           disp_width     = 6;
  394           disp_precision = 0;
  395           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
  396           break;
  397 
  398         case UI_KEY_X:
  399           name           = "X";
  400           unit           = "pixel";
  401           ocomment       = "Flux weighted center (FITS axis 1).";
  402           ccomment       = ocomment;
  403           otype          = GAL_TYPE_FLOAT32;
  404           ctype          = GAL_TYPE_FLOAT32;
  405           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  406           disp_width     = 10;
  407           disp_precision = 3;
  408           oiflag[ OCOL_VX     ] = ciflag[ CCOL_VX     ] = 1;
  409           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
  410           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
  411           oiflag[ OCOL_NUMWHT ] = ciflag[ CCOL_NUMWHT ] = 1;
  412           break;
  413 
  414         case UI_KEY_Y:
  415           name           = "Y";
  416           unit           = "pixel";
  417           ocomment       = "Flux weighted center (FITS axis 2).";
  418           ccomment       = ocomment;
  419           otype          = GAL_TYPE_FLOAT32;
  420           ctype          = GAL_TYPE_FLOAT32;
  421           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  422           disp_width     = 10;
  423           disp_precision = 3;
  424           oiflag[ OCOL_VY     ] = ciflag[ CCOL_VY     ] = 1;
  425           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
  426           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
  427           oiflag[ OCOL_NUMWHT ] = ciflag[ CCOL_NUMWHT ] = 1;
  428           break;
  429 
  430         case UI_KEY_GEOX:
  431           name           = "GEO_X";
  432           unit           = "pixel";
  433           ocomment       = "Geometric center (FITS axis 1).";
  434           ccomment       = ocomment;
  435           otype          = GAL_TYPE_FLOAT32;
  436           ctype          = GAL_TYPE_FLOAT32;
  437           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  438           disp_width     = 10;
  439           disp_precision = 3;
  440           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
  441           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
  442           break;
  443 
  444         case UI_KEY_GEOY:
  445           name           = "GEO_Y";
  446           unit           = "pixel";
  447           ocomment       = "Geometric center (FITS axis 2).";
  448           ccomment       = ocomment;
  449           otype          = GAL_TYPE_FLOAT32;
  450           ctype          = GAL_TYPE_FLOAT32;
  451           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  452           disp_width     = 10;
  453           disp_precision = 3;
  454           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
  455           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
  456           break;
  457 
  458         case UI_KEY_CLUMPSX:
  459           name           = "CLUMPS_X";
  460           unit           = "pixel";
  461           ocomment       = "Flux weighted center of clumps (FITS axis 1).";
  462           ccomment       = NULL;
  463           otype          = GAL_TYPE_FLOAT32;
  464           ctype          = GAL_TYPE_INVALID;
  465           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  466           disp_width     = 10;
  467           disp_precision = 3;
  468           oiflag[ OCOL_C_VX     ] = 1;
  469           oiflag[ OCOL_C_GX     ] = 1;
  470           oiflag[ OCOL_C_SUMWHT ] = 1;
  471           oiflag[ OCOL_C_NUMWHT ] = 1;
  472           break;
  473 
  474         case UI_KEY_CLUMPSY:
  475           name           = "CLUMPS_Y";
  476           unit           = "pixel";
  477           ocomment       = "Flux weighted center of clumps (FITS axis 2).";
  478           ccomment       = NULL;
  479           otype          = GAL_TYPE_FLOAT32;
  480           ctype          = GAL_TYPE_INVALID;
  481           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  482           disp_width     = 10;
  483           disp_precision = 3;
  484           oiflag[ OCOL_C_VY     ] = 1;
  485           oiflag[ OCOL_C_GY     ] = 1;
  486           oiflag[ OCOL_C_SUMWHT ] = 1;
  487           oiflag[ OCOL_C_NUMWHT ] = 1;
  488           break;
  489 
  490         case UI_KEY_CLUMPSGEOX:
  491           name           = "CLUMPS_GEO_X";
  492           unit           = "pixel";
  493           ocomment       = "Geometric center of clumps (FITS axis 1).";
  494           ccomment       = NULL;
  495           otype          = GAL_TYPE_FLOAT32;
  496           ctype          = GAL_TYPE_INVALID;
  497           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  498           disp_width     = 10;
  499           disp_precision = 3;
  500           oiflag[ OCOL_C_GX     ] = 1;
  501           oiflag[ OCOL_C_NUMALL ] = 1;
  502           break;
  503 
  504         case UI_KEY_CLUMPSGEOY:
  505           name           = "CLUMPS_GEO_Y";
  506           unit           = "pixel";
  507           ocomment       = "Geometric center of clumps (FITS axis 2).";
  508           ccomment       = NULL;
  509           otype          = GAL_TYPE_FLOAT32;
  510           ctype          = GAL_TYPE_INVALID;
  511           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  512           disp_width     = 10;
  513           disp_precision = 3;
  514           oiflag[ OCOL_C_GY     ] = 1;
  515           oiflag[ OCOL_C_NUMALL ] = 1;
  516           break;
  517 
  518         case UI_KEY_MINX:
  519           name           = "MIN_X";
  520           unit           = "pixel";
  521           ocomment       = "Minimum X axis pixel position.";
  522           ccomment       = ocomment;
  523           otype          = GAL_TYPE_UINT32;
  524           ctype          = GAL_TYPE_UINT32;
  525           disp_fmt       = 0;
  526           disp_width     = 10;
  527           disp_precision = 0;
  528           ciflag[ CCOL_MINX ] = 1;
  529           break;
  530 
  531         case UI_KEY_MAXX:
  532           name           = "MAX_X";
  533           unit           = "pixel";
  534           ocomment       = "Maximum X axis pixel position.";
  535           ccomment       = ocomment;
  536           otype          = GAL_TYPE_UINT32;
  537           ctype          = GAL_TYPE_UINT32;
  538           disp_fmt       = 0;
  539           disp_width     = 10;
  540           disp_precision = 0;
  541           ciflag[ CCOL_MAXX ] = 1;
  542           break;
  543 
  544         case UI_KEY_MINY:
  545           name           = "MIN_Y";
  546           unit           = "pixel";
  547           ocomment       = "Minimum Y axis pixel position.";
  548           ccomment       = ocomment;
  549           otype          = GAL_TYPE_UINT32;
  550           ctype          = GAL_TYPE_UINT32;
  551           disp_fmt       = 0;
  552           disp_width     = 10;
  553           disp_precision = 0;
  554           ciflag[ CCOL_MINY ] = 1;
  555           break;
  556 
  557         case UI_KEY_MAXY:
  558           name           = "MAX_Y";
  559           unit           = "pixel";
  560           ocomment       = "Maximum Y axis pixel position.";
  561           ccomment       = ocomment;
  562           otype          = GAL_TYPE_UINT32;
  563           ctype          = GAL_TYPE_UINT32;
  564           disp_fmt       = 0;
  565           disp_width     = 10;
  566           disp_precision = 0;
  567           ciflag[ CCOL_MAXY ] = 1;
  568           break;
  569 
  570         case UI_KEY_W1:
  571           name           = p->ctype[0];
  572           unit           = p->objects->wcs->cunit[0];
  573           ocomment       = "Flux weighted center (WCS axis 1).";
  574           ccomment       = ocomment;
  575           otype          = GAL_TYPE_FLOAT64;
  576           ctype          = GAL_TYPE_FLOAT64;
  577           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  578           disp_width     = 13;
  579           disp_precision = 7;
  580           columns_alloc_radec(p);
  581           oiflag[ OCOL_VX     ] = ciflag[ CCOL_VX     ] = 1;
  582           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
  583           oiflag[ OCOL_VY     ] = ciflag[ CCOL_VY     ] = 1;
  584           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
  585           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
  586           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
  587           break;
  588 
  589         case UI_KEY_W2:
  590           name           = p->ctype[1];
  591           unit           = p->objects->wcs->cunit[1];
  592           ocomment       = "Flux weighted center (WCS axis 2).";
  593           ccomment       = ocomment;
  594           otype          = GAL_TYPE_FLOAT64;
  595           ctype          = GAL_TYPE_FLOAT64;
  596           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  597           disp_width     = 13;
  598           disp_precision = 7;
  599           columns_alloc_radec(p);
  600           oiflag[ OCOL_VX     ] = ciflag[ CCOL_VX     ] = 1;
  601           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
  602           oiflag[ OCOL_VY     ] = ciflag[ CCOL_VY     ] = 1;
  603           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
  604           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
  605           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
  606           break;
  607 
  608         case UI_KEY_GEOW1:
  609           name           = gal_checkset_malloc_cat("GEO_", p->ctype[0]);
  610           unit           = p->objects->wcs->cunit[0];
  611           ocomment       = "Geometric center (WCS axis 1).";
  612           ccomment       = ocomment;
  613           otype          = GAL_TYPE_FLOAT64;
  614           ctype          = GAL_TYPE_FLOAT64;
  615           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  616           disp_width     = 13;
  617           disp_precision = 7;
  618           columns_alloc_georadec(p);
  619           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
  620           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
  621           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
  622           break;
  623 
  624         case UI_KEY_GEOW2:
  625           name           = gal_checkset_malloc_cat("GEO_", p->ctype[1]);
  626           unit           = p->objects->wcs->cunit[1];
  627           ocomment       = "Geometric center (WCS axis 2).";
  628           ccomment       = ocomment;
  629           otype          = GAL_TYPE_FLOAT64;
  630           ctype          = GAL_TYPE_FLOAT64;
  631           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  632           disp_width     = 13;
  633           disp_precision = 7;
  634           columns_alloc_georadec(p);
  635           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
  636           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
  637           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
  638           break;
  639 
  640         case UI_KEY_CLUMPSW1:
  641           name           = gal_checkset_malloc_cat("CLUMPS_", p->ctype[0]);
  642           unit           = p->objects->wcs->cunit[0];
  643           ocomment       = "Flux.wht center of all clumps (WCS axis 1).";
  644           ccomment       = NULL;
  645           otype          = GAL_TYPE_FLOAT64;
  646           ctype          = GAL_TYPE_INVALID;
  647           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  648           disp_width     = 13;
  649           disp_precision = 7;
  650           columns_alloc_clumpsradec(p);
  651           oiflag[ OCOL_C_VX     ] = 1;
  652           oiflag[ OCOL_C_VY     ] = 1;
  653           oiflag[ OCOL_C_GX     ] = 1;
  654           oiflag[ OCOL_C_GY     ] = 1;
  655           oiflag[ OCOL_C_SUMWHT ] = 1;
  656           oiflag[ OCOL_C_NUMALL ] = 1;
  657           break;
  658 
  659         case UI_KEY_CLUMPSW2:
  660           name           = gal_checkset_malloc_cat("CLUMPS_", p->ctype[1]);
  661           unit           = p->objects->wcs->cunit[1];
  662           ocomment       = "Flux.wht center of all clumps (WCS axis 2).";
  663           ccomment       = NULL;
  664           otype          = GAL_TYPE_FLOAT64;
  665           ctype          = GAL_TYPE_INVALID;
  666           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  667           disp_width     = 15;
  668           disp_precision = 7;
  669           columns_alloc_clumpsradec(p);
  670           oiflag[ OCOL_C_VX     ] = 1;
  671           oiflag[ OCOL_C_VY     ] = 1;
  672           oiflag[ OCOL_C_GX     ] = 1;
  673           oiflag[ OCOL_C_GY     ] = 1;
  674           oiflag[ OCOL_C_SUMWHT ] = 1;
  675           oiflag[ OCOL_C_NUMALL ] = 1;
  676           break;
  677 
  678         case UI_KEY_CLUMPSGEOW1:
  679           name           = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[0]);
  680           unit           = p->objects->wcs->cunit[0];
  681           ocomment       = "Geometric center of all clumps (WCS axis 1).";
  682           ccomment       = NULL;
  683           otype          = GAL_TYPE_FLOAT64;
  684           ctype          = GAL_TYPE_INVALID;
  685           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  686           disp_width     = 13;
  687           disp_precision = 7;
  688           columns_alloc_clumpsgeoradec(p);
  689           oiflag[ OCOL_C_GX     ] = 1;
  690           oiflag[ OCOL_C_GY     ] = 1;
  691           oiflag[ OCOL_C_NUMALL ] = 1;
  692           break;
  693 
  694         case UI_KEY_CLUMPSGEOW2:
  695           name           = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[1]);
  696           unit           = p->objects->wcs->cunit[1];
  697           ocomment       = "Geometric center of all clumps (WCS axis 2).";
  698           ccomment       = NULL;
  699           otype          = GAL_TYPE_FLOAT64;
  700           ctype          = GAL_TYPE_INVALID;
  701           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  702           disp_width     = 13;
  703           disp_precision = 7;
  704           columns_alloc_clumpsgeoradec(p);
  705           oiflag[ OCOL_C_GX     ] = 1;
  706           oiflag[ OCOL_C_GY     ] = 1;
  707           oiflag[ OCOL_C_NUMALL ] = 1;
  708           break;
  709 
  710         case UI_KEY_BRIGHTNESS:
  711           name           = "BRIGHTNESS";
  712           unit           = MKCATALOG_NO_UNIT;
  713           ocomment       = "Brightness (sum of sky subtracted values).";
  714           ccomment       = "Brightness (sum of pixels subtracted by rivers).";
  715           otype          = GAL_TYPE_FLOAT32;
  716           ctype          = GAL_TYPE_FLOAT32;
  717           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  718           disp_width     = 10;
  719           disp_precision = 5;
  720           oiflag[ OCOL_NUM     ] = ciflag[ CCOL_NUM     ] = 1;
  721           oiflag[ OCOL_SUM     ] = ciflag[ CCOL_SUM     ] = 1;
  722                                    ciflag[ CCOL_RIV_NUM ] = 1;
  723                                    ciflag[ CCOL_RIV_SUM ] = 1;
  724           break;
  725 
  726         case UI_KEY_BRIGHTNESSERR:
  727           name           = "BRIGHTNESS_ERROR";
  728           unit           = MKCATALOG_NO_UNIT;
  729           ocomment       = "Error (1-sigma) in measuring brightness.";
  730           ccomment       = ocomment;
  731           otype          = GAL_TYPE_FLOAT32;
  732           ctype          = GAL_TYPE_FLOAT32;
  733           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  734           disp_width     = 10;
  735           disp_precision = 5;
  736           oiflag[ OCOL_NUM     ] = ciflag[ CCOL_NUM         ] = 1;
  737           oiflag[ OCOL_SUM_VAR ] = ciflag[ CCOL_SUM_VAR     ] = 1;
  738                                    ciflag[ CCOL_RIV_NUM     ] = 1;
  739                                    ciflag[ CCOL_RIV_SUM_VAR ] = 1;
  740           break;
  741 
  742         case UI_KEY_CLUMPSBRIGHTNESS:
  743           name           = "CLUMPS_BRIGHTNESS";
  744           unit           = MKCATALOG_NO_UNIT;
  745           ocomment       = "Brightness (sum of pixel values) in clumps.";
  746           ccomment       = NULL;
  747           otype          = GAL_TYPE_FLOAT32;
  748           ctype          = GAL_TYPE_INVALID;
  749           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  750           disp_width     = 10;
  751           disp_precision = 5;
  752           oiflag[ OCOL_C_NUM ] = 1;
  753           oiflag[ OCOL_C_SUM ] = 1;
  754           break;
  755 
  756         case UI_KEY_BRIGHTNESSNORIVER:
  757           name           = "NO_RIVER_BRIGHTNESS";
  758           unit           = MKCATALOG_NO_UNIT;
  759           ocomment       = NULL;
  760           ccomment       = "Brightness (sum of sky subtracted values).";
  761           otype          = GAL_TYPE_INVALID;
  762           ctype          = GAL_TYPE_FLOAT32;
  763           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  764           disp_width     = 10;
  765           disp_precision = 5;
  766           ciflag[ CCOL_NUM ] = 1;
  767           ciflag[ CCOL_SUM ] = 1;
  768           break;
  769 
  770         case UI_KEY_MEAN:
  771           name           = "MEAN";
  772           unit           = MKCATALOG_NO_UNIT;
  773           ocomment       = "Mean of sky subtracted values.";
  774           ccomment       = "Mean of pixels subtracted by rivers.";
  775           otype          = GAL_TYPE_FLOAT32;
  776           ctype          = GAL_TYPE_FLOAT32;
  777           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  778           disp_width     = 10;
  779           disp_precision = 5;
  780           oiflag[ OCOL_NUM     ] = ciflag[ CCOL_NUM     ] = 1;
  781           oiflag[ OCOL_SUM     ] = ciflag[ CCOL_SUM     ] = 1;
  782                                    ciflag[ CCOL_RIV_NUM ] = 1;
  783                                    ciflag[ CCOL_RIV_SUM ] = 1;
  784           break;
  785 
  786         case UI_KEY_MEDIAN:
  787           name           = "MEDIAN";
  788           unit           = MKCATALOG_NO_UNIT;
  789           ocomment       = "Median of sky subtracted values.";
  790           ccomment       = "Median of pixels subtracted by rivers.";
  791           otype          = GAL_TYPE_FLOAT32;
  792           ctype          = GAL_TYPE_FLOAT32;
  793           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  794           disp_width     = 10;
  795           disp_precision = 5;
  796           oiflag[ OCOL_NUM     ] = ciflag[ CCOL_NUM     ] = 1;
  797           oiflag[ OCOL_MEDIAN  ] = ciflag[ CCOL_MEDIAN  ] = 1;
  798                                    ciflag[ CCOL_RIV_NUM ] = 1;
  799                                    ciflag[ CCOL_RIV_SUM ] = 1;
  800           break;
  801 
  802         case UI_KEY_MAGNITUDE:
  803           name           = "MAGNITUDE";
  804           unit           = "log";
  805           ocomment       = "Magnitude.";
  806           ccomment       = ocomment;
  807           otype          = GAL_TYPE_FLOAT32;
  808           ctype          = GAL_TYPE_FLOAT32;
  809           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  810           disp_width     = 8;
  811           disp_precision = 3;
  812           p->hasmag      = 1;
  813           oiflag[ OCOL_NUM     ] = ciflag[ CCOL_NUM     ] = 1;
  814           oiflag[ OCOL_SUM     ] = ciflag[ CCOL_SUM     ] = 1;
  815                                    ciflag[ CCOL_SUM     ] = 1;
  816                                    ciflag[ CCOL_RIV_NUM ] = 1;
  817                                    ciflag[ CCOL_RIV_SUM ] = 1;
  818                                    ciflag[ CCOL_RIV_NUM ] = 1;
  819           break;
  820 
  821         case UI_KEY_MAGNITUDEERR:
  822           name           = "MAGNITUDE_ERROR";
  823           unit           = "log";
  824           ocomment       = "Error in measuring magnitude.";
  825           ccomment       = ocomment;
  826           otype          = GAL_TYPE_FLOAT32;
  827           ctype          = GAL_TYPE_FLOAT32;
  828           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  829           disp_width     = 8;
  830           disp_precision = 3;
  831           oiflag[ OCOL_SUM         ] = ciflag[ CCOL_SUM         ] = 1;
  832           oiflag[ OCOL_SUM_VAR     ] = ciflag[ CCOL_SUM_VAR     ] = 1;
  833                                        ciflag[ CCOL_RIV_SUM     ] = 1;
  834                                        ciflag[ CCOL_RIV_SUM_VAR ] = 1;
  835           break;
  836 
  837         case UI_KEY_CLUMPSMAGNITUDE:
  838           name           = "CLUMPS_MAGNITUDE";
  839           unit           = "log";
  840           ocomment       = "Magnitude in all clumps.";
  841           ccomment       = NULL;
  842           otype          = GAL_TYPE_FLOAT32;
  843           ctype          = GAL_TYPE_INVALID;
  844           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  845           disp_width     = 8;
  846           disp_precision = 3;
  847           p->hasmag      = 1;
  848           oiflag[ OCOL_SUM         ] = ciflag[ CCOL_SUM         ] = 1;
  849           oiflag[ OCOL_SUM_VAR     ] = ciflag[ CCOL_SUM_VAR     ] = 1;
  850                                        ciflag[ CCOL_RIV_NUM     ] = 1;
  851                                        ciflag[ CCOL_RIV_SUM     ] = 1;
  852                                        ciflag[ CCOL_RIV_SUM_VAR ] = 1;
  853           break;
  854 
  855         case UI_KEY_UPPERLIMIT:
  856           name           = "UPPERLIMIT";
  857           unit           = MKCATALOG_NO_UNIT;
  858           ocomment       = "Upper limit value (random positionings).";
  859           ccomment       = ocomment;
  860           otype          = GAL_TYPE_FLOAT32;
  861           ctype          = GAL_TYPE_FLOAT32;
  862           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  863           disp_width     = 10;
  864           disp_precision = 5;
  865           p->upperlimit  = 1;
  866           oiflag[ OCOL_UPPERLIMIT_B ] = ciflag[ CCOL_UPPERLIMIT_B ] = 1;
  867 
  868           break;
  869 
  870         case UI_KEY_UPPERLIMITMAG:
  871           name           = "UPPERLIMIT_MAG";
  872           unit           = "log";
  873           ocomment       = "Upper limit magnitude (random positionings).";
  874           ccomment       = ocomment;
  875           otype          = GAL_TYPE_FLOAT32;
  876           ctype          = GAL_TYPE_FLOAT32;
  877           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  878           disp_width     = 8;
  879           disp_precision = 3;
  880           p->hasmag      = 1;
  881           p->upperlimit  = 1;
  882           oiflag[ OCOL_UPPERLIMIT_B ] = ciflag[ CCOL_UPPERLIMIT_B ] = 1;
  883           break;
  884 
  885         case UI_KEY_UPPERLIMITONESIGMA:
  886           name           = "UPPERLIMIT_ONE_SIGMA";
  887           unit           = MKCATALOG_NO_UNIT;
  888           ocomment       = "One sigma value of all random measurements.";
  889           ccomment       = ocomment;
  890           otype          = GAL_TYPE_FLOAT32;
  891           ctype          = GAL_TYPE_FLOAT32;
  892           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  893           disp_width     = 10;
  894           disp_precision = 5;
  895           p->upperlimit  = 1;
  896           oiflag[ OCOL_UPPERLIMIT_S ] = ciflag[ CCOL_UPPERLIMIT_S ] = 1;
  897           break;
  898 
  899         case UI_KEY_UPPERLIMITSIGMA:
  900           name           = "UPPERLIMIT_SIGMA";
  901           unit           = "frac";
  902           ocomment       = "Place in upperlimit distribution (sigma "
  903                            "multiple).";
  904           ccomment       = ocomment;
  905           otype          = GAL_TYPE_FLOAT32;
  906           ctype          = GAL_TYPE_FLOAT32;
  907           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  908           disp_width     = 10;
  909           disp_precision = 5;
  910           p->upperlimit  = 1;
  911           oiflag[ OCOL_NUM          ] = ciflag[ CCOL_NUM          ] = 1;
  912           oiflag[ OCOL_SUM          ] = ciflag[ CCOL_SUM          ] = 1;
  913           oiflag[ OCOL_UPPERLIMIT_S ] = ciflag[ CCOL_UPPERLIMIT_S ] = 1;
  914                                         ciflag[ CCOL_RIV_NUM      ] = 1;
  915                                         ciflag[ CCOL_RIV_SUM      ] = 1;
  916           break;
  917 
  918         case UI_KEY_UPPERLIMITQUANTILE:
  919           name           = "UPPERLIMIT_QUANTILE";
  920           unit           = "quantile";
  921           ocomment       = "Quantile of brightness in random distribution.";
  922           ccomment       = ocomment;
  923           otype          = GAL_TYPE_FLOAT32;
  924           ctype          = GAL_TYPE_FLOAT32;
  925           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  926           disp_width     = 10;
  927           disp_precision = 5;
  928           p->upperlimit  = 1;
  929           oiflag[ OCOL_UPPERLIMIT_Q ] = ciflag[ CCOL_UPPERLIMIT_Q ] = 1;
  930           break;
  931 
  932         case UI_KEY_UPPERLIMITSKEW:
  933           name           = "UPPERLIMIT_SKEW";
  934           unit           = "frac";
  935           ocomment       = "(Mean-Median)/STD of random distribution.";
  936           ccomment       = ocomment;
  937           otype          = GAL_TYPE_FLOAT32;
  938           ctype          = GAL_TYPE_FLOAT32;
  939           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  940           disp_width     = 8;
  941           disp_precision = 3;
  942           p->upperlimit  = 1;
  943           oiflag[ OCOL_UPPERLIMIT_SKEW ] = oiflag[ CCOL_UPPERLIMIT_SKEW ] = 1;
  944           break;
  945 
  946         case UI_KEY_RIVERAVE:
  947           name           = "RIVER_AVE";
  948           unit           = MKCATALOG_NO_UNIT;
  949           ocomment       = NULL;
  950           ccomment       = "Average river value surrounding this clump.";
  951           otype          = GAL_TYPE_INVALID;
  952           ctype          = GAL_TYPE_FLOAT32;
  953           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  954           disp_width     = 10;
  955           disp_precision = 5;
  956           ciflag[ CCOL_RIV_NUM ] = ciflag[ CCOL_RIV_SUM ] = 1;
  957           break;
  958 
  959         case UI_KEY_RIVERNUM:
  960           name           = "RIVER_NUM";
  961           unit           = "counter";
  962           ocomment       = NULL;
  963           ccomment       = "Number of river pixels around this clump.";
  964           otype          = GAL_TYPE_INVALID;
  965           ctype          = GAL_TYPE_INT32;
  966           disp_fmt       = 0;
  967           disp_width     = 5;
  968           disp_precision = 0;
  969           ciflag[ CCOL_RIV_NUM ] = 1;
  970           break;
  971 
  972         case UI_KEY_SN:
  973           name           = "SN";
  974           unit           = "ratio";
  975           ocomment       = "Signal to noise ratio.";
  976           ccomment       = ocomment;
  977           otype          = GAL_TYPE_FLOAT32;
  978           ctype          = GAL_TYPE_FLOAT32;
  979           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
  980           disp_width     = 10;
  981           disp_precision = 3;
  982           oiflag[ OCOL_SUM         ] = ciflag[ CCOL_SUM         ] = 1;
  983           oiflag[ OCOL_SUM_VAR     ] = ciflag[ CCOL_SUM_VAR     ] = 1;
  984                                        ciflag[ CCOL_RIV_NUM     ] = 1;
  985                                        ciflag[ CCOL_RIV_SUM     ] = 1;
  986                                        ciflag[ CCOL_RIV_SUM_VAR ] = 1;
  987           break;
  988 
  989         case UI_KEY_SKY:
  990           name           = "SKY";
  991           unit           = MKCATALOG_NO_UNIT;
  992           ocomment       = "Average input sky value.";
  993           ccomment       = ocomment;
  994           otype          = GAL_TYPE_FLOAT32;
  995           ctype          = GAL_TYPE_FLOAT32;
  996           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
  997           disp_width     = 10;
  998           disp_precision = 4;
  999           oiflag[ OCOL_NUM    ] = ciflag[ CCOL_NUM    ] = 1;
 1000           oiflag[ OCOL_SUMSKY ] = ciflag[ CCOL_SUMSKY ] = 1;
 1001           break;
 1002 
 1003         case UI_KEY_STD:
 1004           name           = "STD";
 1005           unit           = MKCATALOG_NO_UNIT;
 1006           ocomment       = "Sky standard deviation under object.";
 1007           ccomment       = ocomment;
 1008           otype          = GAL_TYPE_FLOAT32;
 1009           ctype          = GAL_TYPE_FLOAT32;
 1010           disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
 1011           disp_width     = 10;
 1012           disp_precision = 4;
 1013           oiflag[ OCOL_NUM    ] = ciflag[ CCOL_NUM    ] = 1;
 1014           oiflag[ OCOL_SUMVAR ] = ciflag[ CCOL_SUMVAR ] = 1;
 1015           break;
 1016 
 1017         case UI_KEY_SEMIMAJOR:
 1018           name           = "SEMI_MAJOR";
 1019           unit           = "pixel";
 1020           ocomment       = "Flux weighted semi-major axis.";
 1021           ccomment       = ocomment;
 1022           otype          = GAL_TYPE_FLOAT32;
 1023           ctype          = GAL_TYPE_FLOAT32;
 1024           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1025           disp_width     = 10;
 1026           disp_precision = 3;
 1027           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
 1028           oiflag[ OCOL_VX     ] = ciflag[ CCOL_VX     ] = 1;
 1029           oiflag[ OCOL_VY     ] = ciflag[ CCOL_VY     ] = 1;
 1030           oiflag[ OCOL_VXX    ] = ciflag[ CCOL_VXX    ] = 1;
 1031           oiflag[ OCOL_VYY    ] = ciflag[ CCOL_VYY    ] = 1;
 1032           oiflag[ OCOL_VXY    ] = ciflag[ CCOL_VXY    ] = 1;
 1033           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1034           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1035           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1036           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1037           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1038           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1039           break;
 1040 
 1041         case UI_KEY_SEMIMINOR:
 1042           name           = "SEMI_MINOR";
 1043           unit           = "pixel";
 1044           ocomment       = "Flux weighted semi-minor axis.";
 1045           ccomment       = ocomment;
 1046           otype          = GAL_TYPE_FLOAT32;
 1047           ctype          = GAL_TYPE_FLOAT32;
 1048           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1049           disp_width     = 10;
 1050           disp_precision = 3;
 1051           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
 1052           oiflag[ OCOL_VX     ] = ciflag[ CCOL_VX     ] = 1;
 1053           oiflag[ OCOL_VY     ] = ciflag[ CCOL_VY     ] = 1;
 1054           oiflag[ OCOL_VXX    ] = ciflag[ CCOL_VXX    ] = 1;
 1055           oiflag[ OCOL_VYY    ] = ciflag[ CCOL_VYY    ] = 1;
 1056           oiflag[ OCOL_VXY    ] = ciflag[ CCOL_VXY    ] = 1;
 1057           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1058           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1059           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1060           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1061           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1062           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1063           break;
 1064 
 1065         case UI_KEY_AXISRATIO:
 1066           name           = "AXIS_RATIO";
 1067           unit           = "ratio";
 1068           ocomment       = "Flux weighted axis ratio (minor/major).";
 1069           ccomment       = ocomment;
 1070           otype          = GAL_TYPE_FLOAT32;
 1071           ctype          = GAL_TYPE_FLOAT32;
 1072           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1073           disp_width     = 7;
 1074           disp_precision = 3;
 1075           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
 1076           oiflag[ OCOL_VX     ] = ciflag[ CCOL_VX     ] = 1;
 1077           oiflag[ OCOL_VY     ] = ciflag[ CCOL_VY     ] = 1;
 1078           oiflag[ OCOL_VXX    ] = ciflag[ CCOL_VXX    ] = 1;
 1079           oiflag[ OCOL_VYY    ] = ciflag[ CCOL_VYY    ] = 1;
 1080           oiflag[ OCOL_VXY    ] = ciflag[ CCOL_VXY    ] = 1;
 1081           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1082           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1083           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1084           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1085           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1086           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1087           break;
 1088 
 1089         case UI_KEY_POSITIONANGLE:
 1090           name           = "POSITION_ANGLE";
 1091           unit           = "degrees";
 1092           ocomment       = "Position angle.";
 1093           ccomment       = ocomment;
 1094           otype          = GAL_TYPE_FLOAT32;
 1095           ctype          = GAL_TYPE_FLOAT32;
 1096           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1097           disp_width     = 10;
 1098           disp_precision = 3;
 1099           oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
 1100           oiflag[ OCOL_VX     ] = ciflag[ CCOL_VX     ] = 1;
 1101           oiflag[ OCOL_VY     ] = ciflag[ CCOL_VY     ] = 1;
 1102           oiflag[ OCOL_VXX    ] = ciflag[ CCOL_VXX    ] = 1;
 1103           oiflag[ OCOL_VYY    ] = ciflag[ CCOL_VYY    ] = 1;
 1104           oiflag[ OCOL_VXY    ] = ciflag[ CCOL_VXY    ] = 1;
 1105           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1106           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1107           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1108           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1109           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1110           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1111           break;
 1112 
 1113         case UI_KEY_GEOSEMIMAJOR:
 1114           name           = "GEO_SEMI_MAJOR";
 1115           unit           = "pixel";
 1116           ocomment       = "Geometric semi-major axis.";
 1117           ccomment       = ocomment;
 1118           otype          = GAL_TYPE_FLOAT32;
 1119           ctype          = GAL_TYPE_FLOAT32;
 1120           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1121           disp_width     = 10;
 1122           disp_precision = 3;
 1123           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1124           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1125           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1126           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1127           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1128           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1129           break;
 1130 
 1131         case UI_KEY_GEOSEMIMINOR:
 1132           name           = "GEO_SEMI_MINOR";
 1133           unit           = "pixel";
 1134           ocomment       = "Geometric semi-minor axis.";
 1135           ccomment       = ocomment;
 1136           otype          = GAL_TYPE_FLOAT32;
 1137           ctype          = GAL_TYPE_FLOAT32;
 1138           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1139           disp_width     = 10;
 1140           disp_precision = 3;
 1141           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1142           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1143           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1144           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1145           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1146           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1147           break;
 1148 
 1149         case UI_KEY_GEOAXISRATIO:
 1150           name           = "GEO_AXIS_RATIO";
 1151           unit           = "ratio";
 1152           ocomment       = "Geometric axis ratio (minor/major).";
 1153           ccomment       = ocomment;
 1154           otype          = GAL_TYPE_FLOAT32;
 1155           ctype          = GAL_TYPE_FLOAT32;
 1156           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1157           disp_width     = 7;
 1158           disp_precision = 3;
 1159           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1160           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1161           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1162           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1163           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1164           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1165           break;
 1166 
 1167         case UI_KEY_GEOPOSITIONANGLE:
 1168           name           = "GEO_POSITION_ANGLE";
 1169           unit           = "degrees";
 1170           ocomment       = "Geometric Position angle.";
 1171           ccomment       = ocomment;
 1172           otype          = GAL_TYPE_FLOAT32;
 1173           ctype          = GAL_TYPE_FLOAT32;
 1174           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
 1175           disp_width     = 10;
 1176           disp_precision = 3;
 1177           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
 1178           oiflag[ OCOL_GX     ] = ciflag[ CCOL_GX     ] = 1;
 1179           oiflag[ OCOL_GY     ] = ciflag[ CCOL_GY     ] = 1;
 1180           oiflag[ OCOL_GXX    ] = ciflag[ CCOL_GXX    ] = 1;
 1181           oiflag[ OCOL_GYY    ] = ciflag[ CCOL_GYY    ] = 1;
 1182           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
 1183           break;
 1184 
 1185         default:
 1186           error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix "
 1187                 "the problem. The code %d is not an internally recognized "
 1188                 "column code", __func__, PACKAGE_BUGREPORT, colcode->v);
 1189         }
 1190 
 1191 
 1192       /* If this is an objects column, add it to the list of columns. We
 1193          will be using the `status' element to keep the MakeCatalog code
 1194          for the columns. */
 1195       if(otype!=GAL_TYPE_INVALID)
 1196         {
 1197           gal_list_data_add_alloc(&p->objectcols, NULL, otype, 1,
 1198                                   &p->numobjects, NULL, 0, p->cp.minmapsize,
 1199                                   name, unit, ocomment);
 1200           p->objectcols->status         = colcode->v;
 1201           p->objectcols->disp_fmt       = disp_fmt;
 1202           p->objectcols->disp_width     = disp_width;
 1203           p->objectcols->disp_precision = disp_precision;
 1204         }
 1205 
 1206 
 1207       /* Similar to the objects column above but for clumps, but since the
 1208          clumps image is optional, we need a further check before actually
 1209          allocating the column. */
 1210       if(ctype!=GAL_TYPE_INVALID)
 1211         {
 1212           /* If a clumps labeled image, add this column for the output. */
 1213           if(p->clumps)
 1214             {
 1215               gal_list_data_add_alloc(&p->clumpcols, NULL, ctype, 1,
 1216                                       &p->numclumps, NULL, 0,
 1217                                       p->cp.minmapsize, name, unit, ccomment);
 1218               p->clumpcols->status         = colcode->v;
 1219               p->clumpcols->disp_fmt       = disp_fmt;
 1220               p->clumpcols->disp_width     = disp_width;
 1221               p->clumpcols->disp_precision = disp_precision;
 1222             }
 1223 
 1224 
 1225           /* If this is a clumps-only column and no clumps image was
 1226              given. Add the column to the list of similar columns to inform
 1227              the user.
 1228 
 1229              We'll just ignore the clump-specific ID-related columns,
 1230              because the `--ids' (generic for both objects and clumps) is a
 1231              simple generic solution for identifiers. Also, ultimately,
 1232              they aren't measurements. */
 1233           else if( otype==GAL_TYPE_INVALID
 1234                    && colcode->v!=UI_KEY_HOSTOBJID
 1235                    && colcode->v!=UI_KEY_IDINHOSTOBJ )
 1236             gal_list_str_add(&noclumpimg, name, 1);
 1237         }
 1238     }
 1239 
 1240 
 1241   /* If the user has asked for clump-only columns, but no clumps catalog is
 1242      to be created (the `--clumpscat' option was not given or there were no
 1243      clumps in the specified image), then print an informative message that
 1244      the columns in question will be ignored. */
 1245   if(noclumpimg)
 1246     {
 1247       gal_list_str_reverse(&noclumpimg);
 1248       fprintf(stderr, "WARNING: the following column(s) are unique to "
 1249               "clumps (not objects), but the `--clumpscat' option has not "
 1250               "been called, or there were no clumps in the clumps labeled "
 1251               "image. Hence, these columns will be ignored in the "
 1252               "output.\n\n");
 1253       for(strtmp=noclumpimg; strtmp!=NULL; strtmp=strtmp->next)
 1254         fprintf(stderr, "\t%s\n", strtmp->v);
 1255       gal_list_str_free(noclumpimg, 1);
 1256       fprintf(stderr, "\n-------\n");
 1257     }
 1258 
 1259 
 1260   /* Free the general columns information because it is no longe needed,
 1261      we'll set it back to NULL afterwards so it is not mistakenly used. */
 1262   gal_list_i32_free(p->columnids);
 1263   p->columnids=NULL;
 1264 }
 1265 
 1266 
 1267 
 1268 
 1269 
 1270 
 1271 
 1272 
 1273 
 1274 
 1275 
 1276 
 1277 
 1278 
 1279 
 1280 
 1281 
 1282 
 1283 
 1284 
 1285 /******************************************************************/
 1286 /**********            Column calculation           ***************/
 1287 /******************************************************************/
 1288 #define MKC_RATIO(TOP,BOT) ( (BOT)!=0.0f ? (TOP)/(BOT) : NAN )
 1289 #define MKC_MAG(B)         ( ((B)>0) ? -2.5f * log10(B) + p->zeropoint : NAN )
 1290 
 1291 
 1292 
 1293 
 1294 
 1295 /* Calculate the error in brightness. */
 1296 static double
 1297 columns_brightness_error(struct mkcatalogparams *p, double *row, int o0c1)
 1298 {
 1299   double V = row[ o0c1 ? CCOL_SUM_VAR : OCOL_SUM_VAR ];
 1300   double OV = (o0c1 && row[ CCOL_RIV_NUM ]) ? row[ CCOL_RIV_SUM_VAR ] : 0.0;
 1301   return sqrt(V+OV);
 1302 }
 1303 
 1304 
 1305 
 1306 
 1307 
 1308 /* Calculate the Signal to noise ratio for the object. */
 1309 static double
 1310 columns_sn(struct mkcatalogparams *p, double *row, int o0c1)
 1311 {
 1312   double I = row[ o0c1 ? CCOL_SUM     : OCOL_SUM     ];
 1313 
 1314   /* When grown clumps are requested from NoiseChisel, some "clumps" will
 1315      completely cover their objects and there will be no rivers. So if this
 1316      is a clump, and the river area is 0, we should treat the S/N as a an
 1317      object. */
 1318   double O = (o0c1 && row[ CCOL_RIV_NUM ]) ? row[ CCOL_RIV_SUM ] : 0.0 ;
 1319 
 1320   /* Return the derived value. */
 1321   return sqrt(1/p->cpscorr) * (I-O) / columns_brightness_error(p, row, o0c1);
 1322 }
 1323 
 1324 
 1325 
 1326 
 1327 
 1328 /* Do the second order calculations, see "Measuring elliptical parameters"
 1329    section of the book/manual for a thorough explanation of the
 1330    derivation. */
 1331 static double
 1332 columns_second_order(struct mkcatalog_passparams *pp, double *row,
 1333                      int key, int o0c1)
 1334 {
 1335   double x=NAN, y=NAN, xx=NAN, yy=NAN, xy=NAN;
 1336   double denom, kx=pp->shift[1], ky=pp->shift[0];
 1337 
 1338   /* Preparations. */
 1339   switch(key)
 1340     {
 1341     /* Brightness weighted. */
 1342     case UI_KEY_SEMIMAJOR:
 1343     case UI_KEY_SEMIMINOR:
 1344     case UI_KEY_POSITIONANGLE:
 1345 
 1346       /* Denominator (to be divided). */
 1347       denom = row[ o0c1 ? CCOL_SUMWHT : OCOL_SUMWHT ];
 1348 
 1349       /* First order. */
 1350       x  = MKC_RATIO( row[ o0c1 ? CCOL_VX : OCOL_VX ], denom );
 1351       y  = MKC_RATIO( row[ o0c1 ? CCOL_VY : OCOL_VY ], denom );
 1352 
 1353       /* Second order. */
 1354       xx = ( MKC_RATIO( row[ o0c1 ? CCOL_VXX : OCOL_VXX ], denom )
 1355              - (x-kx) * (x-kx) );
 1356       yy = ( MKC_RATIO( row[ o0c1 ? CCOL_VYY : OCOL_VYY ], denom )
 1357              - (y-ky) * (y-ky) );
 1358       xy = ( MKC_RATIO( row[ o0c1 ? CCOL_VXY : OCOL_VXY ], denom )
 1359              - (x-kx) * (y-ky) );
 1360       break;
 1361 
 1362     /* Geometric. */
 1363     case UI_KEY_GEOSEMIMAJOR:
 1364     case UI_KEY_GEOSEMIMINOR:
 1365     case UI_KEY_GEOPOSITIONANGLE:
 1366 
 1367       /* Denominator (to be divided). */
 1368       denom = row[ o0c1 ? CCOL_NUMALL : OCOL_NUMALL ];
 1369 
 1370       /* First order. */
 1371       x  = MKC_RATIO( row[ o0c1 ? CCOL_GX : OCOL_GX  ], denom );
 1372       y  = MKC_RATIO( row[ o0c1 ? CCOL_GY : OCOL_GY  ], denom );
 1373 
 1374       /* Second order. */
 1375       xx = ( MKC_RATIO( row[ o0c1 ? CCOL_GXX : OCOL_GXX ], denom )
 1376              - (x-kx) * (x-kx) );
 1377       yy = ( MKC_RATIO( row[ o0c1 ? CCOL_GYY : OCOL_GYY ], denom )
 1378              - (y-ky) * (y-ky) );
 1379       xy = ( MKC_RATIO( row[ o0c1 ? CCOL_GXY : OCOL_GXY ], denom )
 1380              - (x-kx) * (y-ky) );
 1381       break;
 1382 
 1383     /* Error. */
 1384     default:
 1385       error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
 1386             "problem. %d is not a recognized key", __func__, PACKAGE_BUGREPORT,
 1387             key);
 1388     }
 1389 
 1390   /* Return the output. */
 1391   switch(key)
 1392     {
 1393     /* Semi-major axis. */
 1394     case UI_KEY_SEMIMAJOR:
 1395     case UI_KEY_GEOSEMIMAJOR:
 1396       return sqrt( ( xx + yy ) / 2
 1397                    + sqrt( (xx - yy)/2 * (xx - yy)/2 + xy * xy ) );
 1398 
 1399     /* Semi-minor axis. */
 1400     case UI_KEY_SEMIMINOR:
 1401     case UI_KEY_GEOSEMIMINOR:
 1402       return sqrt( ( xx + yy )/2
 1403                    - sqrt( (xx - yy)/2 * (xx - yy)/2 + xy * xy ) );
 1404 
 1405     /* Position angle. */
 1406     case UI_KEY_POSITIONANGLE:
 1407     case UI_KEY_GEOPOSITIONANGLE:
 1408       return 0.5f * atan2(2 * xy, xx - yy) * 180/M_PI;
 1409     }
 1410 
 1411 
 1412   /* Control should not reach here! If it does, its a bug, so abort and let
 1413      the user know. */
 1414   error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s, so we can "
 1415         "address the problem. Control should not have reached the end of "
 1416         "this function", __func__, PACKAGE_BUGREPORT);
 1417   return NAN;
 1418 }
 1419 
 1420 
 1421 
 1422 
 1423 
 1424 /* The clump brightness is needed in several places, so we've defined this
 1425    function to have an easier code. */
 1426 static double
 1427 columns_clump_brightness(double *ci)
 1428 {
 1429   double tmp;
 1430   /* Calculate the river flux over the clump area. But only when rivers are
 1431      present. When grown clumps are requested, the clumps can fully cover a
 1432      detection (that has one or no clumps). */
 1433   tmp = ( ci[ CCOL_RIV_NUM ]>0.0f
 1434           ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]*ci[ CCOL_NUM ]
 1435           : 0 );
 1436 
 1437   /* Subtract it from the clump's brightness. */
 1438   return ci[ CCOL_NUM ]>0.0f ? (ci[ CCOL_SUM ] - tmp) : NAN;
 1439 }
 1440 
 1441 
 1442 
 1443 
 1444 
 1445 /* Measure the minimum and maximum positions. */
 1446 static uint32_t
 1447 columns_xy_extrema(struct mkcatalog_passparams *pp, size_t *coord, int key)
 1448 {
 1449   gal_data_t *tile=pp->tile, *block=tile->block;
 1450 
 1451   /* We only want to do the coordinate estimation once: in `columns_fill',
 1452      we initialized the coordinates with `GAL_BLANK_SIZE_T'. When the
 1453      coordinate has already been measured already, it won't have this value
 1454      any more. */
 1455   if(coord[0]==GAL_BLANK_SIZE_T)
 1456     gal_dimension_index_to_coord(gal_pointer_num_between(block->array,
 1457                                                          tile->array,
 1458                                                          block->type),
 1459                                  block->ndim, block->dsize, coord);
 1460 
 1461   /* Return the proper value: note that `coord' is in C standard: starting
 1462      from the slowest dimension and counting from zero.*/
 1463   switch(key)
 1464     {
 1465     case UI_KEY_MINX: return coord[1] + 1;              break;
 1466     case UI_KEY_MAXX: return coord[1] + tile->dsize[1]; break;
 1467     case UI_KEY_MINY: return coord[0] + 1;              break;
 1468     case UI_KEY_MAXY: return coord[0] + tile->dsize[0]; break;
 1469     default:
 1470       error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
 1471             "problem. The value %d is not a recognized value", __func__,
 1472             PACKAGE_BUGREPORT, key);
 1473     }
 1474 
 1475   /* Control should not reach here. */
 1476   error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix the "
 1477         "problem. Control should not reach the end of this function",
 1478         __func__, PACKAGE_BUGREPORT);
 1479   return GAL_BLANK_UINT32;
 1480 }
 1481 
 1482 
 1483 
 1484 
 1485 
 1486 /* The magnitude error is directly derivable from the S/N:
 1487 
 1488    To derive the error in measuring the magnitude from the S/N, let's take
 1489    `F' as the flux, `Z' is the zeropoint, `M' is the magnitude, `S' is the
 1490    S/N, and `D' to stand for capital delta (or error in a value) then from
 1491 
 1492       `M = -2.5*log10(F) + Z'
 1493 
 1494    we get the following equation after calculating the derivative with
 1495    respect to F.
 1496 
 1497       `dM/df = -2.5 * ( 1 / ( F * ln(10) ) )'
 1498 
 1499    From the Tailor series, `DM' can be written as:
 1500 
 1501       `DM = dM/dF * DF'
 1502 
 1503    So
 1504 
 1505       `DM = |-2.5/ln(10)| * DF/F'
 1506 
 1507    But `DF/F' is just the inverse of the Signal to noise ratio, or
 1508   `1/S'. So
 1509 
 1510       `DM = 2.5 / ( S * ln(10) )'               */
 1511 #define MAG_ERROR(P,ROW,O0C1) ( 2.5f                                    \
 1512                                 / ( ( columns_sn((P),(ROW),(O0C1)) > 0  \
 1513                                       ? columns_sn((P),(ROW),(O0C1))    \
 1514                                       : NAN )                           \
 1515                                     * log(10) ) )
 1516 
 1517 
 1518 
 1519 
 1520 
 1521 
 1522 /* All the raw first and second pass information has been collected, now
 1523    write them into the output columns. The list of columns here is in the
 1524    same order as `columns_alloc_set_out_cols', see there for the type of
 1525    each column. */
 1526 #define POS_V_G(ARRAY, SUMWHT_COL, NUMALL_COL, V_COL, G_COL)            \
 1527   ( (ARRAY)[ SUMWHT_COL ]>0                                             \
 1528     ? MKC_RATIO( (ARRAY)[ V_COL ], (ARRAY)[ SUMWHT_COL ] )              \
 1529     : MKC_RATIO( (ARRAY)[ G_COL ], (ARRAY)[ NUMALL_COL ] ) )
 1530 void
 1531 columns_fill(struct mkcatalog_passparams *pp)
 1532 {
 1533   struct mkcatalogparams *p=pp->p;
 1534 
 1535   int key;
 1536   double tmp;
 1537   void *colarr;
 1538   gal_data_t *column;
 1539   double *ci, *oi=pp->oi;
 1540   size_t coord[2]={GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T};
 1541 
 1542   size_t sr=pp->clumpstartindex, cind, coind;
 1543   size_t oind=pp->object-1; /* IDs start from 1, indexs from 0. */
 1544   double **vo=NULL, **vc=NULL, **go=NULL, **gc=NULL, **vcc=NULL, **gcc=NULL;
 1545 
 1546   /* If a WCS column is requested (check will be done inside the function),
 1547      then set the pointers. */
 1548   columns_set_wcs_pointers(p, &vo, &vc, &go, &gc, &vcc, &gcc);
 1549 
 1550   /* Go over all the object columns and fill in the information. */
 1551   for(column=p->objectcols; column!=NULL; column=column->next)
 1552     {
 1553       /* For easy reading. */
 1554       key=column->status;
 1555       colarr=column->array;
 1556 
 1557       /* Put the number of clumps in the internal array which we will need
 1558          later to order the clump table by object ID. */
 1559       if(p->numclumps_c)
 1560         p->numclumps_c[oind]=pp->clumpsinobj;
 1561 
 1562       /* Go over all the columns. */
 1563       switch(key)
 1564         {
 1565         case UI_KEY_OBJID:
 1566           ((int32_t *)colarr)[oind] = pp->object;
 1567           break;
 1568 
 1569         case UI_KEY_NUMCLUMPS:
 1570           ((int32_t *)colarr)[oind] = pp->clumpsinobj;
 1571           break;
 1572 
 1573         case UI_KEY_AREA:
 1574           ((int32_t *)colarr)[oind] = oi[OCOL_NUM];
 1575           break;
 1576 
 1577         case UI_KEY_CLUMPSAREA:
 1578           ((int32_t *)colarr)[oind] = oi[OCOL_C_NUM];
 1579           break;
 1580 
 1581         case UI_KEY_WEIGHTAREA:
 1582           ((int32_t *)colarr)[oind] = oi[OCOL_NUMWHT];
 1583           break;
 1584 
 1585         case UI_KEY_GEOAREA:
 1586           ((int32_t *)colarr)[oind] = oi[OCOL_NUMALL];
 1587           break;
 1588 
 1589         case UI_KEY_X:
 1590           ((float *)colarr)[oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMWHT,
 1591                                             OCOL_VX, OCOL_GX);
 1592           break;
 1593 
 1594         case UI_KEY_Y:
 1595           ((float *)colarr)[oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMWHT,
 1596                                             OCOL_VY, OCOL_GY);
 1597           break;
 1598 
 1599         case UI_KEY_GEOX:
 1600           ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_GX], oi[OCOL_NUMALL] );
 1601           break;
 1602 
 1603         case UI_KEY_GEOY:
 1604           ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_GY], oi[OCOL_NUMALL] );
 1605           break;
 1606 
 1607         case UI_KEY_CLUMPSX:
 1608           ((float *)colarr)[oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMWHT,
 1609                                             OCOL_C_VX, OCOL_C_GX);
 1610           break;
 1611 
 1612         case UI_KEY_CLUMPSY:
 1613           ((float *)colarr)[oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMWHT,
 1614                                             OCOL_C_VY, OCOL_C_GY);
 1615           break;
 1616 
 1617         case UI_KEY_CLUMPSGEOX:
 1618           ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_C_GX],
 1619                                                oi[OCOL_C_NUMALL] );
 1620           break;
 1621 
 1622         case UI_KEY_CLUMPSGEOY:
 1623           ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_C_GY],
 1624                                                oi[OCOL_C_NUMALL] );
 1625           break;
 1626 
 1627         case UI_KEY_MINX:
 1628         case UI_KEY_MAXX:
 1629         case UI_KEY_MINY:
 1630         case UI_KEY_MAXY:
 1631           ((uint32_t *)colarr)[oind]=columns_xy_extrema(pp, coord, key);
 1632           break;
 1633 
 1634         case UI_KEY_W1:
 1635         case UI_KEY_W2:
 1636           vo[0][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL, OCOL_VX,
 1637                                 OCOL_GX);
 1638           vo[1][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL, OCOL_VY,
 1639                                 OCOL_GY);
 1640           break;
 1641 
 1642         case UI_KEY_GEOW1:
 1643         case UI_KEY_GEOW2:
 1644           go[0][oind] = MKC_RATIO( oi[OCOL_GX], oi[OCOL_NUMALL] );
 1645           go[1][oind] = MKC_RATIO( oi[OCOL_GY], oi[OCOL_NUMALL] );
 1646           break;
 1647 
 1648         case UI_KEY_CLUMPSW1:
 1649         case UI_KEY_CLUMPSW2:
 1650           vcc[0][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL, OCOL_C_VX,
 1651                                  OCOL_C_GX);
 1652           vcc[1][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL, OCOL_C_VY,
 1653                                  OCOL_C_GY);
 1654           break;
 1655 
 1656         case UI_KEY_CLUMPSGEOW1:
 1657         case UI_KEY_CLUMPSGEOW2:
 1658           gcc[0][oind] = MKC_RATIO( oi[OCOL_C_GX], oi[OCOL_C_NUMALL] );
 1659           gcc[1][oind] = MKC_RATIO( oi[OCOL_C_GY], oi[OCOL_C_NUMALL] );
 1660           break;
 1661 
 1662         case UI_KEY_BRIGHTNESS:
 1663           ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
 1664                                       ? oi[ OCOL_SUM ]
 1665                                       : NAN );
 1666           break;
 1667 
 1668         case UI_KEY_BRIGHTNESSERR:
 1669           ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
 1670                                       ? columns_brightness_error(p, oi, 0)
 1671                                       : NAN );
 1672           break;
 1673 
 1674         case UI_KEY_CLUMPSBRIGHTNESS:
 1675           ((float *)colarr)[oind] = ( oi[ OCOL_C_NUM ]>0.0f
 1676                                       ? oi[ OCOL_C_SUM ]
 1677                                       : NAN );
 1678           break;
 1679 
 1680         case UI_KEY_MEAN:
 1681           ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
 1682                                       ? oi[ OCOL_SUM ] / oi[ OCOL_NUM ]
 1683                                       : NAN );
 1684           break;
 1685 
 1686         case UI_KEY_MEDIAN:
 1687           ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
 1688                                       ? oi[ OCOL_MEDIAN ]
 1689                                       : NAN );
 1690           break;
 1691 
 1692         case UI_KEY_MAGNITUDE:
 1693           ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
 1694                                       ? MKC_MAG(oi[ OCOL_SUM ])
 1695                                       : NAN );
 1696           break;
 1697 
 1698         case UI_KEY_MAGNITUDEERR:
 1699           ((float *)colarr)[oind] = MAG_ERROR(p, oi, 0);
 1700           break;
 1701 
 1702         case UI_KEY_CLUMPSMAGNITUDE:
 1703           ((float *)colarr)[oind] = MKC_MAG(oi[ OCOL_C_SUM ]);
 1704           break;
 1705 
 1706         case UI_KEY_UPPERLIMIT:
 1707           ((float *)colarr)[oind] = oi[ OCOL_UPPERLIMIT_B ];
 1708           break;
 1709 
 1710         case UI_KEY_UPPERLIMITMAG:
 1711           ((float *)colarr)[oind] = MKC_MAG(oi[ OCOL_UPPERLIMIT_B ]);
 1712           break;
 1713 
 1714         case UI_KEY_UPPERLIMITONESIGMA:
 1715           ((float *)colarr)[oind] = oi[ OCOL_UPPERLIMIT_S ];
 1716           break;
 1717 
 1718         case UI_KEY_UPPERLIMITSIGMA:
 1719           ((float *)colarr)[oind] = ( ( oi[ OCOL_NUM ]>0.0f
 1720                                         ? oi[ OCOL_SUM ] : NAN )
 1721                                       / oi[ OCOL_UPPERLIMIT_S ] );
 1722           break;
 1723 
 1724         case UI_KEY_UPPERLIMITQUANTILE:
 1725           ((float *)colarr)[oind] = oi[ OCOL_UPPERLIMIT_Q ];
 1726           break;
 1727 
 1728         case UI_KEY_UPPERLIMITSKEW:
 1729           ((float *)colarr)[oind] = oi[ OCOL_UPPERLIMIT_SKEW ];
 1730           break;
 1731 
 1732         case UI_KEY_SN:
 1733           ((float *)colarr)[oind] = columns_sn(p, oi, 0);
 1734           break;
 1735 
 1736         case UI_KEY_SKY:
 1737           ((float *)colarr)[oind] = MKC_RATIO(oi[OCOL_SUMSKY], oi[OCOL_NUM]);
 1738           break;
 1739 
 1740         case UI_KEY_STD:
 1741           ((float *)colarr)[oind] = sqrt( MKC_RATIO( oi[ OCOL_SUMVAR ],
 1742                                                      oi[ OCOL_NUM    ]) );
 1743           break;
 1744 
 1745         case UI_KEY_SEMIMAJOR:
 1746           ((float *)colarr)[oind] = columns_second_order(pp, oi, key, 0);
 1747           break;
 1748 
 1749         case UI_KEY_SEMIMINOR:
 1750           ((float *)colarr)[oind] = columns_second_order(pp, oi, key, 0);
 1751           break;
 1752 
 1753         case UI_KEY_AXISRATIO:
 1754           ((float *)colarr)[oind]
 1755             = ( columns_second_order(pp, oi, UI_KEY_SEMIMINOR, 0)
 1756                 / columns_second_order(pp, oi, UI_KEY_SEMIMAJOR, 0) );
 1757           break;
 1758 
 1759         case UI_KEY_POSITIONANGLE:
 1760           ((float *)colarr)[oind] = columns_second_order(pp, oi, key, 0);
 1761           break;
 1762 
 1763         case UI_KEY_GEOSEMIMAJOR:
 1764           ((float *)colarr)[oind] = columns_second_order(pp, oi, key, 0);
 1765           break;
 1766 
 1767         case UI_KEY_GEOSEMIMINOR:
 1768           ((float *)colarr)[oind] = columns_second_order(pp, oi, key, 0);
 1769           break;
 1770 
 1771         case UI_KEY_GEOAXISRATIO:
 1772           ((float *)colarr)[oind]
 1773             = ( columns_second_order(pp, oi, UI_KEY_GEOSEMIMINOR, 0)
 1774                 / columns_second_order(pp, oi, UI_KEY_GEOSEMIMAJOR, 0) );
 1775           break;
 1776 
 1777         case UI_KEY_GEOPOSITIONANGLE:
 1778           ((float *)colarr)[oind] = columns_second_order(pp, oi, key, 0);
 1779           break;
 1780 
 1781         default:
 1782           error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
 1783                 "solve the problem. the output column code %d not recognized "
 1784                 "(for objects). ", __func__, PACKAGE_BUGREPORT, key);
 1785         }
 1786     }
 1787 
 1788   /* Go over the clump columns and fill the information. */
 1789   for(column=p->clumpcols; column!=NULL; column=column->next)
 1790     for(coind=0;coind<pp->clumpsinobj;++coind)
 1791       {
 1792         /* `coind': clump-in-object-index.
 1793            `cind': clump-index (over all the catalog). */
 1794         cind   = sr + coind;
 1795         colarr = column->array;
 1796         key    = column->status;
 1797         ci     = &pp->ci[ coind * CCOL_NUMCOLS ];
 1798 
 1799         /* Put the object ID of this clump into the temporary array that we
 1800            will later need to sort the final clumps catalog. */
 1801         if(p->hostobjid_c)
 1802           p->hostobjid_c[cind]=pp->object;
 1803 
 1804         /* Parse columns */
 1805         switch(key)
 1806           {
 1807           case UI_KEY_HOSTOBJID:
 1808             ((int32_t *)colarr)[cind]=pp->object;
 1809             break;
 1810 
 1811           case UI_KEY_IDINHOSTOBJ:
 1812             ((int32_t *)colarr)[cind]=coind+1;
 1813             break;
 1814 
 1815           case UI_KEY_AREA:
 1816             ((int32_t *)colarr)[cind]=ci[CCOL_NUM];
 1817             break;
 1818 
 1819           case UI_KEY_WEIGHTAREA:
 1820             ((int32_t *)colarr)[cind]=ci[CCOL_NUMWHT];
 1821             break;
 1822 
 1823           case UI_KEY_GEOAREA:
 1824             ((int32_t *)colarr)[cind]=ci[CCOL_NUMALL];
 1825             break;
 1826 
 1827           case UI_KEY_X:
 1828             ((float *)colarr)[cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
 1829                                               CCOL_VX, CCOL_GX);
 1830             break;
 1831 
 1832           case UI_KEY_Y:
 1833             ((float *)colarr)[cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
 1834                                               CCOL_VY, CCOL_GY);
 1835             break;
 1836 
 1837           case UI_KEY_GEOX:
 1838             ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_GX],
 1839                                                  ci[CCOL_NUMALL] );
 1840             break;
 1841 
 1842           case UI_KEY_GEOY:
 1843             ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_GY],
 1844                                                  ci[CCOL_NUMALL] );
 1845             break;
 1846 
 1847          case UI_KEY_MINX:
 1848             ((uint32_t *)colarr)[cind] = ci[CCOL_MINX];
 1849             break;
 1850 
 1851          case UI_KEY_MAXX:
 1852             ((uint32_t *)colarr)[cind] = ci[CCOL_MAXX];
 1853             break;
 1854 
 1855          case UI_KEY_MINY:
 1856             ((uint32_t *)colarr)[cind] = ci[CCOL_MINY];
 1857             break;
 1858 
 1859          case UI_KEY_MAXY:
 1860             ((uint32_t *)colarr)[cind] = ci[CCOL_MAXY];
 1861             break;
 1862 
 1863           case UI_KEY_W1:
 1864           case UI_KEY_W2:
 1865             vc[0][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL, CCOL_VX,
 1866                                   CCOL_GX);
 1867             vc[1][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL, CCOL_VY,
 1868                                   CCOL_GY);
 1869             break;
 1870 
 1871           case UI_KEY_GEOW1:
 1872           case UI_KEY_GEOW2:
 1873             gc[0][cind] = MKC_RATIO( ci[CCOL_GX], ci[CCOL_NUMALL] );
 1874             gc[1][cind] = MKC_RATIO( ci[CCOL_GY], ci[CCOL_NUMALL] );
 1875             break;
 1876 
 1877           case UI_KEY_BRIGHTNESS:
 1878             ((float *)colarr)[cind] = columns_clump_brightness(ci);
 1879             break;
 1880 
 1881           case UI_KEY_BRIGHTNESSERR:
 1882             ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
 1883                                         ? columns_brightness_error(p, ci, 1)
 1884                                         : NAN );
 1885             break;
 1886 
 1887           case UI_KEY_BRIGHTNESSNORIVER:
 1888             ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
 1889                                         ? ci[ CCOL_SUM ] : NAN );
 1890             break;
 1891 
 1892           case UI_KEY_MEAN:
 1893             ((float *)colarr)[cind] = ( columns_clump_brightness(ci)
 1894                                         /ci[CCOL_NUM] );
 1895             break;
 1896 
 1897           case UI_KEY_MEDIAN:
 1898             ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
 1899                                         ? ci[ CCOL_MEDIAN ] : NAN );
 1900             break;
 1901 
 1902           case UI_KEY_MAGNITUDE: /* Similar: brightness for clumps */
 1903             tmp = ( ci[ CCOL_RIV_NUM ]
 1904                     ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]*ci[ CCOL_NUM ]
 1905                     : 0 );
 1906             ((float *)colarr)[cind] = MKC_MAG(ci[ CCOL_SUM ]-tmp);
 1907             break;
 1908 
 1909           case UI_KEY_MAGNITUDEERR:
 1910             ((float *)colarr)[cind] = MAG_ERROR(p, ci, 1);
 1911             break;
 1912 
 1913           case UI_KEY_UPPERLIMIT:
 1914             ((float *)colarr)[cind] = ci[ CCOL_UPPERLIMIT_B ];
 1915             break;
 1916 
 1917           case UI_KEY_UPPERLIMITMAG:
 1918             ((float *)colarr)[cind] = MKC_MAG(ci[ CCOL_UPPERLIMIT_B ]);
 1919             break;
 1920 
 1921           case UI_KEY_UPPERLIMITONESIGMA:
 1922             ((float *)colarr)[cind] = ci[ CCOL_UPPERLIMIT_S ];
 1923             break;
 1924 
 1925           case UI_KEY_UPPERLIMITSIGMA:
 1926             ((float *)colarr)[cind] = ( columns_clump_brightness(ci)
 1927                                         / ci[ CCOL_UPPERLIMIT_S ] );
 1928             break;
 1929 
 1930           case UI_KEY_UPPERLIMITQUANTILE:
 1931             ((float *)colarr)[cind] = ci[ CCOL_UPPERLIMIT_Q ];
 1932             break;
 1933 
 1934           case UI_KEY_UPPERLIMITSKEW:
 1935             ((float *)colarr)[cind] = ci[ CCOL_UPPERLIMIT_SKEW ];
 1936             break;
 1937 
 1938           case UI_KEY_RIVERAVE:
 1939             ((float *)colarr)[cind] = ( ci[ CCOL_RIV_NUM]
 1940                                         ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM]
 1941                                         : NAN );
 1942             break;
 1943 
 1944           case UI_KEY_RIVERNUM:
 1945             ((int32_t *)colarr)[cind] = ci[ CCOL_RIV_NUM ];
 1946             break;
 1947 
 1948           case UI_KEY_SN:
 1949             ((float *)colarr)[cind] = columns_sn(p, ci, 1);
 1950             break;
 1951 
 1952           case UI_KEY_SKY:
 1953             ((float *)colarr)[cind] = MKC_RATIO( ci[ CCOL_SUMSKY],
 1954                                                  ci[ CCOL_NUM] );
 1955             break;
 1956 
 1957           case UI_KEY_STD:
 1958             ((float *)colarr)[cind] = sqrt( MKC_RATIO( ci[ CCOL_SUMVAR ],
 1959                                                        ci[ CCOL_NUM    ] ));
 1960             break;
 1961 
 1962           case UI_KEY_SEMIMAJOR:
 1963             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
 1964             break;
 1965 
 1966           case UI_KEY_SEMIMINOR:
 1967             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
 1968             break;
 1969 
 1970           case UI_KEY_AXISRATIO:
 1971             ((float *)colarr)[cind]
 1972               = ( columns_second_order(  pp, ci, UI_KEY_SEMIMINOR, 1)
 1973                   / columns_second_order(pp, ci, UI_KEY_SEMIMAJOR, 1) );
 1974             break;
 1975 
 1976           case UI_KEY_POSITIONANGLE:
 1977             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
 1978             break;
 1979 
 1980           case UI_KEY_GEOSEMIMAJOR:
 1981             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
 1982             break;
 1983 
 1984           case UI_KEY_GEOSEMIMINOR:
 1985             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
 1986             break;
 1987 
 1988           case UI_KEY_GEOAXISRATIO:
 1989             ((float *)colarr)[cind]
 1990               = ( columns_second_order(  pp, ci, UI_KEY_GEOSEMIMINOR, 1)
 1991                   / columns_second_order(pp, ci, UI_KEY_GEOSEMIMAJOR, 1) );
 1992             break;
 1993 
 1994           case UI_KEY_GEOPOSITIONANGLE:
 1995             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
 1996             break;
 1997 
 1998           default:
 1999             error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
 2000                   "solve the problem. The output column code %d not "
 2001                   "recognized (for clumps). ", __func__, PACKAGE_BUGREPORT,
 2002                   key);
 2003           }
 2004       }
 2005 
 2006   /* Clean up. */
 2007   if(vo)  free(vo);
 2008   if(vc)  free(vc);
 2009   if(go)  free(go);
 2010   if(gc)  free(gc);
 2011   if(vcc) free(vcc);
 2012   if(gcc) free(gcc);
 2013 }