[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 99f3c71: Mode functions moved to statistics.h
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 99f3c71: Mode functions moved to statistics.h |
|
Date: |
Tue, 20 Sep 2016 23:44:03 +0000 (UTC) |
branch: master
commit 99f3c71d4711384545c9132e33f13ed1215947f4
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Mode functions moved to statistics.h
The user-callable mode calculation functions are now defined in
`lib/statistics.c' instead of `lib/mode.c'. `lib/gnuastro/mode.h' is also
no longer installed, it is just distributed and internally included.
---
lib/Makefile.am | 17 ++--
lib/gnuastro/mode.h | 79 +++------------
lib/gnuastro/statistics.h | 42 ++++++++
lib/mesh.c | 1 -
lib/mode.c | 213 ++---------------------------------------
lib/statistics.c | 213 ++++++++++++++++++++++++++++++++++++++++-
src/imgstat/imgstat.c | 34 +++----
src/noisechisel/thresh.c | 6 +-
src/subtractsky/subtractsky.c | 6 +-
9 files changed, 312 insertions(+), 299 deletions(-)
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 965390c..8a8145d 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -39,12 +39,11 @@ libgnuastro_la_SOURCES = array.c box.c checkset.c
configfiles.c fits.c \
# Installed headers, note that we are not blindly including all `.h' files
# in the $(headersdir) directory. Some of the header files don't need to be
# installed.
-pkginclude_HEADERS = gnuastro/gnuastro.h $(headersdir)/array.h \
- $(headersdir)/box.h $(headersdir)/fits.h $(headersdir)/linkedlist.h \
- $(headersdir)/mesh.h $(headersdir)/mode.h $(headersdir)/polygon.h \
- $(headersdir)/qsort.h $(headersdir)/spatialconvolve.h \
- $(headersdir)/statistics.h $(headersdir)/threads.h $(headersdir)/wcs.h \
- $(headersdir)/txtarray.h
+pkginclude_HEADERS = gnuastro/gnuastro.h $(headersdir)/array.h \
+ $(headersdir)/box.h $(headersdir)/fits.h $(headersdir)/linkedlist.h \
+ $(headersdir)/mesh.h $(headersdir)/polygon.h $(headersdir)/qsort.h \
+ $(headersdir)/spatialconvolve.h $(headersdir)/statistics.h \
+ $(headersdir)/threads.h $(headersdir)/wcs.h $(headersdir)/txtarray.h
# Files to distribute in the tarball. For the headers: these are internal
@@ -53,9 +52,9 @@ pkginclude_HEADERS = gnuastro/gnuastro.h
$(headersdir)/array.h \
# the Makefile, they will not distributed, so we need to explicitly tell
# Automake to distribute them here.
headersdir=$(top_srcdir)/lib/gnuastro
-EXTRA_DIST = gnuastro.pc.in gnuastro.h.in $(headersdir)/commonargs.h \
- $(headersdir)/commonparams.h $(headersdir)/fixedstringmacros.h \
- $(headersdir)/neighbors.h $(headersdir)/checkset.h \
+EXTRA_DIST = gnuastro.pc.in gnuastro.h.in $(headersdir)/commonargs.h \
+ $(headersdir)/commonparams.h $(headersdir)/fixedstringmacros.h \
+ $(headersdir)/mode.h $(headersdir)/neighbors.h $(headersdir)/checkset.h \
$(headersdir)/configfiles.h $(headersdir)/timing.h $(headersdir)/README
diff --git a/lib/gnuastro/mode.h b/lib/gnuastro/mode.h
index 1c38afa..3838f4c 100644
--- a/lib/gnuastro/mode.h
+++ b/lib/gnuastro/mode.h
@@ -20,76 +20,27 @@ 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 __GAL_MODE_H__
-#define __GAL_MODE_H__
+#ifndef __MODE_H__
+#define __MODE_H__
-/* Include other headers if necessary here. Note that other header files
- must be included before the C++ preparations below */
+#include <gnuastro/statistics.h>
+#define MODE_GOLDEN_RATIO 1.618034f
+#define MODE_TWO_TAKE_GOLDEN_RATIO 0.38197f
+#define MODE_MIRROR_IS_ABOVE_RESULT (size_t)(-1)
-
-/* C++ Preparations */
-#undef __BEGIN_C_DECLS
-#undef __END_C_DECLS
-#ifdef __cplusplus
-# define __BEGIN_C_DECLS extern "C" {
-# define __END_C_DECLS }
-#else
-# define __BEGIN_C_DECLS /* empty */
-# define __END_C_DECLS /* empty */
-#endif
-/* End of C++ preparations */
-
-
-
-/* Actual header contants (the above were for the Pre-processor). */
-__BEGIN_C_DECLS /* From C++ preparations */
-
-
-
-#define GAL_MODE_LOW_QUANTILE 0.01f
-#define GAL_MODE_HIGH_QUANTILE 0.51f
-
-#define GAL_MODE_SYM_GOOD 0.2f
-#define GAL_MODE_LOW_QUANT_GOOD 0.02f
-
-#define GAL_MODE_SYMMETRICITY_LOW_QUANT 0.01f
-
-#define GAL_MODE_GOLDEN_RATIO 1.618034f
-#define GAL_MODE_TWO_TAKE_GOLDEN_RATIO 0.38197f
-
-#define GAL_MODE_MIRROR_IS_ABOVE_RESULT (size_t)(-1)
-
-struct gal_mode_params
-{
- float *sorted; /* Sorted array to be used. */
- size_t size; /* Number of elements in the sorted array. */
- size_t lowi; /* Lower quantile of interval. */
- size_t midi; /* Middle quantile of interval. */
- size_t midd; /* Middle index of inteval. */
- size_t highi; /* Higher quantile of interval. */
- float tolerance; /* Tolerance level to terminate search. */
- size_t numcheck; /* Number of pixels after mode to check. */
- size_t interval; /* Interval to check pixels. */
- float errorstdm; /* Multiple of standard deviation. */
-};
+size_t
+mirrormaxdiff(float *a, size_t size, size_t m,
+ size_t numcheck, size_t interval, size_t stdm);
void
-gal_mode_make_mirror_plots(float *sorted, size_t size, size_t mirrorindex,
- float min, float max, size_t numbins,
- char *histsname, char *cfpsname,
- float mirrorplotdist);
+modesymmetricity(float *a, size_t size, size_t mi, float errorstdm,
+ float *sym);
-float
-gal_mode_value_from_sym(float *sorted, size_t size, size_t modeindex,
- float sym);
+size_t
+modegoldenselection(struct gal_statistics_mode_params *mp);
void
-gal_mode_index_in_sorted(float *sorted, size_t size, float errorstdm,
- size_t *modeindex, float *modesym);
+makemirrored(float *in, size_t mi, float **outmirror, size_t *outsize);
-
-
-__END_C_DECLS /* From C++ preparations */
-
-#endif /* __GAL_MODE_H__ */
+#endif
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index 8bdee9b..5147675 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -267,6 +267,48 @@ gal_statistics_remove_outliers_flat_cdf(float *sorted,
size_t *outsize);
+
+
+/****************************************************************/
+/************* Mode ****************/
+/****************************************************************/
+#define GAL_STATISTICS_MODE_LOW_QUANTILE 0.01f
+#define GAL_STATISTICS_MODE_HIGH_QUANTILE 0.51f
+
+#define GAL_STATISTICS_MODE_SYM_GOOD 0.2f
+#define GAL_STATISTICS_MODE_LOW_QUANT_GOOD 0.02f
+
+#define GAL_STATISTICS_MODE_SYMMETRICITY_LOW_QUANT 0.01f
+
+struct gal_statistics_mode_params
+{
+ float *sorted; /* Sorted array to be used. */
+ size_t size; /* Number of elements in the sorted array. */
+ size_t lowi; /* Lower quantile of interval. */
+ size_t midi; /* Middle quantile of interval. */
+ size_t midd; /* Middle index of inteval. */
+ size_t highi; /* Higher quantile of interval. */
+ float tolerance; /* Tolerance level to terminate search. */
+ size_t numcheck; /* Number of pixels after mode to check. */
+ size_t interval; /* Interval to check pixels. */
+ float errorstdm; /* Multiple of standard deviation. */
+};
+
+void
+gal_statistics_mode_mirror_plots(float *sorted, size_t size, size_t
mirrorindex,
+ float min, float max, size_t numbins,
+ char *histsname, char *cfpsname,
+ float mirrorplotdist);
+
+float
+gal_statistics_mode_value_from_sym(float *sorted, size_t size, size_t
modeindex,
+ float sym);
+
+void
+gal_statistics_mode_index_in_sorted(float *sorted, size_t size, float
errorstdm,
+ size_t *modeindex, float *modesym);
+
+
__END_C_DECLS /* From C++ preparations */
#endif /* __GAL_STATISTICS_H__ */
diff --git a/lib/mesh.c b/lib/mesh.c
index aab457d..7144654 100644
--- a/lib/mesh.c
+++ b/lib/mesh.c
@@ -32,7 +32,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/fits.h>
#include <gnuastro/mesh.h>
-#include <gnuastro/mode.h>
#include <gnuastro/qsort.h>
#include <gnuastro/neighbors.h>
#include <gnuastro/linkedlist.h>
diff --git a/lib/mode.c b/lib/mode.c
index b5f67e1..d98b041 100644
--- a/lib/mode.c
+++ b/lib/mode.c
@@ -72,7 +72,7 @@ the functions below.
/* This is used for the plots, it will allocate an array and put the
mirrored array values in it. `mi` is the index the mirror is to
be placed on. */
-static void
+void
makemirrored(float *in, size_t mi, float **outmirror, size_t *outsize)
{
size_t i, size;
@@ -100,136 +100,6 @@ makemirrored(float *in, size_t mi, float **outmirror,
size_t *outsize)
-void
-gal_mode_make_mirror_plots(float *sorted, size_t size, size_t mirrorindex,
- float min, float max, size_t numbins,
- char *histsname, char *cfpsname,
- float mirrorplotdist)
-{
- FILE *fp;
- size_t i, msize;
- float *out, maxhist=-FLT_MAX, maxcfp, d;
- int normhist=0, maxhistone=0, normcfp=0;
- float *bins, *mirror, *actual, mf, onebinvalue=0.0f;
-
-
- /* Find the index of the mirror and make the mirror array: */
- mf=sorted[mirrorindex];
- gal_array_float_copy(sorted, size, &actual);
- makemirrored(sorted, mirrorindex, &mirror, &msize);
-
-
- /* Set the best range if asked for, such that the mirror is on the
- 1/3 of the image scale. */
- if(mirrorplotdist!=0.0f)
- {
- min=actual[gal_statistics_index_from_quantile(size, 0.001f)];
- max=mf+mirrorplotdist*(mf-min);
- }
-
-
- /* set the mirror pixel to have the value of zero.*/
- min-=mf;
- max-=mf;
- gal_array_fsum_const(actual, size, -1.0f*mf);
- gal_array_fsum_const(mirror, msize, -1.0f*mf);
-
-
- /* Allocate space for the 3 column array keeping the histograms:*/
- errno=0;
- out=malloc(numbins*3*sizeof *out);
- if(out==NULL)
- error(EXIT_FAILURE, errno, "%lu bytes for out in "
- "gal_mode_make_mirror_plots (mode.c)", numbins*3*sizeof *out);
-
-
- /* Define the bin sides: */
- gal_statistics_set_bins(actual, size, numbins, min,
- max, onebinvalue, 0, &bins);
-
-
- /* Find the histogram of the actual data and put it in out. Note
- that maxhistone=0, because here we want to use one value for both
- histograms so they are comparable. */
- gal_statistics_histogram(actual, size, bins, numbins,
- normhist, maxhistone);
- for(i=0;i<numbins;++i)
- if(bins[i*2+1]>maxhist) maxhist=bins[i*2+1];
- for(i=0;i<numbins;++i)
- { out[i*3]=bins[i*2]; out[i*3+1]=bins[i*2+1]/maxhist; bins[i*2+1]=0.0f;}
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
- d=(out[3]-out[0])/2;
-
-
- /* Find the histogram of the mirrired distribution and put it in
- out: */
- gal_statistics_histogram(mirror, msize, bins, numbins, normhist,
- maxhistone);
- for(i=0;i<numbins;++i)
- { out[i*3+2]=bins[i*2+1]/maxhist; bins[i*2+1]=0.0f;}
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
-
-
- /* Print out the histogram: */
- errno=0;
- fp=fopen(histsname, "w");
- if(fp==NULL)
- error(EXIT_FAILURE, errno, "could not open file %s", histsname);
- fprintf(fp, "# Histogram of actual and mirrored distributions.\n");
- fprintf(fp, "# Column 0: Value in the middle of this bin.\n");
- fprintf(fp, "# Column 1: Input data.\n");
- fprintf(fp, "# Column 2: Mirror distribution.\n");
- for(i=0;i<numbins;++i)
- fprintf(fp, "%-25.6f%-25.6f%-25.6f\n", out[i*3]+d,
- out[i*3+1], out[i*3+2]);
- fclose(fp);
-
-
-
-
- /* Find the cumulative frequency plots of the two distributions: */
- gal_statistics_cumulative_fp(actual, size, bins, numbins, normcfp);
- for(i=0;i<numbins;++i)
- { out[i*3+1]=bins[i*2+1]; bins[i*2+1]=0.0f; }
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
- gal_statistics_cumulative_fp(mirror, msize, bins, numbins, normcfp);
- for(i=0;i<numbins;++i)
- { out[i*3+2]=bins[i*2+1]; bins[i*2+1]=0.0f;}
- bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
-
-
- /* Since the two cumultiave frequency plots have to be on scale, see
- which one is larger and divided both CFPs by the size of the
- larger one. Then print the CFPs. */
- if(size>msize) maxcfp=size;
- else maxcfp=msize;
- errno=0;
- fp=fopen(cfpsname, "w");
- if(fp==NULL)
- error(EXIT_FAILURE, errno, "could not open file %s", cfpsname);
- fprintf(fp, "# Cumulative frequency plot (average index in bin) of\n"
- "# Actual and mirrored distributions.\n");
- fprintf(fp, "# Column 0: Value in the middle of this bin.\n");
- fprintf(fp, "# Column 1: Actual data.\n");
- fprintf(fp, "# Column 2: Mirror distribution.\n");
- for(i=0;i<numbins;++i)
- fprintf(fp, "%-25.6f%-25.6f%-25.6f\n", out[i*3],
- out[i*3+1]/maxcfp, out[i*3+2]/maxcfp);
- fclose(fp);
-
-
-
-
- free(out);
- free(bins);
- free(mirror);
- free(actual);
-}
-
-
-
-
-
@@ -275,7 +145,7 @@ gal_mode_make_mirror_plots(float *sorted, size_t size,
size_t mirrorindex,
`j`. So, if we keep the previous `j` in `prevj` then, all we have to
do is to start incrementing `j` from `prevj`. This will really help
in speeding up the job :-D. Only for the first element, `prevj=0`.*/
-static size_t
+size_t
mirrormaxdiff(float *a, size_t size, size_t m,
size_t numcheck, size_t interval, size_t stdm)
{
@@ -328,7 +198,7 @@ mirrormaxdiff(float *a, size_t size, size_t m,
result is not acceptable! */
if(i>j+errordiff)
{
- maxdiff=GAL_MODE_MIRROR_IS_ABOVE_RESULT;
+ maxdiff = MODE_MIRROR_IS_ABOVE_RESULT;
break;
}
absdiff = i>j ? i-j : j-i;
@@ -361,8 +231,8 @@ mirrormaxdiff(float *a, size_t size, size_t m,
We need a fourth point to be placed between. With this configuration,
the probing point is located at: */
-static size_t
-modegoldenselection(struct gal_mode_params *mp)
+size_t
+modegoldenselection(struct gal_statistics_mode_params *mp)
{
size_t di, dd;
/*static int counter=1;*/
@@ -373,9 +243,9 @@ modegoldenselection(struct gal_mode_params *mp)
/* Find the probing point in the larger interval. */
if(mp->highi-mp->midi > mp->midi-mp->lowi)
- di = mp->midi + GAL_MODE_TWO_TAKE_GOLDEN_RATIO
*(float)(mp->highi-mp->midi);
+ di = mp->midi + MODE_TWO_TAKE_GOLDEN_RATIO *(float)(mp->highi-mp->midi);
else
- di = mp->midi - GAL_MODE_TWO_TAKE_GOLDEN_RATIO *
(float)(mp->midi-mp->lowi);
+ di = mp->midi - MODE_TWO_TAKE_GOLDEN_RATIO * (float)(mp->midi-mp->lowi);
/* Since these are all indexs (and positive) we don't need an
absolute value, highi is also always larger than lowi! In some
@@ -408,7 +278,7 @@ modegoldenselection(struct gal_mode_params *mp)
section minimization is not going to give us what we want. So I have
added this modification to it. In such cases, we want the search to go
to the lower intervals.*/
- if(dd==GAL_MODE_MIRROR_IS_ABOVE_RESULT)
+ if(dd==MODE_MIRROR_IS_ABOVE_RESULT)
{
if(mp->midi < di)
{
@@ -473,7 +343,7 @@ modegoldenselection(struct gal_mode_params *mp)
the symmetricity of the mode can be quantified as: (b-m)/(m-a). For
a completly symmetric mode, this should be 1. Note that the search
for `b` only goes to the 95% of the distribution. */
-static void
+void
modesymmetricity(float *a, size_t size, size_t mi, float errorstdm,
float *sym)
{
@@ -484,7 +354,7 @@ modesymmetricity(float *a, size_t size, size_t mi, float
errorstdm,
errdiff=errorstdm*sqrt(mi);
topi = 2*mi>size-1 ? size-1 : 2*mi;
af=a[gal_statistics_index_from_quantile(2*mi+1,
- GAL_MODE_SYMMETRICITY_LOW_QUANT)];
+ GAL_STATISTICS_MODE_SYMMETRICITY_LOW_QUANT)];
/* This loop is very similar to that of mirrormaxdiff(). It will
find the index where the difference between the two cumulative
@@ -529,66 +399,3 @@ modesymmetricity(float *a, size_t size, size_t mi, float
errorstdm,
printf("symmetricity: %f\n", sym);
*/
}
-
-
-
-
-
-/* It happens that you have the symmetricity and you want the flux
- value at that point, this function will do that job. In practice,
- it just finds bf from the equation to calculate symmetricity in
- modesymmetricity. */
-float
-gal_mode_value_from_sym(float *sorted, size_t size, size_t modeindex,
- float sym)
-{
- float mf=sorted[modeindex];
- float af=
- sorted[gal_statistics_index_from_quantile(2*modeindex+1,
-
GAL_MODE_SYMMETRICITY_LOW_QUANT)];
- return sym*(mf-af)+mf;
-}
-
-
-
-
-
-/* Find the quantile of the mode of a sorted distribution. The return
- value is either 0 (not accurate) or 1 (accurate). Accuracy is
- defined based on the difference between the maximum and minimum
- maxdiffs that were found during the golden section search.
-
- A good mode will have:
-
- modequant=(float)(modeindex)/(float)size;
- modesym>GAL_MODE_SYM_GOOD && modequant>GAL_MODE_LOW_QUANT_GOOD
-*/
-void
-gal_mode_index_in_sorted(float *sorted, size_t size, float errorstdm,
- size_t *modeindex, float *modesym)
-{
- struct gal_mode_params mp;
-
- /* Initialize the gal_mode_params structure: */
- mp.size=size;
- mp.sorted=sorted;
- mp.tolerance=0.01;
- mp.numcheck=size/2;
- mp.errorstdm=errorstdm;
- if(mp.numcheck>1000)
- mp.interval=mp.numcheck/1000;
- else mp.interval=1;
- mp.lowi = gal_statistics_index_from_quantile(size,
- GAL_MODE_LOW_QUANTILE);
- mp.highi = gal_statistics_index_from_quantile(size,
- GAL_MODE_HIGH_QUANTILE);
- mp.midi = ((float)mp.highi
- +
GAL_MODE_GOLDEN_RATIO*(float)mp.lowi)/(1+GAL_MODE_GOLDEN_RATIO);
- mp.midd = mirrormaxdiff(mp.sorted, mp.size, mp.midi, mp.numcheck,
- mp.interval, mp.errorstdm);
-
- /* Do the golden section search and find the resulting
- symmetricity. */
- *modeindex=modegoldenselection(&mp);
- modesymmetricity(sorted, size, *modeindex, errorstdm, modesym);
-}
diff --git a/lib/statistics.c b/lib/statistics.c
index 21e0aa8..888472d 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -34,7 +34,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/array.h>
#include <gnuastro/statistics.h>
-
+#include "gnuastro/mode.h"
@@ -1382,3 +1382,214 @@ gal_statistics_remove_outliers_flat_cdf(float *sorted,
size_t *outsize)
free(slopes);
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/****************************************************************/
+/************* Mode calculation ****************/
+/****************************************************************/
+void
+gal_statistics_mode_mirror_plots(float *sorted, size_t size, size_t
mirrorindex,
+ float min, float max, size_t numbins,
+ char *histsname, char *cfpsname,
+ float mirrorplotdist)
+{
+ FILE *fp;
+ size_t i, msize;
+ float *out, maxhist=-FLT_MAX, maxcfp, d;
+ int normhist=0, maxhistone=0, normcfp=0;
+ float *bins, *mirror, *actual, mf, onebinvalue=0.0f;
+
+
+ /* Find the index of the mirror and make the mirror array: */
+ mf=sorted[mirrorindex];
+ gal_array_float_copy(sorted, size, &actual);
+ makemirrored(sorted, mirrorindex, &mirror, &msize);
+
+
+ /* Set the best range if asked for, such that the mirror is on the
+ 1/3 of the image scale. */
+ if(mirrorplotdist!=0.0f)
+ {
+ min=actual[gal_statistics_index_from_quantile(size, 0.001f)];
+ max=mf+mirrorplotdist*(mf-min);
+ }
+
+
+ /* set the mirror pixel to have the value of zero.*/
+ min-=mf;
+ max-=mf;
+ gal_array_fsum_const(actual, size, -1.0f*mf);
+ gal_array_fsum_const(mirror, msize, -1.0f*mf);
+
+
+ /* Allocate space for the 3 column array keeping the histograms:*/
+ errno=0;
+ out=malloc(numbins*3*sizeof *out);
+ if(out==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for out in "
+ "gal_mode_make_mirror_plots (mode.c)", numbins*3*sizeof *out);
+
+
+ /* Define the bin sides: */
+ gal_statistics_set_bins(actual, size, numbins, min,
+ max, onebinvalue, 0, &bins);
+
+
+ /* Find the histogram of the actual data and put it in out. Note
+ that maxhistone=0, because here we want to use one value for both
+ histograms so they are comparable. */
+ gal_statistics_histogram(actual, size, bins, numbins,
+ normhist, maxhistone);
+ for(i=0;i<numbins;++i)
+ if(bins[i*2+1]>maxhist) maxhist=bins[i*2+1];
+ for(i=0;i<numbins;++i)
+ { out[i*3]=bins[i*2]; out[i*3+1]=bins[i*2+1]/maxhist; bins[i*2+1]=0.0f;}
+ bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
+ d=(out[3]-out[0])/2;
+
+
+ /* Find the histogram of the mirrired distribution and put it in
+ out: */
+ gal_statistics_histogram(mirror, msize, bins, numbins, normhist,
+ maxhistone);
+ for(i=0;i<numbins;++i)
+ { out[i*3+2]=bins[i*2+1]/maxhist; bins[i*2+1]=0.0f;}
+ bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
+
+
+ /* Print out the histogram: */
+ errno=0;
+ fp=fopen(histsname, "w");
+ if(fp==NULL)
+ error(EXIT_FAILURE, errno, "could not open file %s", histsname);
+ fprintf(fp, "# Histogram of actual and mirrored distributions.\n");
+ fprintf(fp, "# Column 0: Value in the middle of this bin.\n");
+ fprintf(fp, "# Column 1: Input data.\n");
+ fprintf(fp, "# Column 2: Mirror distribution.\n");
+ for(i=0;i<numbins;++i)
+ fprintf(fp, "%-25.6f%-25.6f%-25.6f\n", out[i*3]+d,
+ out[i*3+1], out[i*3+2]);
+ fclose(fp);
+
+
+
+
+ /* Find the cumulative frequency plots of the two distributions: */
+ gal_statistics_cumulative_fp(actual, size, bins, numbins, normcfp);
+ for(i=0;i<numbins;++i)
+ { out[i*3+1]=bins[i*2+1]; bins[i*2+1]=0.0f; }
+ bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
+ gal_statistics_cumulative_fp(mirror, msize, bins, numbins, normcfp);
+ for(i=0;i<numbins;++i)
+ { out[i*3+2]=bins[i*2+1]; bins[i*2+1]=0.0f;}
+ bins[i*2+1]=0.0f; /* bins[] actually has numbins+1 elements. */
+
+
+ /* Since the two cumultiave frequency plots have to be on scale, see
+ which one is larger and divided both CFPs by the size of the
+ larger one. Then print the CFPs. */
+ if(size>msize) maxcfp=size;
+ else maxcfp=msize;
+ errno=0;
+ fp=fopen(cfpsname, "w");
+ if(fp==NULL)
+ error(EXIT_FAILURE, errno, "could not open file %s", cfpsname);
+ fprintf(fp, "# Cumulative frequency plot (average index in bin) of\n"
+ "# Actual and mirrored distributions.\n");
+ fprintf(fp, "# Column 0: Value in the middle of this bin.\n");
+ fprintf(fp, "# Column 1: Actual data.\n");
+ fprintf(fp, "# Column 2: Mirror distribution.\n");
+ for(i=0;i<numbins;++i)
+ fprintf(fp, "%-25.6f%-25.6f%-25.6f\n", out[i*3],
+ out[i*3+1]/maxcfp, out[i*3+2]/maxcfp);
+ fclose(fp);
+
+
+ /* Clean up. */
+ free(out);
+ free(bins);
+ free(mirror);
+ free(actual);
+}
+
+
+
+
+
+/* It happens that you have the symmetricity and you want the flux
+ value at that point, this function will do that job. In practice,
+ it just finds bf from the equation to calculate symmetricity in
+ modesymmetricity. */
+float
+gal_statistics_mode_value_from_sym(float *sorted, size_t size,
+ size_t modeindex, float sym)
+{
+ float mf=sorted[modeindex];
+ float af=
+ sorted[gal_statistics_index_from_quantile(2*modeindex+1,
+ GAL_STATISTICS_MODE_SYMMETRICITY_LOW_QUANT)];
+ return sym*(mf-af)+mf;
+}
+
+
+
+
+
+/* Find the quantile of the mode of a sorted distribution. The return
+ value is either 0 (not accurate) or 1 (accurate). Accuracy is
+ defined based on the difference between the maximum and minimum
+ maxdiffs that were found during the golden section search.
+
+ A good mode will have:
+
+ modequant=(float)(modeindex)/(float)size;
+ modesym>GAL_MODE_SYM_GOOD && modequant>GAL_MODE_LOW_QUANT_GOOD
+*/
+void
+gal_statistics_mode_index_in_sorted(float *sorted, size_t size, float
errorstdm,
+ size_t *modeindex, float *modesym)
+{
+ struct gal_statistics_mode_params mp;
+
+ /* Initialize the gal_mode_params structure: */
+ mp.size=size;
+ mp.sorted=sorted;
+ mp.tolerance=0.01;
+ mp.numcheck=size/2;
+ mp.errorstdm=errorstdm;
+ if(mp.numcheck>1000)
+ mp.interval=mp.numcheck/1000;
+ else mp.interval=1;
+ mp.lowi = gal_statistics_index_from_quantile(size,
+ GAL_STATISTICS_MODE_LOW_QUANTILE);
+ mp.highi = gal_statistics_index_from_quantile(size,
+ GAL_STATISTICS_MODE_HIGH_QUANTILE);
+ mp.midi = (((float)mp.highi
+ + MODE_GOLDEN_RATIO*(float)mp.lowi)
+ /(1+MODE_GOLDEN_RATIO));
+ mp.midd = mirrormaxdiff(mp.sorted, mp.size, mp.midi, mp.numcheck,
+ mp.interval, mp.errorstdm);
+
+ /* Do the golden section search and find the resulting
+ symmetricity. */
+ *modeindex=modegoldenselection(&mp);
+ modesymmetricity(sorted, size, *modeindex, errorstdm, modesym);
+}
diff --git a/src/imgstat/imgstat.c b/src/imgstat/imgstat.c
index e07b5b6..8caf171 100644
--- a/src/imgstat/imgstat.c
+++ b/src/imgstat/imgstat.c
@@ -30,7 +30,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <string.h>
#include <stdlib.h>
-#include <gnuastro/mode.h>
#include <gnuastro/timing.h>
#include <gnuastro/statistics.h>
@@ -60,27 +59,30 @@ reportsimplestats(struct imgstatparams *p)
"Standard deviation", std, "Median", med);
/* The mode: */
- gal_mode_index_in_sorted(p->sorted, p->size, p->mirrordist, &modeindex,
- &modesym);
+ gal_statistics_mode_index_in_sorted(p->sorted, p->size, p->mirrordist,
+ &modeindex, &modesym);
modequant=(float)(modeindex)/(float)(p->size);
/* Report the values: */
printf(" -- %-45s%.4f %g\n", "Mode (quantile, value)",
modequant, p->sorted[modeindex]);
- symvalue=gal_mode_value_from_sym(p->sorted, p->size, modeindex, modesym);
+ symvalue=gal_statistics_mode_value_from_sym(p->sorted, p->size,
+ modeindex, modesym);
printf(" -- %-45s%.4f %g\n", "Mode symmetricity and its cutoff"
" value", modesym, symvalue);
- if(modesym<GAL_MODE_SYM_GOOD)
+ if(modesym<GAL_STATISTICS_MODE_SYM_GOOD)
printf(" ## MODE SYMMETRICITY IS TOO LOW ##\n");
/* Save the mode histogram and cumulative frequency plot. Note
that if the histograms are to be built, then
mhistname!=NULL. */
if(p->mhistname)
- gal_mode_make_mirror_plots(p->sorted, p->size, modeindex, p->histmin,
- p->histmax, p->histnumbins, p->mhistname,
- p->mcfpname, p->histrangeformirror ? 0.0f
- : p->mirrorplotdist);
+ gal_statistics_mode_mirror_plots(p->sorted, p->size, modeindex,
+ p->histmin, p->histmax,
+ p->histnumbins, p->mhistname,
+ p->mcfpname, (p->histrangeformirror
+ ? 0.0f
+ : p->mirrorplotdist) );
}
@@ -266,13 +268,13 @@ imgstat(struct imgstatparams *p)
/* Make the mirror distribution if asked for: */
if(isnan(p->mirror)==0)
- gal_mode_make_mirror_plots(p->sorted, p->size,
- gal_statistics_index_from_quantile(p->size,
- p->mirror),
- p->histmin, p->histmax, p->histnumbins,
- p->mirrorhist, p->mirrorcfp,
- p->histrangeformirror ? 0.0f
- : p->mirrorplotdist);
+ gal_statistics_mode_mirror_plots(p->sorted, p->size,
+
gal_statistics_index_from_quantile(p->size,
+
p->mirror),
+ p->histmin, p->histmax, p->histnumbins,
+ p->mirrorhist, p->mirrorcfp,
+ (p->histrangeformirror ? 0.0f
+ : p->mirrorplotdist));
/* Print out the Sigma clippings: */
if(p->sigclip && p->cp.verb)
diff --git a/src/noisechisel/thresh.c b/src/noisechisel/thresh.c
index 48105ed..a18cf87 100644
--- a/src/noisechisel/thresh.c
+++ b/src/noisechisel/thresh.c
@@ -99,9 +99,9 @@ qthreshonmesh(void *inparam)
{
qsort(oneforall, num, sizeof *oneforall,
gal_qsort_float_increasing);
- gal_mode_index_in_sorted(oneforall, num, mirrordist,
- &modeindex, &modesym);
- if( modesym>GAL_MODE_SYM_GOOD
+ gal_statistics_mode_index_in_sorted(oneforall, num, mirrordist,
+ &modeindex, &modesym);
+ if( modesym>GAL_STATISTICS_MODE_SYM_GOOD
&& (float)modeindex/(float)num>minmodeq)
{
mp->garray1[ind] = oneforall[
diff --git a/src/subtractsky/subtractsky.c b/src/subtractsky/subtractsky.c
index 2c8402c..d443110 100644
--- a/src/subtractsky/subtractsky.c
+++ b/src/subtractsky/subtractsky.c
@@ -107,8 +107,10 @@ avestdonthread(void *inparam)
/* Do the desired operation on the mesh: */
qsort(sorted, num, sizeof *oneforall, gal_qsort_float_increasing);
- gal_mode_index_in_sorted(sorted, num, mirrordist, &modeindex, &modesym);
- if( modesym>GAL_MODE_SYM_GOOD && (float)modeindex/(float)num>minmodeq )
+ gal_statistics_mode_index_in_sorted(sorted, num, mirrordist,
+ &modeindex, &modesym);
+ if( modesym>GAL_STATISTICS_MODE_SYM_GOOD
+ && (float)modeindex/(float)num>minmodeq )
{
/* If cofa was defined, then oneforall was not sorted. */
if(cofa)
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 99f3c71: Mode functions moved to statistics.h,
Mohammad Akhlaghi <=