[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 9c628b1: Statistics: 2D histograms can now be
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 9c628b1: Statistics: 2D histograms can now be made as FITS images |
|
Date: |
Tue, 29 Dec 2020 23:01:02 -0500 (EST) |
branch: master
commit 9c628b1085419bf2ffb1f4494ada03bc11c1bb7a
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Statistics: 2D histograms can now be made as FITS images
Until now 2D histograms were built as tables, that relied on visualization
tools to plot. But the visualization step was slow and too
versatile. However by definition, a 2D histogram is actually a 2D image!
With this commit, the '--histogram2d' option of Statistics now takes an
argument: either '--histogram2d=table' (same old table format) or
'--histogram2d=image' (which is new). When an image 2D histogram is
requested, a WCS is also added to the FITS image, enabling easy viewing of
the plot. See the description in the newly added "2D histogram as an image"
section of the book for more.
---
NEWS | 15 +++
bin/statistics/args.h | 10 +-
bin/statistics/main.h | 2 +-
bin/statistics/statistics.c | 120 ++++++++++++++++-----
bin/statistics/ui.c | 16 ++-
doc/gnuastro.texi | 249 +++++++++++++++++++++++++++++++++++---------
lib/gnuastro/wcs.h | 4 +
lib/statistics.c | 5 +-
lib/wcs.c | 48 +++++++++
9 files changed, 379 insertions(+), 90 deletions(-)
diff --git a/NEWS b/NEWS
index 9940cda..dcaffa2 100644
--- a/NEWS
+++ b/NEWS
@@ -101,6 +101,15 @@ See the end of the file for license conditions.
value to use for that radial interval. See the description of
'--customtable' for more.
+ Statistics:
+ - 2D histograms can now be built as a FITS image with a linear WCS that
+ contains axis values and box size. This allows using the power of FITS
+ viewers for plotting/inspecting distributions of two columns in a
+ table relative to each other (for example color-magnitude plots). You
+ can also convert these 2D histogram images to PDF or JPEG with
+ Gnuastro's ConvertType to directly use in your papers/reports (see the
+ newly added "2D histogram as an image" section of the book for more).
+
Table:
- New '--noblank' option will remove all rows in output table that have
at least one blank value in the specified columns. For example if
@@ -122,6 +131,7 @@ See the end of the file for license conditions.
- gal_pointer_allocate_ram_or_mmap: allocate space either in RAM or as a
memory-mapped file.
- gal_pointer_mmap_free: "free" (actually delete) the memory-mapped file.
+ - gal_wcs_create: create WCSLIB-compatible WCS from raw values.
- gal_wcs_coverage: Return the sky coverage of given image HDU.
- gal_wcs_dimension_name: return the name of the requested WCS dim.
@@ -170,6 +180,11 @@ See the end of the file for license conditions.
robustly, the default value of 'outliersigma' has been decreased to
'5' (previously it was 10).
+ Statistics:
+ - The '--histogram2d' now takes a string argument: either 'image' or
+ 'table'. For the old behavior please use '--histogram2d=table'. See
+ the new features above for the 'image' mode.
+
Table:
- Column arithmetic operators 'degree-to-ra' and 'degree-to-dec' will
return the sexagesimal format of '_h_m_s' and '_d_m_s'
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 710a753..a95dcc7 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -100,6 +100,8 @@ struct argp_option program_options[] =
+
+
/* Tessellation */
{
"interpolate",
@@ -426,13 +428,13 @@ struct argp_option program_options[] =
{
"histogram2d",
UI_KEY_HISTOGRAM2D,
+ "STR",
0,
- 0,
- "Save the 2D histogram in output.",
+ "2D histogram (as 'table' or 'image').",
UI_GROUP_PARTICULAR_STAT,
&p->histogram2d,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 488ae36..df9a0e2 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -72,7 +72,7 @@ struct statisticsparams
uint8_t asciihist; /* Print an ASCII histogram. */
uint8_t asciicfp; /* Print an ASCII cumulative frequency plot.*/
uint8_t histogram; /* Save histogram in output. */
- uint8_t histogram2d; /* Save 2D-histogram in output. */
+ char *histogram2d; /* Save 2D-histogram as image or table. */
uint8_t cumulative; /* Save cumulative distibution in output. */
double mirror; /* Mirror value for hist and CFP. */
uint8_t sky; /* Find the Sky value over the image. */
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index d28dfa3..06a7473 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <string.h>
#include <stdlib.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
#include <gnuastro/tile.h>
#include <gnuastro/blank.h>
@@ -586,37 +587,55 @@ ascii_plots(struct statisticsparams *p)
/*******************************************************************/
/******* Histogram and cumulative frequency tables ***********/
/*******************************************************************/
-void
-write_output_table(struct statisticsparams *p, gal_data_t *table,
- char *suf, char *contents)
+static char *
+statistics_output_name(struct statisticsparams *p, char *suf, int *isfits)
{
- char *output;
int use_auto_output=0;
- char *fix, *suffix=NULL, *tmp;
- gal_list_str_t *comments=NULL;
-
+ char *out, *fix, *suffix=NULL;
/* Automatic output should be used when no output name was specified or
we have more than one output file. */
use_auto_output = p->cp.output ? (p->numoutfiles>1 ? 1 : 0) : 1;
-
/* Set the 'fix' and 'suffix' strings. Note that 'fix' is necessary in
every case, even when no automatic output is to be used. Since it is
used to determine the format of the output. */
- fix = ( p->cp.output
- ? gal_fits_name_is_fits(p->cp.output) ? "fits" : "txt"
- : "txt" );
+ fix = ( *isfits
+ ? "fits"
+ : ( p->cp.output
+ ? gal_fits_name_is_fits(p->cp.output) ? "fits" : "txt"
+ : "txt" ) );
if(use_auto_output)
if( asprintf(&suffix, "%s.%s", suf, fix)<0 )
error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
-
/* Make the output name. */
- output = ( use_auto_output
- ? gal_checkset_automatic_output(&p->cp, p->inputname, suffix)
- : p->cp.output );
+ out = ( use_auto_output
+ ? gal_checkset_automatic_output(&p->cp, p->inputname, suffix)
+ : p->cp.output );
+
+ /* See if it is a FITS file. */
+ *isfits=strcmp(fix, "fits")==0;
+
+ /* Clean up and return. */
+ if(suffix) free(suffix);
+ return out;
+}
+
+
+
+
+void
+write_output_table(struct statisticsparams *p, gal_data_t *table,
+ char *suf, char *contents)
+{
+ int isfits=0;
+ char *tmp, *output;
+ gal_list_str_t *comments=NULL;
+
+ /* Set the output. */
+ output=statistics_output_name(p, suf, &isfits);
/* Write the comments, NOTE: we are writing the first two in reverse of
the order we want them. They will later be freed as part of the list's
@@ -628,7 +647,7 @@ write_output_table(struct statisticsparams *p, gal_data_t
*table,
error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
gal_list_str_add(&comments, tmp, 0);
- if(strcmp(fix, "fits")) /* The intro info will be in FITS files anyway.*/
+ if(!isfits) /* The intro info will be in FITS files anyway.*/
gal_table_comments_add_intro(&comments, PROGRAM_STRING, &p->rawtime);
@@ -639,7 +658,7 @@ write_output_table(struct statisticsparams *p, gal_data_t
*table,
/* Write the configuration information if we have a FITS output. */
- if(!strcmp(fix, "fits"))
+ if(isfits)
{
gal_fits_key_write_filename("input", p->inputname, &p->cp.okeys, 1);
gal_fits_key_write_config(&p->cp.okeys, "Statistics configuration",
@@ -653,7 +672,6 @@ write_output_table(struct statisticsparams *p, gal_data_t
*table,
/* Clean up. */
- if(suffix) free(suffix);
gal_list_str_free(comments, 1);
if(output!=p->cp.output) free(output);
}
@@ -728,26 +746,74 @@ save_hist_and_or_cfp(struct statisticsparams *p)
static void
histogram_2d(struct statisticsparams *p)
{
- char *suf, *contents;
- gal_data_t *hist2d, *bins;
+ int isfits=1;
+ int32_t *imgarr;
+ double *d1, *d2;
+ uint32_t *histarr;
gal_data_t *range1, *range2;
+ gal_data_t *img, *hist2d, *bins;
+ size_t i, j, dsize[2], nb1=p->numbins, nb2=p->numbins2;
+ char *output, suf[]="-hist2d", contents[]="2D Histogram";
+
+ /* WCS-related arrays. */
+ char *cunit[2], *ctype[2];
+ double crpix[2], crval[2], cdelt[2], pc[4]={1,0,0,1};
/* Set the bins for each dimension */
range1=set_bin_range_params(p, 1);
range2=set_bin_range_params(p, 2);
- bins=gal_statistics_regular_bins(p->input, range1, p->numbins,
+ bins=gal_statistics_regular_bins(p->input, range1, nb1,
p->onebinstart);
bins->next=gal_statistics_regular_bins(p->input->next, range2,
- p->numbins2,
- p->onebinstart2);
+ nb2, p->onebinstart2);
/* Build the 2D histogram. */
hist2d=gal_statistics_histogram2d(p->input, bins);
- /* Write the output. */
- suf="-hist2d";
- contents="2D Histogram";
- write_output_table(p, hist2d, suf, contents);
+ /* Write the histogram into a 2D FITS image. Note that in the FITS image
+ standard, the first axis is the fastest array (unlike the default
+ format we have adopted for the tables). */
+ if(!strcmp(p->histogram2d,"image"))
+ {
+ /* Allocate the 2D image array. */
+ dsize[0]=nb1;
+ dsize[1]=nb2;
+ histarr=hist2d->next->next->array;
+ img=gal_data_alloc(NULL, GAL_TYPE_INT32, 2, dsize, NULL, 0,
+ p->cp.minmapsize, p->cp.quietmmap,
+ NULL, NULL, NULL);
+
+ /* Fill the array values. */
+ imgarr=img->array;
+ for(i=0;i<nb1;++i)
+ for(j=0;j<nb2;++j)
+ imgarr[i*nb2+j]=histarr[j*nb1+i];
+
+ /* Set the WCS. */
+ d1=bins->array;
+ d2=bins->next->array;
+ crpix[0] = 1; crpix[1] = 1;
+ crval[0] = d1[0]; crval[1] = d2[0];
+ cdelt[0] = d1[1]-d1[0]; cdelt[1] = d2[1]-d2[0];
+ cunit[0] = p->input->unit; cunit[1] = p->input->next->unit;
+ ctype[0] = p->input->name ? p->input->name : "X";
+ ctype[1] = p->input->next->name ? p->input->next->name : "Y";
+ img->wcs=gal_wcs_create(crpix, crval, cdelt, pc, cunit, ctype, 2);
+
+ /* Write the output. */
+ output=statistics_output_name(p, suf, &isfits);
+ gal_fits_img_write(img, output, NULL, PROGRAM_STRING);
+ gal_fits_key_write_filename("input", p->inputname, &p->cp.okeys, 1);
+ gal_fits_key_write_config(&p->cp.okeys, "Statistics configuration",
+ "STATISTICS-CONFIG", output, "0");
+
+ /* Clean up and let the user know that the histogram is built. */
+ gal_data_free(img);
+ if(!p->cp.quiet) printf("%s created.\n", output);
+ }
+ /* Write 2D histogram as a table. */
+ else
+ write_output_table(p, hist2d, suf, contents);
/* Clean up. */
gal_list_data_free(hist2d);
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 25b4cc4..b576ef7 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -26,6 +26,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <errno.h>
#include <error.h>
#include <stdio.h>
+#include <string.h>
#include <gnuastro/txt.h>
#include <gnuastro/wcs.h>
@@ -503,11 +504,16 @@ ui_read_check_only_options(struct statisticsparams *p)
error(EXIT_FAILURE, 0, "'--numbins' isn't set. When the histogram or "
"cumulative frequency plots are requested, the number of bins "
"('--numbins') is necessary");
- if( p->histogram2d && p->numbins2==0 )
- error(EXIT_FAILURE, 0, "'--numbins2' isn't set. When a 2D histogram "
- "is requested, the number of bins in the second dimension "
- "('--numbins2') is also necessary");
-
+ if( p->histogram2d )
+ {
+ if( p->numbins2==0 )
+ error(EXIT_FAILURE, 0, "'--numbins2' isn't set. When a 2D histogram "
+ "is requested, the number of bins in the second dimension "
+ "('--numbins2') is also necessary");
+ if( strcmp(p->histogram2d,"table") && strcmp(p->histogram2d,"image") )
+ error(EXIT_FAILURE, 0, "the value to '--histogram2d' can either be "
+ "'table' or 'image'");
+ }
/* If an ascii plot is requested, check if the ascii number of bins and
height are given. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 93abd9d..1e2717b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -491,6 +491,11 @@ Statistics
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
+2D Histograms
+
+* 2D histogram as a table:: Format and usage in table format.
+* 2D histogram as an image:: Format and usage in image format
+
Sky value
* Sky value definition:: Definition of the Sky/reference value.
@@ -13696,19 +13701,44 @@ But when the range is determined from the actual
dataset (none of these options
@node 2D Histograms, Sigma clipping, Histogram and Cumulative Frequency Plot,
Statistics
@subsection 2D Histograms
-
+@cindex 2D histogram
+@cindex Histogram, 2D
In @ref{Histogram and Cumulative Frequency Plot} the concept of histograms
were introduced on a single dataset.
-However, especially when doing high-level science on tables, the distribution
in a 2D space may be of interest (for example a color-magnitude diagram).
-But the number of points may be too large for a simple scatter plot to show
the concentration of the points: they will all fall over each other and just
make a large connected region that will hide potentially interesting behaviors.
+But they are only useful for viewing the distribution of a single variable
(column in a table).
+In many contexts, the distribution of two variables in relation to each other
may be of interest.
+For example the color-magnitude diagrams in astronomy, where the horizontal
axis is the luminosity or magnitude of an object, and the vertical axis is the
color.
+Scatter plots are useful to see these relations between the objects of
interest when the number of the objects is small.
+
+As the density of points in the scatter plot increases, the points will fall
over each other and just make a large connected region hide potentially
interesting behaviors/correlations in the densest regions.
This is where 2D histograms can become very useful.
-The desired 2D region is broken up into 2D bins (boxes) and the number of
points falling within each box is returned.
-Added with a color-bar, you can now clearly see the distribution.
+A 2D histogram is composed of 2D bins (boxes or pixels), just as a 1D
histogram consists of 1D bins (lines).
+The number of points falling within each box/pixel will then be the value of
that box.
+Added with a color-bar, you can now clearly see the distribution independent
of the density of points (for example you can even visualize it in log-scale if
you want).
Gnuastro's Statistics program has the @option{--histogram2d} option for this
task.
-Its output will be three columns that have the centers of every box in both
dimensions.
-The first column is the central box coordinates in the first dimension, the
second has values along the second dimension and the third has the number of
input points that fall within each box.
-You can specify the number of bins along each dimension through the
@option{--numbins} (for first input column) an @option{--numbins2} (for second
input column).
-The output file from this command can then be given to any plotting tool to
visualize the distribution.
+It takes a single argument (either @code{table} or @code{image}) that
specifies the format of the output 2D histogram.
+The two formats will be reviewed separately in the sub-sections below.
+But let's start with the generalities that are common to both (related to the
input, not the output).
+
+You can specify the two columns to be shown using the @option{--column} (or
@option{-c}) option.
+So if you want to plot the color-magnitude diagram from a table with the
@code{MAG-R} column on the horizontal and @code{COLOR-G-R} on the vertical
column, you can use @option{--column=MAG-r,COLOR-G-r}.
+The number of bins along each dimension can be set with @option{--numbins}
(for first input column) and @option{--numbins2} (for second input column).
+
+Without specifying any range, the full range of values will be used in each
dimension.
+If you only want to focus on a certain interval of the values in the columns
in any dimension you can use the @option{--greaterequal} and
@option{--lessthan} options to limit the values along the first/horizontal
dimension and @option{--greaterequal2} and @option{--lessthan2} options for the
second/vertical dimension.
+
+@menu
+* 2D histogram as a table:: Format and usage in table format.
+* 2D histogram as an image:: Format and usage in image format
+@end menu
+
+@node 2D histogram as a table, 2D histogram as an image, 2D Histograms, 2D
Histograms
+@subsubsection 2D histogram as a table
+
+When called with the @option{--histogram=table} option, Statistics will output
a table file with three columns that have the information of every box as a
column.
+If you asked for @option{--numbins=N} and @option{--numbins2=M}, all three
columns will have @mymath{M\times N} rows (one row for every box/pixel of the
2D histogram).
+The first and second columns are the position of the box along the first and
second dimensions.
+The third column has the number of input points that fall within that
box/pixel.
For example, you can make high-quality plots within your paper (using the same
@LaTeX{} engine, thus blending very nicely with your text) using
@url{https://ctan.org/pkg/pgfplots, PGFPlots}.
Below you can see one such minimal example, using your favorite text editor,
save it into a file, make the two small corrections in it, then run the
commands shown at the top.
@@ -13779,6 +13809,121 @@ The second one (@code{FILE.txt}) should be replaced
with the name of the file ge
\end@{document@}
@end example
+@node 2D histogram as an image, , 2D histogram as a table, 2D Histograms
+@subsubsection 2D histogram as an image
+
+When called with the @option{--histogram=image} option, Statistics will output
a FITS file with an image/array extension.
+If you asked for @option{--numbins=N} and @option{--numbins2=M} the image will
have a size of @mymath{N\times M} pixels (one pixel per 2D bin).
+Also, the FITS image will have a linear WCS that is scaled to the 2D bin size
along each dimension.
+So when you hover your mouse over any part of the image with a FITS viewer
(for example SAO DS9), besides the number of points in each pixel, you can
directly also see ``coordinates'' of the pixels along the two axes.
+You can also use the optimized and fast FITS viewer features for many aspects
of inspecting the distributions (which we won't go into futher here).
+
+@cindex Color-magnitude diagram
+@cindex Diagaram, Color-magnitude
+For example let's assume you want to derive the color-magnitude diagram (CMD)
of the @url{http://uvudf.ipac.caltech.edu, UVUDF survey}.
+You can run the first command below to download the table with magnitudes of
objects in many filters and run the second command to see general column
metadata after it is downloaded.
+
+@example
+$ wget http://asd.gsfc.nasa.gov/UVUDF/uvudf_rafelski_2015.fits.gz
+$ asttable uvudf_rafelski_2015.fits.gz -i
+@end example
+
+Let's assume you want to find the color to be between the @code{F775W} and
@code{F105W} filters (roughly corresponding to the g and r filters in
ground-based imaging).
+Because the original table doesn't have a color columns, you can use the
@ref{Column arithmetic} feature of the Table program for deriving a new table
with the @code{F775W} magnitude in one column and the difference between the
@code{F606W} and @code{F775W} on the other column.
+With the second command, you can see the actual values if you like.
+
+@example
+$ asttable uvudf_rafelski_2015.fits.gz -cMAG_F775W \
+ -c'arith MAG_F606W MAG_F775W -' \
+ --colmetadata=ARITH_1,F606W-F775W,"AB mag" -ocmd.fits
+$ asttable cmd.fits
+@end example
+
+@noindent
+You can construct your 2D histogram image as a FITS file with this command
(assuming you want @code{F775W} magnitudes between 20 and 30, colors between -1
and 3 and 100 bins in each dimension).
+
+@example
+aststatistics cmd.fits -cF775W,F606W-F775W --histogram2d=image \
+ --numbins=100 --greaterequal=22 --lessthan=30 \
+ --numbins2=100 --greaterequal2=-0.5 --lessthan2=3 \
+ --manualbinrange --output=cmd-2d-hist.fits
+@end example
+
+@noindent
+If you have SAO DS9, you can now open this FITS file as a normal FITS image,
for example with the command below.
+Try hovering/zooming over the pixels: not only will you see the number of
objects in the UVUDF catalog that fall in each bin, but you also see the
@code{F775W} magnitude and color of that pixel also.
+
+@example
+$ ds9 cmd-2d-hist.fits -cmap sls -zoom to fit
+@end example
+
+@noindent
+With the first command below, you can activate the grid feature of DS9 to
actually see the coordinate grid, as well as values on each line.
+With the second command, DS9 will even read the labels of the axises and use
them to generate an almost publication-ready plot.
+
+@example
+$ ds9 cmd-2d-hist.fits -cmap sls -zoom to fit -grid yes
+$ ds9 cmd-2d-hist.fits -cmap sls -zoom to fit -grid yes \
+ -grid type publication
+@end example
+
+If you are happy with the grid and coloring and etc, you can also save this as
a JPEG image to use in your documents with these extra DS9 options (DS9 will
quit as soon as it builds the JPEG image):
+
+@example
+$ ds9 cmd-2d-hist.fits -cmap sls -zoom 4 -grid yes \
+ -grid type publication -saveimage cmd-2d.jpeg -quit
+@end example
+
+@cindex PGFPlots (@LaTeX{} package)
+While this is good for a fast report in the weekly meeting, for your paper,
you want to show something with higher quality.
+For that, you can use the @LaTeX{} PGFPlots to add axises in the same font as
your text, sharp grids and many other elegant/powerful features.
+But to load the 2D histogram into PGFPlots first you need to convert the FITS
image into a more standard format, for example PDF.
+We'll use Gnuastro's @ref{ConvertType} for this, and use the
@code{sls-inverse} color map (which will set the color white for pixles without
a value):
+
+@example
+$ astconvertt cmd-2d-hist.fits --colormap=sls-inverse \
+ --borderwidth=0 -ocmd-2d-hist.pdf
+@end example
+
+@noindent
+Below you can see a minimally working example of how to add axis numbers,
labels and a grid to the PDF generated above.
+Notice the @code{xmin}, @code{xmax}, @code{ymin}, @code{ymax} values and how
they are the same as the range specified above.
+
+@example
+\documentclass@{article@}
+\usepackage@{pgfplots@}
+\dimendef\prevdepth=0
+\begin@{document@}
+
+You can write all you want here...
+
+\begin@{tikzpicture@}
+ \begin@{axis@}[
+ enlargelimits=false,
+ grid,
+ axis on top,
+ width=\linewidth,
+ height=\linewidth,
+ xlabel=@{Magnitude (F775W)@},
+ ylabel=@{Color (F606W-F775W)@}]
+
+ \addplot graphics[xmin=22, xmax=30, ymin=-0.5, ymax=3]
+ @{cmd-2d-hist.pdf@};
+ \end@{axis@}
+\end@{tikzpicture@}
+\end@{document@}
+@end example
+
+@noindent
+Copy and paste the code above into a plain-text file called
@file{cmd-report.tex} and run this command to build your PDF (assuming you have
LaTeX and PGFPlots).
+
+@example
+$ pdflatex cmd-report.tex
+@end example
+
+The improved quality, blending in with the text, vector-graphics resolution
and other features make this plot pleasing to the eye, and let your readers
focus on the main point of your scientific argument.
+
+
@node Sigma clipping, Sky value, 2D Histograms, Statistics
@subsection Sigma clipping
@@ -24024,31 +24169,29 @@ the returned structure with WCSLIB's @code{wcsvfree}
keyword:
status = wcsvfree(&nwcs,&wcs);
@end example
-If you don't want to search the full FITS header for WCS-related FITS
-keywords (for example due to conflicting keywords), but only a specific
-range of the header keywords you can use the @code{hstartwcs} and
-@code{hendwcs} arguments to specify the keyword number range (counting from
-zero). If @code{hendwcs} is larger than @code{hstartwcs}, then only
-keywords in the given range will be checked. Hence, to ignore this feature
-(and search the full FITS header), give both these arguments the same
-value.
+If you don't want to search the full FITS header for WCS-related FITS keywords
(for example due to conflicting keywords), but only a specific range of the
header keywords you can use the @code{hstartwcs} and @code{hendwcs} arguments
to specify the keyword number range (counting from zero).
+If @code{hendwcs} is larger than @code{hstartwcs}, then only keywords in the
given range will be checked.
+Hence, to ignore this feature (and search the full FITS header), give both
these arguments the same value.
-If the WCS information couldn't be read from the FITS file, this function
-will return a @code{NULL} pointer and put a zero in @code{nwcs}. A WCSLIB
-error message will also be printed in @code{stderr} if there was an error.
+If the WCS information couldn't be read from the FITS file, this function will
return a @code{NULL} pointer and put a zero in @code{nwcs}.
+A WCSLIB error message will also be printed in @code{stderr} if there was an
error.
-This function is just a wrapper over WCSLIB's @code{wcspih} function which
-is not thread-safe. Therefore, be sure to not call this function
-simultaneously (over multiple threads).
+This function is just a wrapper over WCSLIB's @code{wcspih} function which is
not thread-safe.
+Therefore, be sure to not call this function simultaneously (over multiple
threads).
@end deftypefun
@deftypefun {struct wcsprm *} gal_wcs_read (char @code{*filename}, char
@code{*hdu}, size_t @code{hstartwcs}, size_t @code{hendwcs}, int @code{*nwcs})
-[@strong{Not thread-safe}] Return the WCSLIB structure that is read from
-the HDU/extension @code{hdu} of the file @code{filename}. Also put the
-number of coordinate representations found into the space that @code{nwcs}
-points to. Please see @code{gal_wcs_read_fitsptr} for more.
+[@strong{Not thread-safe}] Return the WCSLIB structure that is read from the
HDU/extension @code{hdu} of the file @code{filename}.
+Also put the number of coordinate representations found into the space that
@code{nwcs} points to.
+Please see @code{gal_wcs_read_fitsptr} for more.
@end deftypefun
+@deftypefun {struct wcsprm *} gal_wcs_create (double @code{*crpix}, double
@code{*crval}, double @code{*cdelt}, double @code{*pc}, char @code{**cunit},
char @code{**ctype}, size_t @code{ndim})
+Given all the most common standard components of the WCS standard, construct a
@code{struct wcsprm}, initialize and set it for future processing.
+See the FITS WCS standard for more on these keywords.
+All the arrays must have @code{ndim} elements with them except for @code{pc}
which should have @code{ndim*ndim} elements (a square matrix).
+Also, @code{cunit} and @code{ctype} are arrays of strings.
+@end deftypefun
@deftypefun {char *} gal_wcs_dimension_name (struct wcsprm @code{*wcs}, size_t
@code{dimension})
Return an allocated string array (that should be freed later) containing the
first part of the @code{CTYPEi} FITS keyword (which contains the dimension name
in the FITS standard).
@@ -26154,39 +26297,43 @@ shifted to accommodate this request.
@deftypefun {gal_data_t *} gal_statistics_histogram (gal_data_t @code{*input},
gal_data_t @code{*bins}, int @code{normalize}, int @code{maxone})
-Make a histogram of all the elements in the given dataset with bin values
-that are defined in the @code{inbins} structure (see
-@code{gal_statistics_regular_bins}, they currently have to be equally
-spaced). @code{inbins} is not mandatory, if you pass a @code{NULL} pointer,
-the bins structure will be built within this function based on the
-@code{numbins} input. As a result, when you have already defined the bins,
-@code{numbins} is not used.
+@cindex Histogram
+Make a histogram of all the elements in the given dataset with bin values that
are defined in the @code{inbins} structure (see
@code{gal_statistics_regular_bins}, they currently have to be equally spaced).
+@code{inbins} is not mandatory, if you pass a @code{NULL} pointer, the bins
structure will be built within this function based on the @code{numbins} input.
+As a result, when you have already defined the bins, @code{numbins} is not
used.
-Let's write the center of the @mymath{i}th element of the bin array as
-@mymath{b_i}, and the fixed half-bin width as @mymath{h}. Then element
-@mymath{j} of the input array (@mymath{in_j}) will be counted in
-@mymath{b_i} if @mymath{(b_i-h) \le in_j < (b_i+h)}. However, if
-@mymath{in_j} is somewhere in the last bin, the condition changes to
-@mymath{(b_i-h) \le in_j \le (b_i+h)}.
+Let's write the center of the @mymath{i}th element of the bin array as
@mymath{b_i}, and the fixed half-bin width as @mymath{h}.
+Then element @mymath{j} of the input array (@mymath{in_j}) will be counted in
@mymath{b_i} if @mymath{(b_i-h) \le in_j < (b_i+h)}.
+However, if @mymath{in_j} is somewhere in the last bin, the condition changes
to @mymath{(b_i-h) \le in_j \le (b_i+h)}.
+
+If @code{normalize!=0}, the histogram will be ``normalized'' such that the sum
of the counts column will be one. In other words, all the counts in every bin
will be divided by the total number of counts. If @code{maxone!=0}, the
histogram's maximum count will be 1. In other words, the counts in every bin
will be divided by the value of the maximum.
@end deftypefun
+@deftypefun {gal_data_t *} gal_statistics_histogram2d (gal_data_t
@code{*input}, gal_data_t @code{*bins})
+@cindex Histogram, 2D
+@cindex 2D histogram
+This function is very similar to @code{gal_statistics_histogram}, but will
build a 2D histogram (count how many of the elements of @code{input} have a
within a 2D box.
+The bins comprising the first dimension of the 2D box are defined by
@code{bins}.
+The bins of the second dimension are defined by @code{bins->next} (@code{bins}
is a @ref{List of gal_data_t}).
+Both the @code{bin} and @code{bin->next} can be created with
@code{gal_statistics_regular_bins}.
+
+This function returns a list of @code{gal_data_t} with three nodes/columns, so
you can directly write them into a table (see @ref{Table input output}).
+Assuming @code{bins} has @mymath{N1} bins and @code{bins->next} has
@mymath{N2} bins, each node/column of the returned output is a 1D array with
@mymath{N1\times N2} elements.
+The first and second columns are the center of the 2D bin along the first and
second dimensions and have a @code{double} data type.
+The third column is the 2D histogram (the number of input elements that have a
value within that 2D bin) and has a @code{uint32} data type (see @ref{Numeric
data types}).
+@end deftypefun
@deftypefun {gal_data_t *} gal_statistics_cfp (gal_data_t @code{*input},
gal_data_t @code{*bins}, int @code{normalize})
Make a cumulative frequency plot (CFP) of all the elements in @code{input}
with bin values that are defined in the @code{bins} structure (see
@code{gal_statistics_regular_bins}).
-The CFP is built from the histogram: in each bin, the value is the sum of
-all previous bins in the histogram. Thus, if you have already calculated
-the histogram before calling this function, you can pass it onto this
-function as the data structure in @code{bins->next} (see @code{List of
-gal_data_t}). If @code{bin->next!=NULL}, then it is assumed to be the
-histogram. If it is @code{NULL}, then the histogram will be calculated
-internally and freed after the job is finished.
-
-When a histogram is given and it is normalized, the CFP will also be
-normalized (even if the normalized flag is not set here): note that a
-normalized CFP's maximum value is 1.
+The CFP is built from the histogram: in each bin, the value is the sum of all
previous bins in the histogram.
+Thus, if you have already calculated the histogram before calling this
function, you can pass it onto this function as the data structure in
@code{bins->next} (see @code{List of gal_data_t}).
+If @code{bin->next!=NULL}, then it is assumed to be the histogram.
+If it is @code{NULL}, then the histogram will be calculated internally and
freed after the job is finished.
+
+When a histogram is given and it is normalized, the CFP will also be
normalized (even if the normalized flag is not set here): note that a
normalized CFP's maximum value is 1.
@end deftypefun
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 7457411..8fdb0ef 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -85,6 +85,10 @@ struct wcsprm *
gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
size_t hendwcs, int *nwcs);
+struct wcsprm *
+gal_wcs_create(double *crpix, double *crval, double *cdelt,
+ double *pc, char **cunit, char **ctype, size_t ndim);
+
char *
gal_wcs_dimension_name(struct wcsprm *wcs, size_t dimension);
diff --git a/lib/statistics.c b/lib/statistics.c
index 2f55605..ade8602 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1911,9 +1911,10 @@ gal_statistics_histogram(gal_data_t *input, gal_data_t
*bins, int normalize,
gal_data_t *
gal_statistics_histogram2d(gal_data_t *input, gal_data_t *bins)
{
+ uint32_t *h;
double *o1, *o2;
gal_data_t *tmp, *out;
- size_t i, j, *h, bsizea, bsizeb, outsize;
+ size_t i, j, bsizea, bsizeb, outsize;
double *da, *db, binwidtha, binwidthb, mina, minb, maxa, maxb;
/* Basic sanity checks */
@@ -1954,7 +1955,7 @@ gal_statistics_histogram2d(gal_data_t *input, gal_data_t
*bins)
"bin_dim2", input->next->unit,
"Bin centers along second axis.");
out->next=tmp;
- tmp=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1, &outsize,
+ tmp=gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &outsize,
NULL, 1, input->minmapsize, input->quietmmap,
"hist_number", "counts",
"Number of data points within each 2D-bin (box).");
diff --git a/lib/wcs.c b/lib/wcs.c
index 8a401eb..ea9a221 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -376,6 +376,54 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
+struct wcsprm *
+gal_wcs_create(double *crpix, double *crval, double *cdelt,
+ double *pc, char **cunit, char **ctype, size_t ndim)
+{
+ size_t i;
+ int status;
+ struct wcsprm *wcs;
+ double equinox=2000.0f;
+
+ /* Allocate the memory necessary for the wcsprm structure. */
+ errno=0;
+ wcs=malloc(sizeof *wcs);
+ if(wcs==NULL)
+ error(EXIT_FAILURE, errno, "%zu for wcs in preparewcs", sizeof *wcs);
+
+ /* Initialize the structure (allocate all its internal arrays). */
+ wcs->flag=-1;
+ if( (status=wcsini(1, ndim, wcs)) )
+ error(EXIT_FAILURE, 0, "wcsini error %d: %s",
+ status, wcs_errmsg[status]);
+
+ /* Fill in all the important WCS structure parameters. */
+ wcs->altlin = 0x1;
+ wcs->equinox = equinox;
+ for(i=0;i<ndim;++i)
+ {
+ wcs->crpix[i] = crpix[i];
+ wcs->crval[i] = crval[i];
+ wcs->cdelt[i] = cdelt[i];
+ strcpy(wcs->cunit[i], cunit[i]);
+ strcpy(wcs->ctype[i], ctype[i]);
+ }
+ for(i=0;i<ndim*ndim;++i) wcs->pc[i]=pc[i];
+
+ /* Set up the wcs structure with the constants defined above. */
+ status=wcsset(wcs);
+ if(status)
+ error(EXIT_FAILURE, 0, "wcsset error %d: %s", status,
+ wcs_errmsg[status]);
+
+ /* Return the output WCS. */
+ return wcs;
+}
+
+
+
+
+
/* Extract the dimension name from CTYPE. */
char *
gal_wcs_dimension_name(struct wcsprm *wcs, size_t dimension)
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 9c628b1: Statistics: 2D histograms can now be made as FITS images,
Mohammad Akhlaghi <=