"Fossies" - the Fresh Open Source Software Archive

Member "gnuastro-0.9/bin/mkcatalog/mkcatalog.c" (11 Apr 2019, 22166 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 "mkcatalog.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) 2015-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 <float.h>
   30 #include <string.h>
   31 #include <stdlib.h>
   32 #include <inttypes.h>
   33 
   34 #include <gnuastro/git.h>
   35 #include <gnuastro/wcs.h>
   36 #include <gnuastro/data.h>
   37 #include <gnuastro/fits.h>
   38 #include <gnuastro/threads.h>
   39 #include <gnuastro/pointer.h>
   40 #include <gnuastro/dimension.h>
   41 #include <gnuastro/statistics.h>
   42 #include <gnuastro/permutation.h>
   43 
   44 #include <gnuastro-internal/timing.h>
   45 
   46 #include "main.h"
   47 #include "mkcatalog.h"
   48 
   49 #include "ui.h"
   50 #include "parse.h"
   51 #include "columns.h"
   52 #include "upperlimit.h"
   53 
   54 
   55 
   56 
   57 
   58 
   59 
   60 
   61 
   62 
   63 /*********************************************************************/
   64 /**************       Manage a single object       *******************/
   65 /*********************************************************************/
   66 static void
   67 mkcatalog_clump_starting_index(struct mkcatalog_passparams *pp)
   68 {
   69   struct mkcatalogparams *p=pp->p;
   70 
   71   /* Lock the mutex if we are working on more than one thread. NOTE: it is
   72      very important to keep the number of operations within the mutex to a
   73      minimum so other threads don't get delayed. */
   74   if(p->cp.numthreads>1)
   75     pthread_mutex_lock(&p->mutex);
   76 
   77   /* Put the current total number of rows filled into the output, then
   78      increment the total number by the number of clumps. */
   79   pp->clumpstartindex = p->clumprowsfilled;
   80   p->clumprowsfilled += pp->clumpsinobj;
   81 
   82   /* Unlock the mutex (if it was locked). */
   83   if(p->cp.numthreads>1)
   84     pthread_mutex_unlock(&p->mutex);
   85 }
   86 
   87 
   88 
   89 
   90 
   91 /* Each thread will call this function once. It will go over all the
   92    objects that are assigned to it. */
   93 static void *
   94 mkcatalog_single_object(void *in_prm)
   95 {
   96   struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
   97   struct mkcatalogparams *p=(struct mkcatalogparams *)(tprm->params);
   98   size_t ndim=p->objects->ndim;
   99 
  100   size_t i;
  101   uint8_t *oif=p->oiflag;
  102   struct mkcatalog_passparams pp;
  103 
  104 
  105   /* Initialize the mkcatalog_passparams elements. */
  106   pp.p               = p;
  107   pp.clumpstartindex = 0;
  108   pp.rng             = p->rng ? gsl_rng_clone(p->rng) : NULL;
  109   pp.oi              = gal_pointer_allocate(GAL_TYPE_FLOAT64, OCOL_NUMCOLS,
  110                                             0, __func__, "pp.oi");
  111 
  112   /* If we have second order measurements, allocate the array keeping the
  113      temporary shift values for each object of this thread. Note that the
  114      clumps catalog (if requested), will have the same measurements, so its
  115      just enough to check the objects. */
  116   pp.shift = ( ( oif[    OCOL_GXX ]
  117                  || oif[ OCOL_GYY ]
  118                  || oif[ OCOL_GXY ]
  119                  || oif[ OCOL_VXX ]
  120                  || oif[ OCOL_VYY ]
  121                  || oif[ OCOL_VXY ] )
  122                ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
  123                                        "pp.shift")
  124                : NULL );
  125 
  126   /* If we have upper-limit mode, then allocate the container to keep the
  127      values to calculate the standard deviation. */
  128   if(p->upperlimit)
  129     {
  130       /* Allocate the space to keep the upper-limit values. */
  131       pp.up_vals = gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &p->upnum,
  132                                   NULL, 0, p->cp.minmapsize, NULL, NULL,
  133                                   NULL);
  134 
  135       /* Set the blank checked flag to 1. By definition, this dataset won't
  136          have any blank values. Also `flag' is initialized to `0'. So we
  137          just have to set the checked flag (`GAL_DATA_FLAG_BLANK_CH') to
  138          one to inform later steps that there are no blank values. */
  139       pp.up_vals->flag |= GAL_DATA_FLAG_BLANK_CH;
  140     }
  141   else
  142     pp.up_vals=NULL;
  143 
  144 
  145   /* Fill the desired columns for all the objects given to this thread. */
  146   for(i=0; tprm->indexs[i]!=GAL_BLANK_SIZE_T; ++i)
  147     {
  148       /* For easy reading. Note that the object IDs start from one while
  149          the array positions start from 0. */
  150       pp.ci     = NULL;
  151       pp.object = tprm->indexs[i] + 1;
  152       pp.tile   = &p->tiles[ tprm->indexs[i] ];
  153 
  154       /* Initialize the parameters for this object/tile. */
  155       parse_initialize(&pp);
  156 
  157       /* Get the first pass information. */
  158       parse_objects(&pp);
  159 
  160       /* Currently the second pass is only necessary when there is a clumps
  161          image. */
  162       if(p->clumps)
  163         {
  164           /* Allocate space for the properties of each clump. */
  165           pp.ci = gal_pointer_allocate(GAL_TYPE_FLOAT64,
  166                                        pp.clumpsinobj * CCOL_NUMCOLS, 1,
  167                                        __func__, "pp.ci");
  168 
  169           /* Get the starting row of this object's clumps in the final
  170              catalog. This index is also necessary for the unique random
  171              number generator seeds of each clump. */
  172           mkcatalog_clump_starting_index(&pp);
  173 
  174           /* Get the second pass information. */
  175           parse_clumps(&pp);
  176         }
  177 
  178       /* If the median is requested, another pass is necessary. */
  179       if( p->oiflag[ OCOL_MEDIAN ] )
  180         parse_median(&pp);
  181 
  182       /* Calculate the upper limit magnitude (if necessary). */
  183       if(p->upperlimit) upperlimit_calculate(&pp);
  184 
  185       /* Write the pass information into the columns. */
  186       columns_fill(&pp);
  187 
  188       /* Clean up for this object. */
  189       if(pp.ci) free(pp.ci);
  190     }
  191 
  192   /* Clean up. */
  193   free(pp.oi);
  194   free(pp.shift);
  195   gal_data_free(pp.up_vals);
  196   if(pp.rng) gsl_rng_free(pp.rng);
  197 
  198   /* Wait until all the threads finish and return. */
  199   if(tprm->b) pthread_barrier_wait(tprm->b);
  200   return NULL;
  201 }
  202 
  203 
  204 
  205 
  206 
  207 
  208 
  209 
  210 
  211 
  212 
  213 
  214 
  215 
  216 
  217 
  218 
  219 
  220 
  221 
  222 /*********************************************************************/
  223 /********         Processing after threads finish        *************/
  224 /*********************************************************************/
  225 /* Convert internal image coordinates to WCS for table.
  226 
  227    Note that from the beginning (during the passing steps), we saved FITS
  228    coordinates. Also note that we are doing the conversion in place. */
  229 static void
  230 mkcatalog_wcs_conversion(struct mkcatalogparams *p)
  231 {
  232   gal_data_t *c;
  233   gal_data_t *column;
  234 
  235   /* Flux weighted center positions for clumps and objects. */
  236   if(p->wcs_vo)
  237     {
  238       gal_wcs_img_to_world(p->wcs_vo, p->objects->wcs, 1);
  239       if(p->wcs_vc)
  240         gal_wcs_img_to_world(p->wcs_vc, p->objects->wcs, 1);
  241     }
  242 
  243 
  244   /* Geometric center positions for clumps and objects. */
  245   if(p->wcs_go)
  246     {
  247       gal_wcs_img_to_world(p->wcs_go, p->objects->wcs, 1);
  248       if(p->wcs_gc)
  249         gal_wcs_img_to_world(p->wcs_gc, p->objects->wcs, 1);
  250     }
  251 
  252 
  253   /* All clumps flux weighted center. */
  254   if(p->wcs_vcc)
  255     gal_wcs_img_to_world(p->wcs_vcc, p->objects->wcs, 1);
  256 
  257 
  258   /* All clumps geometric center. */
  259   if(p->wcs_gcc)
  260     gal_wcs_img_to_world(p->wcs_gcc, p->objects->wcs, 1);
  261 
  262 
  263   /* Go over all the object columns and fill in the values. */
  264   for(column=p->objectcols; column!=NULL; column=column->next)
  265     {
  266       /* Definitions */
  267       c=NULL;
  268 
  269       /* Set `c' for the columns that must be corrected. Note that this
  270          `switch' statement doesn't need any `default', because there are
  271          probably columns that don't need any correction. */
  272       switch(column->status)
  273         {
  274         case UI_KEY_W1:           c=p->wcs_vo;                break;
  275         case UI_KEY_W2:           c=p->wcs_vo->next;          break;
  276         case UI_KEY_GEOW1:        c=p->wcs_go;                break;
  277         case UI_KEY_GEOW2:        c=p->wcs_go->next;          break;
  278         case UI_KEY_CLUMPSW1:     c=p->wcs_vcc;               break;
  279         case UI_KEY_CLUMPSW2:     c=p->wcs_vcc->next;         break;
  280         case UI_KEY_CLUMPSGEOW1:  c=p->wcs_gcc;               break;
  281         case UI_KEY_CLUMPSGEOW2:  c=p->wcs_gcc->next;         break;
  282         }
  283 
  284       /* Copy the elements into the output column. */
  285       if(c)
  286         memcpy(column->array, c->array,
  287                column->size*gal_type_sizeof(c->type));
  288     }
  289 
  290 
  291   /* Go over all the clump columns and fill in the values. */
  292   for(column=p->clumpcols; column!=NULL; column=column->next)
  293     {
  294       /* Definitions */
  295       c=NULL;
  296 
  297       /* Set `c' for the columns that must be corrected. Note that this
  298          `switch' statement doesn't need any `default', because there are
  299          probably columns that don't need any correction. */
  300       switch(column->status)
  301         {
  302         case UI_KEY_W1:           c=p->wcs_vc;                break;
  303         case UI_KEY_W2:           c=p->wcs_vc->next;          break;
  304         case UI_KEY_GEOW1:        c=p->wcs_gc;                break;
  305         case UI_KEY_GEOW2:        c=p->wcs_gc->next;          break;
  306         }
  307 
  308       /* Copy the elements into the output column. */
  309       if(c)
  310         memcpy(column->array, c->array,
  311                column->size*gal_type_sizeof(c->type));
  312     }
  313 }
  314 
  315 
  316 
  317 
  318 
  319 void
  320 mkcatalog_write_inputs_in_comments(struct mkcatalogparams *p,
  321                                    gal_list_str_t **comments, int withsky,
  322                                    int withstd)
  323 {
  324   char *tmp, *str;
  325 
  326   /* Basic classifiers for plain text outputs. */
  327   if(p->cp.tableformat==GAL_TABLE_FORMAT_TXT)
  328     {
  329       if( asprintf(&str, "--------- Input files ---------")<0 )
  330         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  331       gal_list_str_add(comments, str, 0);
  332     }
  333 
  334   /* Object labels. */
  335   if( asprintf(&str, "Objects: %s (hdu: %s).", p->objectsfile, p->cp.hdu)<0 )
  336     error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  337   gal_list_str_add(comments, str, 0);
  338 
  339   /* Clump labels. */
  340   if(p->clumps)
  341     {
  342       if(asprintf(&str, "Clumps:  %s (hdu: %s).", p->usedclumpsfile,
  343                   p->clumpshdu)<0)
  344         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  345       gal_list_str_add(comments, str, 0);
  346     }
  347 
  348   /* Values dataset. */
  349   if(p->values)
  350     {
  351       if( asprintf(&str, "Values:  %s (hdu: %s).", p->usedvaluesfile,
  352                    p->valueshdu)<0 )
  353         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  354       gal_list_str_add(comments, str, 0);
  355     }
  356 
  357   /* Sky dataset. */
  358   if(withsky && p->sky)
  359     {
  360       if(p->sky->size==1)
  361         {
  362           if( asprintf(&str, "Sky:     %g.", *((float *)(p->sky->array)) )<0 )
  363             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  364         }
  365       else
  366         {
  367           if( asprintf(&str, "Sky:     %s (hdu: %s).", p->usedskyfile,
  368                        p->skyhdu)<0 )
  369             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  370         }
  371       gal_list_str_add(comments, str, 0);
  372     }
  373 
  374   /* Sky standard deviation dataset. */
  375   tmp = p->variance ? "VAR" : "STD";
  376   if(withstd && p->std)
  377     {
  378       if(p->std->size==1)
  379         {
  380           if( asprintf(&str, "Sky %s: %g.", tmp,
  381                        *((float *)(p->std->array)) )<0 )
  382             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  383         }
  384       else
  385         {
  386           if( asprintf(&str, "Sky %s: %s (hdu: %s).", tmp, p->usedstdfile,
  387                        p->stdhdu)<0 )
  388             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  389         }
  390       gal_list_str_add(comments, str, 0);
  391     }
  392 
  393   /* Upper limit mask. */
  394   if(p->upmaskfile)
  395     {
  396       if( asprintf(&str, "Upperlimit mask: %s (hdu: %s).", p->upmaskfile,
  397                    p->upmaskhdu)<0 )
  398         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  399       gal_list_str_add(comments, str, 0);
  400     }
  401 }
  402 
  403 
  404 
  405 
  406 
  407 /* Write the similar information. */
  408 static gal_list_str_t *
  409 mkcatalog_outputs_same_start(struct mkcatalogparams *p, int o0c1,
  410                              char *ObjClump)
  411 {
  412   char *str, *tstr;
  413   double pixarea=NAN;
  414   gal_list_str_t *comments=NULL;
  415 
  416   if( asprintf(&str, "%s catalog of %s", o0c1 ? "Object" : "Clump",
  417                PROGRAM_STRING)<0 )
  418     error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  419   gal_list_str_add(&comments, str, 0);
  420 
  421   /* If in a Git controlled directory and output isn't a FITS file (in
  422      FITS, this will be automatically included). */
  423   if(p->cp.tableformat==GAL_TABLE_FORMAT_TXT && gal_git_describe())
  424     {
  425       if(asprintf(&str, "Working directory commit %s", gal_git_describe())<0)
  426         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  427       gal_list_str_add(&comments, str, 0);
  428     }
  429 
  430   /* Write the date. However, `ctime' is going to put a new-line character
  431      in the end of its string, so we are going to remove it manually. */
  432   if( asprintf(&str, "%s started on %s", PROGRAM_NAME, ctime(&p->rawtime))<0 )
  433     error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  434   str[strlen(str)-1]='\0';
  435   gal_list_str_add(&comments, str, 0);
  436 
  437 
  438   /* Write the basic information. */
  439   mkcatalog_write_inputs_in_comments(p, &comments, 1, 1);
  440 
  441 
  442   /* Write other supplimentary information. */
  443   if(p->cp.tableformat==GAL_TABLE_FORMAT_TXT)
  444     {
  445       if( asprintf(&str, "--------- Supplimentary information ---------")<0 )
  446         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  447       gal_list_str_add(&comments, str, 0);
  448     }
  449 
  450   if(p->objects->wcs)
  451     {
  452       pixarea=gal_wcs_pixel_area_arcsec2(p->objects->wcs);
  453       if( isnan(pixarea)==0 )
  454         {
  455           if( asprintf(&str, "Pixel area (arcsec^2): %g", pixarea)<0 )
  456             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  457           gal_list_str_add(&comments, str, 0);
  458         }
  459     }
  460 
  461   if(p->hasmag)
  462     {
  463       if( asprintf(&str, "Zeropoint magnitude: %.4f", p->zeropoint)<0 )
  464         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  465       gal_list_str_add(&comments, str, 0);
  466     }
  467 
  468   /* Print surface brightness limits. */
  469   if( !isnan(p->medstd) && !isnan(p->zeropoint) &&  !isnan(p->sfmagnsigma) )
  470     {
  471       /* Per pixel. */
  472       if( asprintf(&str, "%g sigma surface brightness (magnitude/pixel): "
  473                    "%.3f", p->sfmagnsigma, ( -2.5f
  474                                              *log10( p->sfmagnsigma
  475                                                      * p->medstd )
  476                                              + p->zeropoint ) )<0 )
  477         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  478       gal_list_str_add(&comments, str, 0);
  479 
  480       /* Requested projected area: if a pixel area could be measured (a WCS
  481          was given), then also estimate the surface brightness over one
  482          arcsecond^2. From the pixel area, we know how many pixels are
  483          necessary to fill the requested projected area (in
  484          arcsecond^2). We also know that as the number of samples (pixels)
  485          increases (to N), the noise increases by sqrt(N), see the full
  486          discussion in the book. */
  487       if(!isnan(pixarea) && !isnan(p->sfmagarea))
  488         {
  489           /* Prepare the comment/information. */
  490           if(p->sfmagarea==1.0f)
  491             tstr=NULL;
  492           else
  493             if( asprintf(&tstr, "%g-", p->sfmagarea)<0 )
  494               error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  495           if( asprintf(&str, "%g sigma surface brightness "
  496                        "(magnitude/%sarcsec^2): %.3f", p->sfmagnsigma,
  497                        tstr ? tstr : "",
  498                        ( -2.5f * log10( p->sfmagnsigma
  499                                         * p->medstd
  500                                         * sqrt( p->sfmagarea / pixarea) )
  501                          + p->zeropoint ) )<0 )
  502             error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  503 
  504           /* Add the final string/line to the catalog comments. */
  505           gal_list_str_add(&comments, str, 0);
  506 
  507           /* Clean up (if necessary). */
  508           if (tstr)
  509             {
  510               free(tstr);
  511               tstr=NULL;
  512             }
  513         }
  514 
  515       /* Notice: */
  516       if( asprintf(&str, "Pixel STD for surface brightness calculation%s: %f",
  517                    (!isnan(pixarea) && !isnan(p->sfmagarea))?"s":"",
  518                    p->medstd)<0 )
  519         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  520       gal_list_str_add(&comments, str, 0);
  521     }
  522 
  523   if(p->cpscorr>1.0f)
  524     {
  525       if( asprintf(&str, "Counts-per-second correction: %.3f", p->cpscorr)<0 )
  526         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  527       gal_list_str_add(&comments, str, 0);
  528     }
  529 
  530   if(p->upperlimit)
  531     upperlimit_write_comments(p, &comments, 1);
  532 
  533 
  534 
  535   if(p->cp.tableformat==GAL_TABLE_FORMAT_TXT)
  536     {
  537       if( asprintf(&str, "--------- Table columns ---------")<0 )
  538         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  539       gal_list_str_add(&comments, str, 0);
  540     }
  541 
  542   /* Return the comments. */
  543   return comments;
  544 }
  545 
  546 
  547 
  548 
  549 
  550 
  551 /* Since all the measurements were done in parallel (and we didn't know the
  552    number of clumps per object a-priori), the clumps informtion is just
  553    written in as they are measured. Here, we'll sort the clump columns by
  554    object ID. There is an option to disable this. */
  555 static void
  556 sort_clumps_by_objid(struct mkcatalogparams *p)
  557 {
  558   gal_data_t *col;
  559   size_t o, i, j, *permute, *rowstart;
  560 
  561   /* Make sure everything is fine. */
  562   if(p->hostobjid_c==NULL || p->numclumps_c==NULL)
  563     error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
  564           "problem. `p->hostobjid_c' and `p->numclumps_c' must not be "
  565           "NULL.", __func__, PACKAGE_BUGREPORT);
  566 
  567 
  568   /* Allocate the necessary arrays. */
  569   rowstart=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->numobjects, 0, __func__,
  570                                  "rowstart");
  571   permute=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->numclumps, 0, __func__,
  572                                 "permute");
  573 
  574 
  575   /* The objects array is already sorted by object ID. So we should just
  576      add up the number of clumps to find the row where each object's clumps
  577      should start from in the final sorted clumps catalog. */
  578   rowstart[0] = 0;
  579   for(i=1;i<p->numobjects;++i)
  580     rowstart[i] = p->numclumps_c[i-1] + rowstart[i-1];
  581 
  582   /* Fill the permutation array. Note that WE KNOW that all the objects for
  583      one clump are after each other.*/
  584   i=0;
  585   while(i<p->numclumps)
  586     {
  587       o=p->hostobjid_c[i]-1;
  588       for(j=0; j<p->numclumps_c[o]; ++j)
  589         permute[i++] = rowstart[o] + j;
  590     }
  591 
  592   /* Permute all the clump columns. */
  593   for(col=p->clumpcols; col!=NULL; col=col->next)
  594     gal_permutation_apply_inverse(col, permute);
  595 
  596   /* Clean up */
  597   free(permute);
  598   free(rowstart);
  599 }
  600 
  601 
  602 
  603 
  604 
  605 /* Write the produced columns into the output */
  606 static void
  607 mkcatalog_write_outputs(struct mkcatalogparams *p)
  608 {
  609   /*char *str;*/
  610   gal_list_str_t *comments;
  611 
  612 
  613   /* OBJECT catalog */
  614   comments=mkcatalog_outputs_same_start(p, 0, "Detection");
  615 
  616   /* Reverse the comments list (so it is printed in the same order here),
  617      write the objects catalog and free the comments. */
  618   gal_list_str_reverse(&comments);
  619   gal_table_write(p->objectcols, comments, p->cp.tableformat,
  620                   p->objectsout, "OBJECTS", 0);
  621   gal_list_str_free(comments, 1);
  622 
  623 
  624   /* CLUMPS catalog */
  625   if(p->clumps)
  626     {
  627       /* Make the comments. */
  628       comments=mkcatalog_outputs_same_start(p, 1, "Clumps");
  629 
  630       /* Reverse the comments and write the catalog. */
  631       gal_list_str_reverse(&comments);
  632       gal_table_write(p->clumpcols, comments, p->cp.tableformat,
  633                       p->clumpsout, "CLUMPS", 0);
  634       gal_list_str_free(comments, 1);
  635     }
  636 
  637 
  638   /* Configuration information. */
  639   if(gal_fits_name_is_fits(p->objectsout))
  640     {
  641       gal_fits_key_write_filename("input", p->objectsfile, &p->cp.okeys, 1);
  642       gal_fits_key_write_config(&p->cp.okeys, "MakeCatalog configuration",
  643                                 "MKCATALOG-CONFIG", p->objectsout, "0");
  644     }
  645 
  646 
  647   /* Inform the user */
  648   if(!p->cp.quiet)
  649     {
  650       if(p->clumpsout && strcmp(p->clumpsout,p->objectsout))
  651         {
  652           printf("  - Output objects catalog: %s\n", p->objectsout);
  653           if(p->clumps)
  654             printf("  - Output clumps catalog: %s\n", p->clumpsout);
  655         }
  656       else
  657         printf("  - Catalog written to %s\n", p->objectsout);
  658     }
  659 }
  660 
  661 
  662 
  663 
  664 
  665 
  666 
  667 
  668 
  669 
  670 
  671 
  672 
  673 
  674 
  675 
  676 
  677 
  678 
  679 
  680 /*********************************************************************/
  681 /*****************       Top-level function        *******************/
  682 /*********************************************************************/
  683 void
  684 mkcatalog(struct mkcatalogparams *p)
  685 {
  686   /* When more than one thread is to be used, initialize the mutex: we need
  687      it to assign a column to the clumps in the final catalog. */
  688   if( p->cp.numthreads > 1 ) pthread_mutex_init(&p->mutex, NULL);
  689 
  690   /* Do the processing on each thread. */
  691   gal_threads_spin_off(mkcatalog_single_object, p, p->numobjects,
  692                        p->cp.numthreads);
  693 
  694   /* Post-thread processing, for example to convert image coordinates to RA
  695      and Dec. */
  696   mkcatalog_wcs_conversion(p);
  697 
  698   /* If the columns need to be sorted (by object ID), then some adjustments
  699      need to be made (possibly to both the objects and clumps catalogs). */
  700   if(p->hostobjid_c)
  701     sort_clumps_by_objid(p);
  702 
  703   /* Write the filled columns into the output. */
  704   mkcatalog_write_outputs(p);
  705 
  706   /* Destroy the mutex. */
  707   if( p->cp.numthreads>1 ) pthread_mutex_destroy(&p->mutex);
  708 }