[gnuastro-commits] master faa4b03e: Library (pool.h): new library and arithmetic operator

[Mohammad Akhlaghi]

branch: master
commit faa4b03e781bc3843787948da93e9e279c679640
Author: Faezeh Bidjarchian <fbidjarchian@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (pool.h): new library and arithmetic operator
    
    Until now, there was no operator in Gnuastro for "pooling" in image
    analysis. The closest operation was scaling (through warp), but that didn't
    allow diffferent statistical criteria, only sum.
    
    With this commit, the Gnuastro library now has a 'pool.h' set of functions
    that can also be called within the Arithmetic program.
---
 NEWS                        |  19 +++
 bin/arithmetic/arithmetic.c |  39 ++++-
 doc/gnuastro.texi           | 174 ++++++++++++++++++++--
 lib/Makefile.am             |   2 +
 lib/arithmetic.c            |  89 ++++++++++++
 lib/gnuastro/arithmetic.h   |   9 ++
 lib/gnuastro/pool.h         |  66 +++++++++
 lib/pool.c                  | 344 ++++++++++++++++++++++++++++++++++++++++++++
 lib/statistics.c            |   2 +-
 9 files changed, 730 insertions(+), 14 deletions(-)

diff --git a/NEWS b/NEWS
index 0dbaf20c..1a6229cf 100644
--- a/NEWS
+++ b/NEWS
@@ -25,6 +25,25 @@ See the end of the file for license conditions.
     the different configuration files of different instances of the same
     program without overwriting them. See the example in the book.
 
+  Arithmetic:
+  - New operators:
+    - pool-min: Min-pooling to reduce the size of the input by calculating
+                the minimum of a the pixels within the pooling window. See
+                the new "Pooling opeators" section of the book for
+                more. This and the rest of the pooling operators introduced
+                here were implemented by Faezeh Bidjarchian.
+    - pool-max: Similar to 'pool-min' but using maximum.
+    - pool-sum: Similar to 'pool-min' but using sum.
+    - pool-mean: Similar to 'pool-min' but using mean.
+    - pool-median: Similar to 'pool-min' but using median.
+
+  Library:
+  -gal_pool_min: min-pooling function, see 'pool-min' above.
+  -gal_pool_max: max-pooling function, see 'pool-min' above.
+  -gal_pool_sum: sum-pooling function, see 'pool-min' above.
+  -gal_pool_mean: mean-pooling function, see 'pool-min' above.
+  -gal_pool_median: median-pooling function, see 'pool-min' above.
+
 ** Removed features
 
 ** Changed features
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 8aaf6263..7e3dfed0 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1507,17 +1507,47 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
                 PACKAGE_BUGREPORT, num_operands, operator_string);
         }
 
-      /* Meta operators (where only the information about the pixels is
-         relevant, not the pixel values themselves, like the 'index'
-         operator). For these (that only accept one operand), we should add
-         the first popped operand back on the stack, then  */
+      /* Operators with special attention. */
       switch(operator)
         {
+         /* Meta operators (where only the information about the pixels is
+            relevant, not the pixel values themselves, like the 'index'
+            operator). For these (that only accept one operand), we should
+            add the first popped operand back on the stack. */
         case GAL_ARITHMETIC_OP_INDEX:
         case GAL_ARITHMETIC_OP_COUNTER:
           operands_add(p, NULL, d1);
           d1=arithmetic_prepare_meta(d1, d2, d3);
           break;
+
+        /* Operators that need/modify the WCS. */
+        case GAL_ARITHMETIC_OP_POOLMAX:
+        case GAL_ARITHMETIC_OP_POOLMIN:
+        case GAL_ARITHMETIC_OP_POOLSUM:
+        case GAL_ARITHMETIC_OP_POOLMEAN:
+        case GAL_ARITHMETIC_OP_POOLMEDIAN:
+
+          /* Print warning in case of existing WCS. */
+          if(p->refdata.wcs)
+            {
+              if(p->cp.quiet==0)
+                error(EXIT_SUCCESS, 0, "WARNING: the WCS is not currently "
+                      "supported for the pooling operators (the output "
+                      "will not contain WCS). If it is necessary in your "
+                      "usage, please get in touch with us at '%s'. You "
+                      "can suppress this warning with the '--quiet' option",
+                      PACKAGE_BUGREPORT);
+              gal_wcs_free(p->refdata.wcs);
+              p->refdata.wcs=NULL;
+            }
+
+          /* In case you want to pass the WCS to the library function ...
+
+          d1->wcs=gal_wcs_copy(p->refdata.wcs);
+
+          ... uncomment the line above and comment the whole
+          'if(p->refdata.wcs)' check above (so the wcs is not freed).*/
+          break;
         }
 
       /* Run the arithmetic operation. Note that 'gal_arithmetic'
@@ -1527,6 +1557,7 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
          arguments will be ignored. */
       operands_add(p, NULL, gal_arithmetic(operator, p->cp.numthreads,
                                            flags, d1, d2, d3, d4));
