[gnuastro-commits] master f96bbdd: Statistics: --contour option calculates contours for PGFPlots

[Mohammad Akhlaghi]

branch: master
commit f96bbdd502a8fce2673207deda17294f9636253a
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Statistics: --contour option calculates contours for PGFPlots
    
    One of the common things we need when doing an 2D analysis of an image is a
    contour plot. Plotting tools like Matplotlib do calculate it internally,
    but others like the LaTeX PGFPlots package doesn't. PGFPlots doesn't need
    all the dependencies of Matplotlib and is much more easier to build within
    LaTeX and creates much higher quality plots.
    
    With this option, it is now very easy to create the necessary inputs to
    PGFPlots and use it to make contours over a dispalyed image.
    
    To do this, a new function has also been added to the library:
    `gal_binary_connected_indexs', it operates very similar to the connected
    component labeling function, but just keeps the labels.
---
 NEWS                        |   8 ++
 bin/statistics/Makefile.am  |   4 +-
 bin/statistics/args.h       |  15 ++++
 bin/statistics/contour.c    | 209 ++++++++++++++++++++++++++++++++++++++++++++
 bin/statistics/contour.h    |  31 +++++++
 bin/statistics/main.h       |   1 +
 bin/statistics/statistics.c |   7 ++
 bin/statistics/ui.c         |   2 +-
 bin/statistics/ui.h         |   3 +-
 doc/gnuastro.texi           |  19 ++++
 lib/binary.c                |  91 +++++++++++++++++++
 lib/gnuastro/binary.h       |   3 +
 12 files changed, 389 insertions(+), 4 deletions(-)

diff --git a/NEWS b/NEWS
index 0aa6a97..40b9cf8 100644
--- a/NEWS
+++ b/NEWS
@@ -25,6 +25,11 @@ See the end of the file for license conditions.
      the redshift given to CosmicCalculator. You can either use known line
      names, or directly give a number as any emitted line's wavelength.
 
+  Statistics:
+   --contour: compute a contour plot which can be directly fed into the
+     PGFPlots package of LaTeX for plotting the contours. Support for more
+     formats will be added based on the need/request.
+
   Table:
    --equal: Output only rows that have a value equal to the given value in
      the given column. For example `--equal=ID,2,4,5' will select only rows
@@ -37,6 +42,9 @@ See the end of the file for license conditions.
        columns, or the distances of all the points in the table rows with a
        fixed point. See the book for examples and better explanation.
 
+  Library:
+   - gal_binary_connected_indexs: store indexs of connected components.
+
 ** Removed features
 
 ** Changed features
diff --git a/bin/statistics/Makefile.am b/bin/statistics/Makefile.am
index 4106976..540757c 100644
--- a/bin/statistics/Makefile.am
+++ b/bin/statistics/Makefile.am
@@ -35,9 +35,9 @@ bin_PROGRAMS = aststatistics
 aststatistics_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la \
                       -lgnuastro $(MAYBE_NORPATH)
 
-aststatistics_SOURCES = main.c ui.c sky.c statistics.c
+aststatistics_SOURCES = main.c ui.c contour.c sky.c statistics.c
 
-EXTRA_DIST = main.h authors-cite.h args.h ui.h sky.h statistics.h
+EXTRA_DIST = main.h authors-cite.h args.h ui.h sky.h statistics.h contour.h
 
 
 
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 3eaed38..9ee8852 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -475,6 +475,21 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "contour",
+      UI_KEY_CONTOUR,
+      "STR",
+      0,
+      "Contour levels, save in PGFPlots format.",
+      UI_GROUP_PARTICULAR_STAT,
+      &p->contour,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      gal_options_parse_csv_float64
+    },
+
 
 
 
