[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master a23d86b 1/2: Arithmetic's filtering operators
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master a23d86b 1/2: Arithmetic's filtering operators also work on blanks |
|
Date: |
Fri, 20 Apr 2018 12:42:03 -0400 (EDT) |
branch: master
commit a23d86b40d7e4bf0f6ac60dcb6b24e53afbe9404
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Arithmetic's filtering operators also work on blanks
Until now, Arithmetic's filtering operations ignored blank valued
pixels. However, in some cases, the user might want to use the filtering to
remove/interpolate over the blank pixels. So the condition to ignore blank
pixels is now removed. A paragraph was added to the `filter-mean' operator,
describing the effect and how to correct for it if it is important to
preserve blank elements.
---
NEWS | 5 ++
bin/arithmetic/arithmetic.c | 181 +++++++++++++++++++++-----------------------
doc/gnuastro.texi | 9 ++-
3 files changed, 98 insertions(+), 97 deletions(-)
diff --git a/NEWS b/NEWS
index beffdd6..45356df 100644
--- a/NEWS
+++ b/NEWS
@@ -107,6 +107,11 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
** Changed features
+ Arithmetic:
+ - filter-mean: a blank value in the input can be non-blank in the
+ output when non-blank elements present in filter.
+ - filter-median: similar to filter-mean.
+
Fits:
--history: can be called/written multiple times in one run.
--comment: can be called/written multiple times in one run.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 87d3a4d..63487bc 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -179,121 +179,110 @@ arithmetic_filter(void *in_prm)
/* For easy reading, put the index in `ind'. */
ind=tprm->indexs[i];
- /* If we are on a blank element, then just set the output to blank
- also. */
- if( afp->hasblank
- && gal_blank_is(gal_data_ptr_increment(input->array, ind,
- input->type), input->type) )
- gal_blank_write(gal_data_ptr_increment(afp->out->array, ind,
- afp->out->type),
- afp->out->type);
- else
- {
- /* Get the coordinate of the pixel. */
- gal_dimension_index_to_coord(ind, ndim, dsize, coord);
+ /* Get the coordinate of the pixel. */
+ gal_dimension_index_to_coord(ind, ndim, dsize, coord);
- /* See which dimensions need trimming. */
- tile->size=1;
- for(j=0;j<ndim;++j)
+ /* See which dimensions need trimming. */
+ tile->size=1;
+ for(j=0;j<ndim;++j)
+ {
+ /* Estimate the coordinate of the filter's starting point. Note
+ that we are dealing with size_t (unsigned int) type here, so
+ there are no negatives. A negative result will produce an
+ extremely large number, so instead of checking for negative,
+ we can just see if the result of a subtraction is less than
+ the width of the input. */
+ if( (coord[j] - hnfsize[j] > dsize[j])
+ || (coord[j] + hpfsize[j] >= dsize[j]) )
{
- /* Estimate the coordinate of the filter's starting
- point. Note that we are dealing with size_t (unsigned int)
- type here, so there are no negatives. A negative result
- will produce an extremely large number, so instead of
- checking for negative, we can just see if the result of a
- subtraction is less than the width of the input. */
- if( (coord[j] - hnfsize[j] > dsize[j])
- || (coord[j] + hpfsize[j] >= dsize[j]) )
- {
- start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
- ? 0 : coord[j] - hnfsize[j] );
- end[j] = ( (coord[j] + hpfsize[j] >= dsize[j])
- ? dsize[j]
- : coord[j] + hpfsize[j] + 1);
- tsize[j] = end[j] - start[j];
- }
- else /* NOT on the edge (given requested filter width). */
- {
- tsize[j] = fsize[j];
- start[j] = coord[j] - hnfsize[j];
- }
- tile->size *= tsize[j];
+ start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
+ ? 0 : coord[j] - hnfsize[j] );
+ end[j] = ( (coord[j] + hpfsize[j] >= dsize[j])
+ ? dsize[j]
+ : coord[j] + hpfsize[j] + 1);
+ tsize[j] = end[j] - start[j];
}
-
- /* For a test.
- printf("coord: %zu, %zu\n", coord[1]+1, coord[0]+1);
- printf("\tstart: %zu, %zu\n", start[1]+1, start[0]+1);
- printf("\ttsize: %zu, %zu\n", tsize[1], tsize[0]);
- */
-
- /* Set the tile's starting pointer. */
- index=gal_dimension_coord_to_index(ndim, dsize, start);
- tile->array=gal_data_ptr_increment(input->array, index,
- input->type);
-
- /* Do the necessary calculation. */
- switch(afp->operator)
+ else /* NOT on the edge (given requested filter width). */
{
- case ARITHMETIC_OP_FILTER_MEDIAN:
- result=gal_statistics_median(tile, 0);
- break;
-
+ tsize[j] = fsize[j];
+ start[j] = coord[j] - hnfsize[j];
+ }
+ tile->size *= tsize[j];
+ }
- case ARITHMETIC_OP_FILTER_MEAN:
- result=gal_statistics_mean(tile);
- break;
+ /* For a test.
+ printf("coord: %zu, %zu\n", coord[1]+1, coord[0]+1);
+ printf("\tstart: %zu, %zu\n", start[1]+1, start[0]+1);
+ printf("\ttsize: %zu, %zu\n", tsize[1], tsize[0]);
+ */
+ /* Set the tile's starting pointer. */
+ index=gal_dimension_coord_to_index(ndim, dsize, start);
+ tile->array=gal_data_ptr_increment(input->array, index,
+ input->type);
- case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
- case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
- /* Find the sigma-clipped results. */
- sigclip=gal_statistics_sigma_clip(tile, afp->sclip_multip,
- afp->sclip_param, 0, 1);
+ /* Do the necessary calculation. */
+ switch(afp->operator)
+ {
+ case ARITHMETIC_OP_FILTER_MEDIAN:
+ result=gal_statistics_median(tile, 0);
+ break;
- /* Set the required index. */
- switch(afp->operator)
- {
- case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN: sind = 2; break;
- case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN: sind = 1; break;
- default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
- "%s to fix the problem. The `afp->operator' value "
- "%d is not recognized as sigma-clipped median or "
- "mean", __func__, PACKAGE_BUGREPORT, afp->operator);
- }
- /* Allocate the output and write the value into it. */
- result=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL,
- 0, -1, NULL, NULL, NULL);
- ((float *)(result->array))[0] =
- ((float *)(sigclip->array))[sind];
+ case ARITHMETIC_OP_FILTER_MEAN:
+ result=gal_statistics_mean(tile);
+ break;
- /* Clean up. */
- gal_data_free(sigclip);
- break;
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
+ /* Find the sigma-clipped results. */
+ sigclip=gal_statistics_sigma_clip(tile, afp->sclip_multip,
+ afp->sclip_param, 0, 1);
+ /* Set the required index. */
+ switch(afp->operator)
+ {
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN: sind = 2; break;
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN: sind = 1; break;
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
- "to fix the problem. `afp->operator' code %d is not "
- "recognized", PACKAGE_BUGREPORT, __func__,
- afp->operator);
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
+ "%s to fix the problem. The `afp->operator' value "
+ "%d is not recognized as sigma-clipped median or "
+ "mean", __func__, PACKAGE_BUGREPORT, afp->operator);
}
- /* Make sure the output array type and result's type are the
- same. */
- if(result->type!=afp->out->type)
- result=gal_data_copy_to_new_type_free(result, afp->out->type);
+ /* Allocate the output and write the value into it. */
+ result=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL,
+ 0, -1, NULL, NULL, NULL);
+ ((float *)(result->array))[0] =
+ ((float *)(sigclip->array))[sind];
+ /* Clean up. */
+ gal_data_free(sigclip);
+ break;
- /* Copy the result into the output array. */
- memcpy(gal_data_ptr_increment(afp->out->array, ind,
- afp->out->type),
- result->array, gal_type_sizeof(afp->out->type));
- /* Clean up for this pixel. */
- gal_data_free(result);
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+ "to fix the problem. `afp->operator' code %d is not "
+ "recognized", PACKAGE_BUGREPORT, __func__,
+ afp->operator);
}
+
+ /* Make sure the output array type and result's type are the
+ same. */
+ if(result->type!=afp->out->type)
+ result=gal_data_copy_to_new_type_free(result, afp->out->type);
+
+
+ /* Copy the result into the output array. */
+ memcpy(gal_data_ptr_increment(afp->out->array, ind,
+ afp->out->type),
+ result->array, gal_type_sizeof(afp->out->type));
+
+ /* Clean up for this pixel. */
+ gal_data_free(result);
}
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8476790..4f67c8a 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10118,6 +10118,13 @@ The final effect of mean filtering is to smooth the
input image, it is
essentially a convolution with a kernel that has identical values for all
its pixels (is flat), see @ref{Convolution process}.
+Note that blank pixels will also be affected by this operator: if there are
+any non-blank elements in the box surrounding a blank pixel, in the
+flitered image, it will have the mean of the non-blank elements, therefore
+it won't be blank any more. If blank elements are important for your
+analysis, you can use the @code{isblank} with the @code{where} operator to
+set them back to blank after filtering.
+
@item filter-median
Apply @url{https://en.wikipedia.org/wiki/Median_filter, median filtering}
on the input dataset. This is very similar to @command{filter-mean}, except
@@ -25375,7 +25382,7 @@ nature. The produced values are also always within the
range of the known
values and strong outliers do not get created. We will hopefully implement
other methods too (wrappers around the GNU Scientific Library's very
complete set of functions), but currently the developers are too busy. So
-if you do have the chance to help your contribution would be very welcome
+if you do have the chance to help, your contribution would be very welcome
and appreciated.
@deftypefun {gal_data_t *} gal_interpolate_close_neighbors (gal_data_t
@code{*input}, struct gal_tile_two_layer_params @code{*tl}, size_t
@code{numneighbors}, size_t @code{numthreads}, int @code{onlyblank}, int
@code{aslinkedlist})