+
     }
 
   /* No need to call the arithmetic library, call the proper wrappers
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 0c6a27c8..f1a3b80a 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -532,6 +532,7 @@ Arithmetic operators
 * Statistical operators::       Statistics of a single dataset (for example, 
mean).
 * Stacking operators::          Coadding or combining multiple datasets into 
one.
 * Filtering operators::         Smoothing a dataset through mixing pixel with 
neighbors.
+* Pooling operators::           Reducing size through statistics of pixels in 
window.
 * Interpolation operators::     Giving blank pixels a value.
 * Dimensionality changing operators::  Collapse or expand a dataset.
 * Conditional operators::       Select certain pixels within the dataset.
@@ -654,7 +655,7 @@ MakeCatalog
 
 Quantifying measurement limits
 
-* Standard deviation vs. error::  The std is not a measure of the error.
+* Standard deviation vs error::  The std is not a measure of the error.
 * Magnitude measurement error of each detection::  Error in measuring 
magnitude.
 * Surface brightness error of each detection::  Error in measuring the Surface 
brightness.
 * Completeness limit of each detection::  Possibility of detecting similar 
objects?
@@ -836,6 +837,7 @@ Gnuastro library
 * Binary datasets::             Datasets that can only have values of 0 or 1.
 * Labeled datasets::            Working with Segmented/labeled datasets.
 * Convolution functions::       Library functions to do convolution.
+* Pooling functions::           Reduce size of input by statistical methods.
 * Interpolation::               Interpolate (over blank values possibly).
 * Warp library::                Warp pixel grid to a new one.
 * Color functions::             Definitions and operations related to colors.
@@ -18057,6 +18059,7 @@ Reading NaN as a floating point number in Gnuastro is 
not case-sensitive.
 * Statistical operators::       Statistics of a single dataset (for example, 
mean).
 * Stacking operators::          Coadding or combining multiple datasets into 
one.
 * Filtering operators::         Smoothing a dataset through mixing pixel with 
neighbors.
+* Pooling operators::           Reducing size through statistics of pixels in 
window.
 * Interpolation operators::     Giving blank pixels a value.
 * Dimensionality changing operators::  Collapse or expand a dataset.
 * Conditional operators::       Select certain pixels within the dataset.
@@ -18759,7 +18762,7 @@ astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-std
 @end example
 @end table
 
-@node Filtering operators, Interpolation operators, Stacking operators, 
Arithmetic operators
+@node Filtering operators, Pooling operators, Stacking operators, Arithmetic 
operators
 @subsubsection Filtering (smoothing) operators
 
 Image filtering is commonly used for smoothing: every pixel value in the 
output image is created by applying a certain statistic to the pixels in its 
vicinity.
@@ -18831,7 +18834,126 @@ Apply a @mymath{\sigma}-clipped median filtering onto 
the input dataset.
 This operator and its necessary operands are almost identical to 
@code{filter-sigclip-mean}, except that after @mymath{\sigma}-clipping, the 
median value (which is less affected by outliers than the mean) is added back 
to the stack.
 @end table
 
-@node Interpolation operators, Dimensionality changing operators, Filtering 
operators, Arithmetic operators
+@node Pooling operators, Interpolation operators, Filtering operators, 
Arithmetic operators
+@subsubsection Pooling operators
+
+@cindex Pooling
+@cindex Convolutional Neural Networks
+Pooling is one way of reducing the complexity of the input image by grouping 
multiple input pixels into one output pixel (using any statistical measure).
+As a result, the output image has fewer pixels (less complexity).
+In Computer Vision, Pooling is commonly used in 
@url{https://en.wikipedia.org/wiki/Convolutional_neural_network, Convolutional 
Neural Networks} (CNNs).
+
+In pooling, the inputs are an image (e.g., a FITS file) and a square window 
pixel size (known as a pooling window).
+The window has to be smaller than the input's number of pixels in both 
dimensions and its width is called the pool size.
+This window slides over all pixels in the input from the top-left corner to 
the bottom-right corner (covering each input pixel only once).
+Currently the ``stride'' (or spacing between the windows as they slide over 
the input) is equal to the window-size in Arithmetic.
+In other words, in pooling, the separate ``windows'' do not overlap with each 
other on the input.
+Therefore there are two major differences with @ref{Spatial domain 
convolution} or @ref{Filtering operators}, but pooling has some similarities to 
the @ref{Warp}.
+@itemize
+@item
+In convolution or filtering the input and output sizes are the same.
+However, the output of pooling has fewer pixels.
+@item
+In convolution or filters, the kernels slide of the input in a pixel-by-pixel 
manner.
+As a result, the same pixel's value will be used in many of the output pixels.
+However, in pooling each input pixel is only used in a single output pixel.
+@item
+Special cases of Warping an image are similar to pooling.
+For example calling @code{pool-sum} with pool size of 2 will give the same 
pixel values (except the outer edges) as giving the same image to 
@command{astwarp} with @option{--scale=1/2 --centeroncorner}.
+However, warping will only provide the sum of the input pixels, there is no 
easy way to generically define something like @code{pool-max} in Warp (which is 
far more general than pooling).
+Also, due to its generic features (for example for non-linear warps), Warp is 
slower than the @code{pool-max} that is introduced here.
+@end itemize
+
+@cartouche
+@noindent
+@strong{No WCS in output:} As of Gnuastro @value{VERSION}, the output of 
pooling will not contain WCS information (primarily due to a lack of time by 
developers).
+Please inform us of your interest in having it, by contacting us at 
@command{bug-gnuastro@@gnu.org}.
+If you need @code{pool-sum}, you can use @ref{Warp} (which also modifies the 
WCS, see note above).
+@end cartouche
+
+If the width or height of input is not divisible by the pool size, the pool 
window will go beyond the input pixel grid.
+In this case, the window pixels that do not overlap with the input are given a 
blank value (and thus ignored in the calculation of the desired statistical 
operation).
+
+The simple ASCII figure below shows the pooling operation where the input is a 
@mymath{3\times3} pixel image with a pool size of 2 pixels.
+In the center of the second row, you see the intermediate input matrix that 
highlights how the input and output pixels relate with each other.
+Since the input is @mymath{3\times3} and we have a pool size of 2, as 
mentioned above blank pseudo-pixels are added with a value of @code{B} (for 
blank).
+
+@example
+        Pool window:                             Input:
+        +-----------+                           +-------------+
+        |     |     |                           | 10   12   9 |
+        | _ _ | _ _ |___________________________| 31   4    1 |
+        |     |     |       ||          ||      | 16   5    8 |
+        |     |     |       ||          ||      +-------------+
+        +-----------+       ||          ||
+    The pooling window 2*2  ||          ||
+           stride 2         \/          \/
+                        +---------------------+
+                        |/ 10   12\|/ 9    B \|
+                        |          |          |
+  +-------+  pool-min   |\ 31   4 /|\ 1    B /|   pool-max  +-------+
+  | 4   1 |   /------   |---------------------|   ------\   |31   9 |
+  | 5   8 |   \------   |/ 16   5 \|/ 8    B \|   ------/   |16   8 |
+  +-------+             |          |          |             +-------+
+                        |\ B    B /.\ B    B /|
+                        +---------------------+
+@end example
+
+The choice of the statistic to use depends on the specific use case, the 
characteristics of the input data, and the desired output.
+Each statistic has its advantages and disadvantages and the choice of which to 
use should be informed by the specific needs of the problem at hand.
+Below, the various pool operators of arithmetic are listed:
+
+@table @command
+
+@item pool-max
+Apply max-pooling on the input dataset.
+This operator takes two operands: the first popped operand is the width of the 
square pooling window (which should be a single integer), and the second should 
be the input image.
+Within the pooling window, this operator will place the largest value in the 
output pixel (any blank pixels will be ignored).
+
+See the ASCII diagram above for a demonstration of how max-pooling works.
+Here is an example of using this operator:
+
+@example
+$ astarithmetic image.fits 2 pool-max
+@end example
+
+Max-pooling retains the largest value of the input window in the output, so 
the returned image is sharper where you have strong signal-to-noise ratio and 
more noisy in regions with no significant signal (only noise).
+It is therefore useful when the background of the image is dark and we are 
interested in only the highest signal-to-noise ratio regions of the image.
+
+@item pool-min
+Apply min-pooling on the input dataset.
+This operator takes two operands: the first popped operand is the width of the 
square pooling window (which should be a single integer), and the second should 
be the input image.
+Except the used statistical measurement, this operator is similar to 
@code{pool-max}, see the description there for more.
+
+Min-pooling is mostly used when the image has a high signal-to-noise ratio and 
a light background: min-pooling will select darker (lower-valued) pixels.
+For low signal-to-noise regions, this operator will increase the noise level 
(similar to the maximum, the scatter in the minimum is very strong).
+
+@item pool-sum
+Apply sum-pooling to the input dataset.
+This operator takes two operands: the first popped operand is the width of the 
square pooling window (which should be a single integer), and the second should 
be the input image.
+Except the used statistical measurement, this operator is similar to 
@code{pool-max}, see the description there for more.
+
+Sum-pooling will increase the signal-to-noise ratio at the cost of having a 
smoother output (less resolution).
+
+@item pool-mean
+Apply mean pooling on the input dataset.
+This operator takes two operands: the first popped operand is the width of the 
square pooling window (which should be a single integer), and the second should 
be the input image.
+Except the used statistical measurement, this operator is similar to 
@code{pool-max}, see the description there for more.
+
+The mean pooling method smooths out the image and hence the sharp features may 
not be identified when this pooling method is used.
+This therefore preserves more information than max-pooling, but may also 
reduces the effect of the most prominent pixels.
+Mean is often used where a more accurate representation of the input is 
required.
+
+@item pool-median
+Apply median pooling on the input dataset.
+This operator takes two operands: the first popped operand is the width of the 
square pooling window (which should be a single integer), and the second should 
be the input image.
+Except the used statistical measurement, this operator is similar to 
@code{pool-max}, see the description there for more.
+
+In general, the mean is mathematically easier to interpret and more 
susceptible to outliers, while the median outputs as being less subject to the 
influence of outliers compared to the mean so we have a smoother image.
+This is therefore better for low signal-to-ratio (noisy) features and extended 
features (where you don't want a single high or low valued pixel to affect the 
output).
+@end table
+
+@node Interpolation operators, Dimensionality changing operators, Pooling 
operators, Arithmetic operators
 @subsubsection Interpolation operators
 
 Interpolation is the process of removing blank pixels from a dataset (by 
giving them a value based on the non-blank neighbors).
@@ -25765,7 +25887,7 @@ In astronomy, it is common to use the magnitude (a 
unit-less scale) and physical
 Therefore the measurements discussed here are commonly used in units of 
magnitudes.
 
 @menu
-* Standard deviation vs. error::  The std is not a measure of the error.
+* Standard deviation vs error::  The std is not a measure of the error.
 * Magnitude measurement error of each detection::  Error in measuring 
magnitude.
 * Surface brightness error of each detection::  Error in measuring the Surface 
brightness.
 * Completeness limit of each detection::  Possibility of detecting similar 
objects?
@@ -25775,8 +25897,8 @@ Therefore the measurements discussed here are commonly 
used in units of magnitud
 * Upper limit magnitude of image::  Measure the noise-level for a certain 
aperture.
 @end menu
 
-@node Standard deviation vs. error, Magnitude measurement error of each 
detection, Quantifying measurement limits, Quantifying measurement limits
-@subsubsection Standard deviation vs. error
+@node Standard deviation vs error, Magnitude measurement error of each 
detection, Quantifying measurement limits, Quantifying measurement limits
+@subsubsection Standard deviation vs error
 The error and the standard deviation are sometimes confused with each other.
 Therefore, before continuing with the various measurement limits below, let's 
review these two fundamental concepts.
 Instead of going into the theoretical defitions of the two (which you can see 
in their resepctive Wikipedia pages), we'll discuss the concepts in a hands-on 
and practical way here.
@@ -25943,7 +26065,7 @@ If this image is not available, it is possible to use 
the  @code{SKY_STD} extens
 For more see @ref{NoiseChisel output}.
 @end table
 
-@node Magnitude measurement error of each detection, Surface brightness error 
of each detection, Standard deviation vs. error, Quantifying measurement limits
+@node Magnitude measurement error of each detection, Surface brightness error 
of each detection, Standard deviation vs error, Quantifying measurement limits
 @subsubsection Magnitude measurement error of each detection
 The raw error in measuring the magnitude is only meaningful when the object's 
magnitude is brighter than the upper-limit magnitude (see below).
 As discussed in @ref{Brightness flux magnitude}, the magnitude (@mymath{M}) of 
an object with brightness @mymath{B} and zero point magnitude @mymath{z} can be 
written as:
@@ -33243,6 +33365,7 @@ If you use the Info version of this manual (see 
@ref{Info}), you do not have to
 * Binary datasets::             Datasets that can only have values of 0 or 1.
 * Labeled datasets::            Working with Segmented/labeled datasets.
 * Convolution functions::       Library functions to do convolution.
+* Pooling functions::           Reduce size of input by statistical methods.
 * Interpolation::               Interpolate (over blank values possibly).
 * Warp library::                Warp pixel grid to a new one.
 * Color functions::             Definitions and operations related to colors.
@@ -40263,7 +40386,7 @@ For example, in a 2D dataset, a connectivity of 
@code{1} and @code{2} correspond
 @end deftypefun
 
 
-@node Convolution functions, Interpolation, Labeled datasets, Gnuastro library
+@node Convolution functions, Pooling functions, Labeled datasets, Gnuastro 
library
 @subsection Convolution functions (@file{convolve.h})
 
 Convolution is a very common operation during data analysis and is
@@ -40317,7 +40440,40 @@ channel borders. Other pixels in the image will not be 
touched. Hence, it
 is much faster.
 @end deftypefun
 
-@node Interpolation, Warp library, Convolution functions, Gnuastro library
+@node Pooling functions, Interpolation, Convolution functions, Gnuastro library
+@subsection Pooling functions (@file{pool.h})
+
+Pooling is the process of reducing the complexity of the input image (its size 
and variation of pixel values).
+Its underlying concepts and an analysis of its usefuless is fully descibed in 
@ref{Pooling operators}.
+The following functions are available pooling in Gnuastro.
+Just note that unlike the Arithmetic operators, the output of these functions 
should contain a correct WCS in their output.
+
+@deftypefun {gal_data_t *} gal_pool_max (gal_data_t @code{*input}, size_t 
@code{psize}, size_t @code{numthreads})
+Return the max-pool of @code{input}, assuming a pool size of @code{psize} 
pixels.
+The number of threads to use can be set with @code{numthreads}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_pool_min (gal_data_t @code{*input}, size_t 
@code{psize}, size_t @code{numthreads})
+Return the min-pool of @code{input}, assuming a pool size of @code{psize} 
pixels.
+The number of threads to use can be set with @code{numthreads}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_pool_sum (gal_data_t @code{*input}, size_t 
@code{psize}, size_t @code{numthreads})
+Return the sum-pool of @code{input}, assuming a pool size of @code{psize} 
pixels.
+The number of threads to use can be set with @code{numthreads}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_pool_mean (gal_data_t @code{*input}, size_t 
@code{psize}, size_t @code{numthreads})
+Return the mean-pool of @code{input}, assuming a pool size of @code{psize} 
pixels.
+The number of threads to use can be set with @code{numthreads}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_pool_median (gal_data_t @code{*input}, size_t 
@code{psize}, size_t @code{numthreads})
+Return the median-pool of @code{input}, assuming a pool size of @code{psize} 
pixels.
+The number of threads to use can be set with @code{numthreads}.
+@end deftypefun
+
+@node Interpolation, Warp library, Pooling functions, Gnuastro library
 @subsection Interpolation (@file{interpolate.h})
 
 @cindex Sky line
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 5e1828c8..05635be1 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -130,6 +130,7 @@ libgnuastro_la_SOURCES = \
   permutation.c \
   pointer.c \
   polygon.c \
+  pool.c \
   qsort.c \
   dimension.c \
   speclines.c \
@@ -182,6 +183,7 @@ pkginclude_HEADERS = gnuastro/config.h \
   $(headersdir)/permutation.h \
   $(headersdir)/pointer.h \
   $(headersdir)/polygon.h \
+  $(headersdir)/pool.h \
   $(headersdir)/qsort.h \
   $(headersdir)/speclines.h \
   $(headersdir)/statistics.h \
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index af4edd37..abc2ebf1 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -5,6 +5,7 @@ This is part of GNU Astronomy Utilities (Gnuastro) package.
 Corresponding author:
      Mohammad Akhlaghi <mohammad@akhlaghi.org>
 Contributing author(s):
+     Faezeh Bidjarchian <fbidjarchian@gmail.com>
      Pedram Ashofteh-Ardakani <pedramardakani@pm.me>
 Copyright (C) 2016-2023 Free Software Foundation, Inc.
 
@@ -39,6 +40,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/box.h>
 #include <gnuastro/wcs.h>
 #include <gnuastro/list.h>
+#include <gnuastro/pool.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/units.h>
 #include <gnuastro/qsort.h>
@@ -2947,6 +2949,64 @@ arithmetic_constant(gal_data_t *input, gal_data_t 
*constant, int operator,
 
 
 
+static gal_data_t *
+arithmetic_pool(gal_data_t *input, gal_data_t *psize, int operator,
+                size_t numthreads, int flags)
+{
+  gal_data_t *out=NULL;
+  size_t *psizearr;
+
+  /* Sanity check. */
+  if(psize->size!=1)
+    error(EXIT_FAILURE, 0, "%s: the pooling operand should only contain "
+          "a single element.", __func__);
+
+  /* This function is only for a integer operand, so make sure the user
+     has given an integer type. */
+  if(psize->type==GAL_TYPE_FLOAT32 || psize->type==GAL_TYPE_FLOAT64)
+    error(EXIT_FAILURE, 0, "lengths of pooling along dimensions "
+          "must be integer values, not floats. The given length "
+          "along dimension %zu is a float.", psize->dsize[0]);
+
+  /* Convert psize to size_t. */
+  psize=gal_data_copy_to_new_type_free(psize, GAL_TYPE_SIZE_T);
+  psizearr=psize->array;
+
+  /* Call the separate functions. */
+  switch(operator)
+    {
+    case GAL_ARITHMETIC_OP_POOLMAX:
+      out=gal_pool_max(input, psizearr[0], numthreads);
+      break;
+    case GAL_ARITHMETIC_OP_POOLMIN:
+      out=gal_pool_min(input, psizearr[0], numthreads);
+      break;
+    case GAL_ARITHMETIC_OP_POOLSUM:
+      out=gal_pool_sum(input, psizearr[0], numthreads);
+      break;
+    case GAL_ARITHMETIC_OP_POOLMEAN:
+      out=gal_pool_mean(input, psizearr[0], numthreads);
+      break;
+    case GAL_ARITHMETIC_OP_POOLMEDIAN:
+      out=gal_pool_median(input, psizearr[0], numthreads);
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
+            "'%s' to fix the problem. The value '%d' to "
+            "'operator' is not recognized", __func__,
+            PACKAGE_BUGREPORT, operator);
+    }
+
+  /* Clean up and return. */
+  if(out!=input && (flags & GAL_ARITHMETIC_FLAG_FREE))
+    gal_data_free(input);
+  return out;
+}
+
+
+
+
+
 gal_data_t *
 gal_arithmetic_load_col(char *str, int searchin, int ignorecase,
                         size_t minmapsize, int quietmmap)