diff --git a/bin/statistics/contour.c b/bin/statistics/contour.c
new file mode 100644
index 0000000..f6e90f8
--- /dev/null
+++ b/bin/statistics/contour.c
@@ -0,0 +1,209 @@
+/*********************************************************************
+Statistics - Statistical analysis on input dataset.
+Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2019, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <stdlib.h>
+
+#include <gnuastro-internal/checkset.h>
+
+#include <gnuastro/wcs.h>
+#include <gnuastro/binary.h>
+#include <gnuastro/arithmetic.h>
+
+#include "main.h"
+#include "contour.h"
+
+/* Pixels containing contour */
+static gal_data_t *
+contour_pixels(gal_data_t *input, double level, size_t minmapsize,
+               int quietmmap)
+{
+  size_t one=1;
+  uint8_t *b, *a, *af;
+  int flags=GAL_ARITHMETIC_NUMOK;
+  gal_data_t *number, *thresh, *eroded;
+
+  /* Allocate the single-element dataset to use in arithmetic.*/
+  number=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &one, NULL, 1,
+                        -1, 1, NULL, NULL, NULL);
+  *(double *)(number->array)=level;
+
+  /* Only keep the pixels above the requested level, we are using the
+     arithmetic library to not have to worry about the type of the
+     input. */
+  thresh=gal_arithmetic(GAL_ARITHMETIC_OP_GT, 1, flags, input, number);
+
+  /* Erode the thresholded image by one. */
+  eroded=gal_binary_erode(thresh, 1, 1, 0);
+
+  /* Only keep the outer pixels. */
+  b=eroded->array;
+  af=(a=thresh->array)+thresh->size;
+  do if(*b++ == 1) *a=0; while(++a<af);
+
+  /* Clean up and return. */
+  gal_data_free(number);
+  gal_data_free(eroded);
+  return thresh;
+}
+
+
+
+
+
+/* Given the indexs of the contours, write them in the proper format. */
+static void
+contour_pgfplots(gal_data_t *edgeindexs, gal_data_t *input, float level,
+                 FILE *fp)
+{
+  double *xa, *ya;
+  gal_data_t *x, *y, *tmp;
+  size_t *s, *sf, w=input->dsize[1];
+
+  /* Go through each connected edge and add the contour positions. */
+  for(tmp=edgeindexs; tmp->next!=NULL; tmp=tmp->next)
+    if(tmp->size>10)
+      {
+        if(input->wcs)
+          {
+            {
+              /* Allocate the coordinate arrays. */
+              x=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &tmp->size,
+                               NULL, 0, -1, 1, NULL, NULL, NULL);
+              y=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &tmp->size,
+                               NULL, 0, -1, 1, NULL, NULL, NULL);
+
+              /* Fill in the coordinates. */
+              xa=x->array;
+              ya=y->array;
+              sf=(s=tmp->array)+tmp->size;
+              do {*xa++=*s%w+1; *ya++=*s/w+1;} while(++s<sf);
+
+              /* Convert the pixel positions to WCS. */
+              x->next=y;
+              gal_wcs_img_to_world(x, input->wcs, 1);
+
+              /* Write them. */
+              xa=x->array;
+              ya=y->array;
+              sf=(s=tmp->array)+tmp->size;
+              do
+                {fprintf(fp, "%.10f  %.10f  %f\n", *xa++, *ya++, level);}
+              while(++s<sf);
+
+              /* Clean up. */
+              gal_data_free(x);
+              gal_data_free(y);
+            }
+          }
+        else
+          {
+            sf=(s=tmp->array)+tmp->size;
+            do
+              {fprintf(fp, "%zu  %zu  %f\n", *s%w+1, *s/w+1, level);}
+            while(++s<sf);
+          }
+
+        /* To separate the different connected regions, we'll need to put
+           an empty line between them. */
+        fprintf(fp, "\n");
+      }
+}
+
+
+
+
+
+/* Contour for each level. */
+static void
+contour_level(gal_data_t *input, double level, FILE *fp,
+              size_t minmapsize, int quietmmap)
+{
+  gal_data_t *edge, *edgeindexs;
+
+  /* Find the edge pixels given this threshold. */
+  edge=contour_pixels(input, level, minmapsize, quietmmap);
+
+  /* Indexs of the edges (separated by groups of connected edges). */
+  edgeindexs=gal_binary_connected_indexs(edge, 2);
+
+  /* Make the PGFPlots contours. */
+  contour_pgfplots(edgeindexs, input, level, fp);
+
+  /* Clean up and return. */
+  gal_list_data_free(edgeindexs);
+  gal_data_free(edge);
+}
+
+
+
+
+
+void
+contour(struct statisticsparams *p)
+{
+  FILE *fp;
+  char *outname;
+  double *d, *df;
+  uint8_t keepinputdir=p->cp.keepinputdir;
+
+  /* Make sure the dataset is 2D. */
+  if(p->input->ndim!=2)
+    error(EXIT_FAILURE, 0, "contours are currently only supported for "
+          "2D datasets (images). The input dataset has %zu dimensions",
+          p->input->ndim);
+
+  /* Make sure the output doesn't exist. */
+  p->cp.keepinputdir = p->cp.output ? 1 : keepinputdir;
+  outname=gal_checkset_automatic_output(&p->cp,
+                                        ( p->cp.output
+                                          ? p->cp.output
+                                          : p->inputname ), "_contour.txt");
+  p->cp.keepinputdir=keepinputdir;
+
+  /* Open the output file. */
+  fp=fopen(outname, "w+");
+
+  /* Print basic information about the columns. */
+  fprintf(fp, "# %s Contour positions\n", PACKAGE_STRING);
+  fprintf(fp, "# Column 1: Coord_1 [position,f64] Position in first axis.\n");
+  fprintf(fp, "# Column 2: Coord_2 [position,f64] Position in second axis.\n");
+  fprintf(fp, "# Column 3: Level   [value,   f32] Contour level.\n");
+  fprintf(fp, "#\n");
+  fprintf(fp, "# Each connected contour is separated by an empty line.\n");
+  fprintf(fp, "# This format is recognized in PGFPlots (package of LaTeX).\n");
+
+  /* Estimate the contours for each level. */
+  df=(d=p->contour->array)+p->contour->size;
+  do
+    contour_level(p->input, *d, fp, p->cp.minmapsize,
+                  p->cp.quietmmap);
+  while(++d<df);
+
+  /* Clean up and cose the file. */
+  free(outname);
+  fclose(fp);
+}
diff --git a/bin/statistics/contour.h b/bin/statistics/contour.h
new file mode 100644
index 0000000..4d59c86
--- /dev/null
+++ b/bin/statistics/contour.h
@@ -0,0 +1,31 @@
+/*********************************************************************
+Statistics - Statistical analysis on input dataset.
+Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2019, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#ifndef CONTOUR_H
+#define CONTOUR_H
+
+void
+contour(struct statisticsparams *p);
+
+
+
+#endif
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 7e835cc..7455dc8 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -74,6 +74,7 @@ struct statisticsparams
   double            mirror;  /* Mirror value for hist and CFP.           */
   uint8_t              sky;  /* Find the Sky value over the image.       */
   uint8_t        sigmaclip;  /* So sigma-clipping over all dataset.      */