@@ -3344,6 +3404,18 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
   else if (!strcmp(string, "constant"))
     { op=GAL_ARITHMETIC_OP_CONSTANT;          *num_operands=2;  }
 
+  /* Pooling operators. */
+  else if (!strcmp(string, "pool-max"))
+    { op=GAL_ARITHMETIC_OP_POOLMAX;           *num_operands=2;  }
+  else if (!strcmp(string, "pool-min"))
+    { op=GAL_ARITHMETIC_OP_POOLMIN;           *num_operands=2;  }
+  else if (!strcmp(string, "pool-sum"))
+    { op=GAL_ARITHMETIC_OP_POOLSUM;           *num_operands=2;  }
+  else if (!strcmp(string, "pool-mean"))
+    { op=GAL_ARITHMETIC_OP_POOLMEAN;          *num_operands=2;  }
+  else if (!strcmp(string, "pool-median"))
+    { op=GAL_ARITHMETIC_OP_POOLMEDIAN;        *num_operands=2;  }
+
   /* Operator not defined. */
   else
     { op=GAL_ARITHMETIC_OP_INVALID; *num_operands=GAL_BLANK_INT; }
@@ -3493,6 +3565,12 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_COUNTERONLY:     return "counteronly";
     case GAL_ARITHMETIC_OP_MAKENEW:         return "makenew";
 
+    case GAL_ARITHMETIC_OP_POOLMAX:         return "pool-max";
+    case GAL_ARITHMETIC_OP_POOLMIN:         return "pool-min";
+    case GAL_ARITHMETIC_OP_POOLSUM:         return "pool-sum";
+    case GAL_ARITHMETIC_OP_POOLMEAN:        return "pool-mean";
+    case GAL_ARITHMETIC_OP_POOLMEDIAN:      return "pool-median";
+
     default:                                return NULL;
     }
   return NULL;
@@ -3768,6 +3846,17 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       out=arithmetic_constant(d1, d2, operator, flags);
       break;
 
+      /* Pooling operators. */
+    case GAL_ARITHMETIC_OP_POOLMAX:
+    case GAL_ARITHMETIC_OP_POOLMIN:
+    case GAL_ARITHMETIC_OP_POOLSUM:
+    case GAL_ARITHMETIC_OP_POOLMEAN:
+    case GAL_ARITHMETIC_OP_POOLMEDIAN:
+      d1 = va_arg(va, gal_data_t *);
+      d2 = va_arg(va, gal_data_t *);
+      out=arithmetic_pool(d1, d2, operator, numthreads, flags);
+      break;
+
     /* When the operator is not recognized. */
     default:
       error(EXIT_FAILURE, 0, "%s: the argument '%d' could not be "
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 94794131..fdaa4f49 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -5,6 +5,7 @@ This is part of GNU Astronomy Utilities (Gnuastro) package.
 Original author:
      Mohammad Akhlaghi <mohammad@akhlaghi.org>
 Contributing author(s):
+     Faezeh Bidjarchian <fbidjarchian@gmail.com>
 Copyright (C) 2017-2023 Free Software Foundation, Inc.
 
 Gnuastro is free software: you can redistribute it and/or modify it
@@ -220,6 +221,14 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_COUNTERONLY,  /* New with the index (counting from 1). */
   GAL_ARITHMETIC_OP_SWAP,         /* Swap the top two operands.            */
 
+  /* Pooling operators. */
+  GAL_ARITHMETIC_OP_POOLMAX,      /* The pool-max of desired pixels.       */
+  GAL_ARITHMETIC_OP_POOLMIN,      /* The pool-min of desired pixels.       */
+  GAL_ARITHMETIC_OP_POOLSUM,      /* The pool-sum of desired pixels.       */
+  GAL_ARITHMETIC_OP_POOLMEAN,     /* The pool-mean of desired pixels.      */
+  GAL_ARITHMETIC_OP_POOLMEDIAN,   /* The pool-median of desired pixels.    */
+
+  /* Counter for number of operators. */
   GAL_ARITHMETIC_OP_LAST_CODE,    /* Last code of the library operands.    */
 };
 