+  gal_data_t      *contour;  /* Levels to show contours.                 */
 
   size_t           numbins;  /* Number of bins in histogram or CFP.      */
   size_t      numasciibins;  /* Number of bins in ASCII plots.           */
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index db87377..3a36134 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -46,6 +46,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include "ui.h"
 #include "sky.h"
+#include "contour.h"
 #include "statistics.h"
 
 
@@ -1024,6 +1025,12 @@ statistics(struct statisticsparams *p)
       print_basic_info=0;
     }
 
+  if(p->contour)
+    {
+      contour(p);
+      print_basic_info=0;
+    }
+
   /* Print the ASCII plots if requested. */
   if(p->asciihist || p->asciicfp)
     {
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 86597e1..e590b51 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -904,7 +904,7 @@ ui_preparations(struct statisticsparams *p)
   ui_out_of_range_to_blank(p);
 
   /* If we are not to work on tiles, then re-order and change the input. */
-  if(p->ontile==0 && p->sky==0)
+  if(p->ontile==0 && p->sky==0 && p->contour==NULL)
     {
       /* Only keep the elements we want. */
       gal_blank_remove(p->input);
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 8b2422b..d1aab3b 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -46,7 +46,7 @@ enum program_args_groups
 /* Available letters for short options:
 
    a b e f j p v w x z
-   B G J L R W X Y
+   B G J L W X Y
 */
 enum option_keys_enum
 {
@@ -70,6 +70,7 @@ enum option_keys_enum
   UI_KEY_INTERPOLATE  = 'i',
   UI_KEY_SKY          = 'y',
   UI_KEY_KERNEL       = 'k',
+  UI_KEY_CONTOUR      = 'R',
 
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 35677c6..233ba84 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12639,6 +12639,19 @@ If the @option{--oneelempertile} option is called, 
then one element/pixel will b
 Otherwise, the output will have the same size as the input, but each element 
will have the value corresponding to that tile's value.
 If multiple single valued operations are called, then for each operation there 
will be one extension in the output FITS file.
 
+@item -R FLT[,FLT[,FLT...]]
+@itemx --contour=FLT[,FLT[,FLT...]]
+@cindex Contour
+@cindex Plot: contour
+Write the contours for the requested levels in a file ending with 
@file{_contour.txt}.
+It will have three columns: the first two are the coordinates of each point 
and the third is the level it belongs to (one of the input values).
+Each disconnected contour region will be separated by a blank line.
+This is the requested format for adding contours with PGFPlots in @LaTeX{}.
+If any other format can be useful for your work please let us know so we can 
add it.
+If the image has World Coordinate System information, the written coordinates 
will be in RA and Dec, otherwise, they will be in pixel coordinates.
+
+Note that currently, this is a very crude/simple implementation, please let us 
know if you find problematic situations so we can fix it.
+
 @item -y
 @itemx --sky
 Estimate the Sky value on each tile as fully described in @ref{Quantifying 
signal in a tile}.
@@ -23518,6 +23531,12 @@ non-zero pixels in @code{binary} will be considered as 
foreground (and will
 be labeled). Blank pixels in the input will also be blank in the output.
 @end deftypefun
 
+@deftypefun {gal_data_t *} gal_binary_connected_indexs(gal_data_t 
@code{*binary}, int @code{connectivity})
+Build a @code{gal_data_t} linked list, where each node of the list contains an 
array with indexs of the connected regions.
+Therefore the arrays of each node can have a different size.
+Note that the indexs will only be calculated on the pixels with a value of 1 
and internally, it will temporarily change the values to 2 (and return them 
back to 1 in the end).
+@end deftypefun
+
 @deftypefun {gal_data_t *} gal_binary_connected_adjacency_matrix (gal_data_t 
@code{*adjacency}, size_t @code{*numconnected})
 @cindex Adjacency matrix
 @cindex Matrix, adjacency
diff --git a/lib/binary.c b/lib/binary.c
index a10e2e4..bfec6c7 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -448,6 +448,97 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
 
 
 
+/* Put the indexs of connected labels in a list of `gal_data_t's, each with
+   a one-dimensional array that has the indexs of that connected
+   component.*/
+#define BINARY_CONINDEX_VAL 2
+gal_data_t *
+gal_binary_connected_indexs(gal_data_t *binary, int connectivity)
+{
+  uint8_t *b, *bf;
+  gal_data_t *lines=NULL;
+  size_t p, i, onelabnum, *onelabarr;
+  gal_list_sizet_t *Q=NULL, *onelab=NULL;
+  size_t *dinc=gal_dimension_increment(binary->ndim, binary->dsize);
+
+  /* Small sanity checks. */
+  if(binary->type!=GAL_TYPE_UINT8)
+    error(EXIT_FAILURE, 0, "%s: the input data set type must be `uint8'",
+          __func__);
+  if(binary->block)
+    error(EXIT_FAILURE, 0, "%s: currently, the input data structure to "
+          "must not be a tile", __func__);
+
+  /* Go over all the pixels and do a breadth-first search. */
+  b=binary->array;
+  for(i=0;i<binary->size;++i)
+    /* A pixel that has already been recorded is given a value of
+       `BINARY_CONINDEX_VAL'. */
+    if( b[i]==1 )
+      {
+        /* Add this pixel to the queue of pixels to work with. */
+       b[i]=BINARY_CONINDEX_VAL;
+        gal_list_sizet_add(&Q, i);
+        gal_list_sizet_add(&onelab, i);
+
+        /* While a pixel remains in the queue, continue labelling and
+           searching for neighbors. */
+        while(Q!=NULL)
+          {
+            /* Pop an element from the queue. */
+            p=gal_list_sizet_pop(&Q);
+
+            /* Go over all its neighbors and add them to the list if they
+               haven't already been labeled. */
+            GAL_DIMENSION_NEIGHBOR_OP(p, binary->ndim, binary->dsize,
+                                      connectivity, dinc,
+              {
+                if( b[nind]==1 )
+                  {
+                   b[nind]=BINARY_CONINDEX_VAL;
+                    gal_list_sizet_add(&Q, nind);
+                   gal_list_sizet_add(&onelab, nind);
+                  }
+              } );
+          }
+
+       /* Parsing has finished, put all the indexs into an array. */
+       onelabarr=gal_list_sizet_to_array(onelab, 1, &onelabnum);
+       gal_list_data_add_alloc(&lines, onelabarr, GAL_TYPE_SIZE_T, 1,
+                               &onelabnum, NULL, 0, -1, 1, NULL, NULL, NULL);
+
+       /* Clean up. */
+       gal_list_sizet_free(onelab);
+       onelab=NULL;
+      }
+
+  /* Reverse the order. */
+  gal_list_data_reverse(&lines);
+
+  /* For a check
+  {
+    gal_data_t *test=lines->next;
+    size_t *b, *bf;
+    bf=(b=test->array)+test->size;
+    do printf("%zu\n", *b++);
+    while(b<bf);
+    exit(0);
+  }
+  */
+
+  /* Set all the `2' values back to `1'. */
+  bf=(b=binary->array)+binary->size;
+  do if(*b==BINARY_CONINDEX_VAL) *b=1; while(++b<bf);
+
+  /* Clean up and return the total number. */
+  free(dinc);
+  return lines;
+}
+
+
+
+
+
 /* Given an adjacency matrix (which should be binary), find the number of
    connected objects and return an array of new labels for each old
    label.  */
diff --git a/lib/gnuastro/binary.h b/lib/gnuastro/binary.h
index 5a4b391..786cfd9 100644
--- a/lib/gnuastro/binary.h
+++ b/lib/gnuastro/binary.h
@@ -90,6 +90,9 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
                                 int connectivity);
 
 gal_data_t *
+gal_binary_connected_indexs(gal_data_t *binary, int connectivity);
+
+gal_data_t *
 gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
                                       size_t *numnewlabs);