diff --git a/lib/gnuastro/pool.h b/lib/gnuastro/pool.h
new file mode 100644
index 00000000..c574ab58
--- /dev/null
+++ b/lib/gnuastro/pool.h
@@ -0,0 +1,66 @@
+/*********************************************************************
+Pool - Pool input data and create a new dataset.
+Pool is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Faezeh Bidjarchian <fbidjarchian@gmail.com>
+Contributing author(s):
+     Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Copyright (C) 2023-2023 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 __GAL_POOL_H__
+#define __GAL_POOL_H__
+
+/* Include other headers if necessary here. Note that other header files
+   must be included before the C++ preparations below */
+#include <gnuastro/data.h>
+
+/* 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 */
+
+
+
+gal_data_t *
+gal_pool_max(gal_data_t *input, size_t psize, size_t numthreads);
+
+gal_data_t *
+gal_pool_min(gal_data_t *input, size_t psize, size_t numthreads);
+
+gal_data_t *
+gal_pool_sum(gal_data_t *input, size_t psize, size_t numthreads);
+
+gal_data_t *
+gal_pool_mean(gal_data_t *input, size_t psize, size_t numthreads);
+
+gal_data_t *
+gal_pool_median(gal_data_t *input, size_t psize, size_t numthreads);
+
+
+__END_C_DECLS    /* From C++ preparations */
+
+#endif
diff --git a/lib/pool.c b/lib/pool.c
new file mode 100644
index 00000000..dcf754f3
--- /dev/null
+++ b/lib/pool.c
@@ -0,0 +1,344 @@
+/*********************************************************************
+Pool - Pool input data and create a new dataset.
+Pool is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Faezeh Bidjarchian <fbidjarchian@gmail.com>
+Contributing author(s):
+     Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Copyright (C) 2023-2023 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 <string.h>
+
+#include <gnuastro/wcs.h>
+#include <gnuastro/type.h>
+#include <gnuastro/fits.h>
+#include <gnuastro/pool.h>
+#include <gnuastro/pointer.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
+
+#include <gnuastro-internal/checkset.h>
+
+
+
+
+
+/* Identifiers for each operator. */
+enum pool_operators
+{
+  POOL_MAX,      /* Maximum the desired pixels.  */
+  POOL_MIN,      /* Minimum the desired pixels.  */
+  POOL_SUM,      /* Sum of desired pixels.       */
+  POOL_MEAN,     /* Mean the desired pixel.      */
+  POOL_MEDIAN,   /* Median the desired pixels.   */
+};
+
+
+
+
+
+/**********************************************************************/
+/****************             Pooling                 *****************/
+/**********************************************************************/
+#define POOLING_DIM 10
+
+/* Main input/output parameters. */
+struct pooling
+{
+  int       operator;     /* The type of pooling.             */
+  size_t    poolsize;     /* The size of pooling.             */
+  size_t      *osize;     /* The output size.                 */
+  gal_data_t  *input;     /* Dataset to print values of.      */
+  gal_data_t    *out;     /* Output data structure.           */
+};
+
+
+
+
+
+/* Do the pooling on each thread.
+
+   Current assumptions:
+
+   - The size of pooling can be every single number (the pooling window is
+     a square).
+   - The width and height of the input are not necessarily divisible
+     by the size of the pooling. In other words, the image can be both
+     square and rectangular.
+   - We apply the pooling to our input with a stride of poolsize. So,
+     for example the stride is 2 when poolsize is equal to 2. */
+static void *
+pool_type_on_thread(void *in_prm)
+{
+  struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+  struct pooling *pp=(struct pooling *)tprm->params;
+  gal_data_t *input=pp->input;
+
+  size_t pools= pp->poolsize;
+  size_t ndim=input->ndim;
+  size_t w=input->dsize[1], wr;
+  size_t h=input->dsize[0], hr;
+  gal_data_t *statv=NULL, *result=NULL;
+  size_t i, a, b, oind, iind, vc, numpixs, coord[POOLING_DIM], index;
+
+  /* All number of pixels that we selected each time par the pooling. */
+  numpixs=pools*pools;
+
+  /* Allocated memory for selected pixels of input to feed into the
+     statisitical operation. */
+  statv=gal_data_alloc(NULL, input->type, 1, &numpixs, NULL, 0,
+                       input->minmapsize, input->quietmmap,
+                       NULL, NULL, NULL);
+
+  /* Go over all the pixels that were assigned to this thread. */
+  for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
+    {
+      /* For easy reading, put the index in 'ind'. */
+      oind=tprm->indexs[i];
+
+      /* Get the coordinate of the pixel. */
+      gal_dimension_index_to_coord(oind, ndim, pp->osize, coord);
+
+      /* Convert the pixel coordinate to the desired pixel that we must
+         select to set the pooling's starting pointer. */
+      iind=pools*w*coord[0]+pools*coord[1];
+      hr=iind%h;
+      wr=iind%w;
+
+      /* In some cases, the pooling window doesn't cover a whole squared
+         window and only has maybe one pixel. So since we initialize the
+         statv by Null (=0 in C), some of these values fill with input
+         values and others remain zero. So these zeros will affect the
+         outputs.  Therefore we initialize the statv by blank value. */
+      gal_blank_initialize(statv);
+
+      /* Set the sorted and blank check flags to 0 so the statistical
+         operators re-check everytime for sorted or blank elements. */
+      statv->flag=0;
+
+      /* Depending on the size of pooling, we condider a set of pixels
+        and fill the temporary 'values' array. Then we do the required
+         operation on them. */
+      vc=0;
+      for(a=0;a<pools && hr+a<=input->dsize[0];++a)
+       for(b=0;b<pools && wr+b<input->dsize[1];++b)
+         {
+           index=iind + a*w + b;
+           if (index<input->size)
+             {
+                memcpy(gal_pointer_increment(statv->array, vc,
+                                             statv->type),
+                       gal_pointer_increment(input->array, index,
+                                             input->type),
+                       gal_type_sizeof(input->type));
+                vc++;
+             }
+         }
+
+      /* Do the necessary calculation. */
+      switch(pp->operator)
+        {
+        case POOL_MAX:    result=gal_statistics_maximum(statv);   break;
+        case POOL_MIN:    result=gal_statistics_minimum(statv);   break;
+        case POOL_SUM:    result=gal_statistics_sum(statv);       break;
+        case POOL_MEAN:   result=gal_statistics_mean(statv);      break;
+        case POOL_MEDIAN: result=gal_statistics_median(statv, 1); break;
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+                "to fix the problem. 'operator' code %d is not "
+                "recognized.", PACKAGE_BUGREPORT, __func__,
+                pp->operator);
+        }
+
+      /* Make sure the output array type and the type of result are the
+         same. */
+      if(result->type!=pp->out->type)
+        result=gal_data_copy_to_new_type_free(result, pp->out->type);
+
+      /* Copy the result into the output array. */
+      memcpy(gal_pointer_increment(pp->out->array, oind, pp->out->type),
+             result->array, gal_type_sizeof(pp->out->type));
+
+      /* Clean up. */
+      gal_data_free(result);
+    }
+
+  /* Clean up. */
+  gal_data_free(statv);
+
+  /* Wait for all the other threads to finish, then return. */
+  if(tprm->b) pthread_barrier_wait(tprm->b);
+  return NULL;
+}
+
+
+
+
+
+static gal_data_t *
+pool_generic(gal_data_t *input, size_t psize, int operator, size_t numthreads)
+{
+  struct pooling pp={0};
+
+  int otype=GAL_TYPE_INVALID;
+  size_t outndim=input->ndim, i, r, diff;
+
+  /* Print a warning if the psize has a wrong value. It happens when the
+     user writes a negative value for the poolsize. */
+  if(psize>(size_t)(-1)/2 || psize==0)
+    error(EXIT_FAILURE, 0, "the value of poolsize must be positive, and "
+          "non zero)");
+
+  /* Make sure the given poolsize is lower than the input's width or
+     hight. */
+  if(psize>input->dsize[0] || psize>input->dsize[1])
+    error(EXIT_FAILURE, 0, "%s: the pool size along dimension must be "
+         "greater than the input's width or hight in that dimension",
+         __func__);
+
+  /* Set the pointers in the structure of the parameter. */
+  pp.input=input;
+  pp.poolsize=psize;
+  pp.operator=operator;
+
+  if(pp.input->size==1) pp.out=pp.input;
+  else
+    {
+      /* Resize output when calculating pooling on it and the remainder
+         is not zero. So we must calculate the pooling one more time for
+         the remaining pixels. */
+      pp.osize=gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
+                                    __func__, "osize");
+
+      /* if the width is not divisible by the poolsize, we have to add one
+         dimension to the output dimension for these remaining pixels. */
+      for(i=0;i<input->ndim;++i)
+        {
+          r=(input->dsize[i])%psize;
+          diff=((r==0)?0:1);
+          pp.osize[i]=(input->dsize[i]/psize)+diff;
+        }
+
+      /* Set the type of the output dataset. */
+      switch(pp.operator)
+        {
+        case POOL_MAX:
+        case POOL_MIN:
+        case POOL_MEDIAN:
+          otype=pp.input->type;
+          break;
+
+        case POOL_SUM:
+        case POOL_MEAN:
+          otype=GAL_TYPE_FLOAT64;
+          break;
+
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix "
+                "the problem. The 'operator' code %d is not recognized",
+                PACKAGE_BUGREPORT, __func__, operator);
+        }
+
+      /* Allocating an struct for the output data. */
+      pp.out=gal_data_alloc(NULL, otype, outndim, pp.osize, NULL, 0,
+                            pp.input->minmapsize, pp.input->quietmmap, NULL,
+                            NULL, NULL);
+
+      /* Spin-off the threads and do the processing on each thread. */
+      gal_threads_spin_off(pool_type_on_thread, &pp, pp.out->size, numthreads,
+                           pp.input->minmapsize, pp.input->quietmmap);
+    }
+
+  /* Correct the WCS (if it has one). */
+  if(input->wcs)
+    {
+      /* We currently assume that a 'cdelt' exists (due to a lack of
+         time)! */
+      if(pp.input->wcs->cdelt==NULL)
+        error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+              "fix the problem. The input WCS has no 'cdelt' component",
+              __func__, PACKAGE_BUGREPORT);
+
+      /* Copy the input WCS to the output and correct the values. */
+      pp.out->wcs=gal_wcs_copy(pp.input->wcs);
+      pp.out->wcs->crpix[0]/=psize;
+      pp.out->wcs->crpix[1]/=psize;
+      pp.out->wcs->cdelt[0]*=psize;
+      pp.out->wcs->cdelt[1]*=psize;
+    }
+
+  /* Clean up and return. */
+  free(pp.osize);
+  return pp.out;
+}
+
+
+
+
+
+gal_data_t *
+gal_pool_max(gal_data_t *input, size_t psize, size_t numthreads)
+{
+  return pool_generic(input, psize, POOL_MAX, numthreads);
+}
+
+
+
+
+
+gal_data_t *
+gal_pool_min(gal_data_t *input, size_t psize, size_t numthreads)
+{
+  return pool_generic(input, psize, POOL_MIN, numthreads);
+}
+
+
+
+
+
+gal_data_t *
+gal_pool_sum(gal_data_t *input, size_t psize, size_t numthreads)
+{
+  return pool_generic(input, psize, POOL_SUM, numthreads);
+}
+
+
+
+
+
+gal_data_t *
+gal_pool_mean(gal_data_t *input, size_t psize, size_t numthreads)
+{
+  return pool_generic(input, psize, POOL_MEAN, numthreads);
+}
+
+
+
+
+
+gal_data_t *
+gal_pool_median(gal_data_t *input, size_t psize, size_t numthreads)
+{
+  return pool_generic(input, psize, POOL_MEDIAN, numthreads);
+}
diff --git a/lib/statistics.c b/lib/statistics.c
index 4990c70f..d255f390 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -368,7 +368,7 @@ gal_data_t *
 gal_statistics_median(gal_data_t *input, int inplace)
 {
   size_t dsize=1;
-  gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);;
+  gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
   gal_data_t *out=gal_data_alloc(NULL, nbs->type, 1, &dsize, NULL, 1, -1,
                                  1, NULL, NULL, NULL);