/* $Id: cpl_image_stats.c,v 1.76 2007/11/22 07:45:13 llundin Exp $ * * This file is part of the ESO Common Pipeline Library * Copyright (C) 2001-2004 European Southern Observatory * * This program 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 2 of the License, or * (at your option) any later version. * * This program 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 this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* * $Author: llundin $ * $Date: 2007/11/22 07:45:13 $ * $Revision: 1.76 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include #include #include #include #include "cpl_memory.h" #include "cpl_stats.h" #include "cpl_image_bpm.h" #include "cpl_image_stats.h" #include "cpl_mask.h" #include "cpl_tools.h" #include "cpl_image_defs.h" /*----------------------------------------------------------------------------- Defines -----------------------------------------------------------------------------*/ #define CPL_IMAGE_STATS_MEDIAN_SUBW 1 #define CPL_IMAGE_STATS_MEDIAN_STAT 2 /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_image_stats Statistics on images * * This module provides functions to compute various statistics in images. * * The bad pixel map stored in the image is taken into account for the * statistics computations. These functions can compute the mean, median, * maximum, minimum, flux, etc... of the good pixels of an image. * * @par Synopsis: * @code * #include "cpl_image_stats.h" * @endcode */ /*----------------------------------------------------------------------------*/ /**@{*/ /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /** @brief computes minimum pixel value over an image sub-window. @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the minimum value The specified bounds are included in the specified region. In case of error, the #_cpl_error_code_ code is set, and the returned double is undefined. Images can be CPL_TYPE_FLOAT, CPL_TYPE_INT or CPL_TYPE_DOUBLE. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ double cpl_image_get_min_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; cpl_ensure(image, CPL_ERROR_NULL_INPUT, 0.0); stats = cpl_stats_new_from_image_window(image, CPL_STATS_MIN, llx, lly, urx, ury) ; if (stats == NULL) return 0; res = cpl_stats_get_min(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief computes minimum pixel value over an image. @param image input image. @return the minimum value In case of error, the #_cpl_error_code_ code is set, and the returned double is undefined. Images can be CPL_TYPE_FLOAT, CPL_TYPE_INT or CPL_TYPE_DOUBLE. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ double cpl_image_get_min(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, 0.0); return cpl_image_get_min_window(image, 1, 1, image->nx, image->ny) ; } /*----------------------------------------------------------------------------*/ /** @brief computes maximum pixel value over an image sub-window. @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the maximum value @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_max_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_MAX, llx, lly, urx, ury) ; if (stats == NULL) return 0; res = cpl_stats_get_max(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief computes maximum pixel value over an image. @param image input image. @return the maximum value @see cpl_image_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_max(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, 0.0); return cpl_image_get_max_window(image, 1, 1, image->nx, image->ny) ; } /*----------------------------------------------------------------------------*/ /** @brief computes mean pixel value over an image sub-window. @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the mean value @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_mean_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_MEAN, llx, lly, urx, ury); if (stats == NULL) return 0; res = cpl_stats_get_mean(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief computes mean pixel value over an image. @param image input image. @return the mean value @see cpl_image_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_mean(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, 0.0); return cpl_image_get_mean_window(image, 1, 1, image->nx, image->ny) ; } #define CPL_OPERATION CPL_IMAGE_STATS_MEDIAN_SUBW /*----------------------------------------------------------------------------*/ /** @brief computes median pixel value over an image sub-window. @param image Input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return The median value, or undefined on error The specified bounds are included in the specified region. In case of error, the #_cpl_error_code_ code is set, and the returned value is undefined. Images can be CPL_TYPE_FLOAT, CPL_TYPE_INT or CPL_TYPE_DOUBLE. For a finite population or sample, the median is the middle value of an odd number of values (arranged in ascending order) or any value between the two middle values of an even number of values. In case of an even number of values in the input array, the median is chosen to be the lower of the two central values - this ensures that the computed median is always among the values of the input arrays. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if the window is outside the image or if there are only bad pixels in the window - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ double cpl_image_get_median_window(const cpl_image * image, int llx, int lly, int urx, int ury) { double median = 0.0; /* Avoid uninit warning */ int ngood = 0; /* Test inputs */ cpl_ensure(image != NULL, CPL_ERROR_NULL_INPUT, 0.0); cpl_ensure(llx <= urx, CPL_ERROR_ILLEGAL_INPUT, 0.0); cpl_ensure(llx > 0, CPL_ERROR_ILLEGAL_INPUT, 0.0); cpl_ensure(urx <= image->nx, CPL_ERROR_ILLEGAL_INPUT, 0.0); cpl_ensure(lly <= ury, CPL_ERROR_ILLEGAL_INPUT, 0.0); cpl_ensure(lly > 0, CPL_ERROR_ILLEGAL_INPUT, 0.0); cpl_ensure(ury <= image->ny, CPL_ERROR_ILLEGAL_INPUT, 0.0) ; /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_stats_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_stats_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_stats_body.h" #undef CPL_CLASS default: cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, 0.0); } /* Verify that there are good pixels */ cpl_ensure(ngood > 0, CPL_ERROR_ILLEGAL_INPUT, 0.0); return median; } #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief computes median pixel value over an image. @param image Input image. @return the median value @see cpl_image_get_median_window() In case of error, the #_cpl_error_code_ code is set, and the returned double is undefined. Images can be CPL_TYPE_FLOAT, CPL_TYPE_INT or CPL_TYPE_DOUBLE. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ double cpl_image_get_median(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, 0.0); return cpl_image_get_median_window(image, 1, 1, image->nx, image->ny) ; } /*----------------------------------------------------------------------------*/ /** @brief computes pixel standard deviation over an image sub-window. @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the standard deviation value @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_stdev_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_STDEV, llx, lly, urx, ury); if (stats == NULL) return 0; res = cpl_stats_get_stdev(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief computes pixel standard deviation over an image. @param image input image. @return the standard deviation value @see cpl_image_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_stdev(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, -1.0); return cpl_image_get_stdev_window(image, 1, 1, image->nx, image->ny) ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the sum of pixel values over an image sub-window @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the flux value @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_flux_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_FLUX, llx, lly, urx, ury); if (stats == NULL) return 0; res = cpl_stats_get_flux(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the sum of pixel values over an image. @param image input image. @return the flux value @see cpl_image_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_flux(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, 0.0); return cpl_image_get_flux_window(image, 1, 1, image->nx, image->ny) ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the sum of absolute values over an image sub-window @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the absolute flux (sum of |pixels|) value @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_absflux_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_ABSFLUX, llx, lly, urx, ury); if (stats == NULL) return -1; res = cpl_stats_get_absflux(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the sum of squared values over an image sub-window @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the square flux @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_sqflux_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_SQFLUX, llx, lly, urx, ury); if (stats == NULL) return -1; res = cpl_stats_get_sqflux(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the sum of absolute values over an image @param image input image. @return the absolute flux (sum of |pixels|) value @see cpl_image_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_absflux(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, -1.0); return cpl_image_get_absflux_window(image, 1, 1, image->nx, image->ny) ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the sum of squared values over an image @param image input image. @return the sqaure flux @see cpl_image_get_min() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_sqflux(const cpl_image * image) { cpl_ensure(image, CPL_ERROR_NULL_INPUT, -1.0); return cpl_image_get_sqflux_window(image, 1, 1, image->nx, image->ny) ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the x centroid value over an image sub-window @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the x centroid value @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_centroid_x_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_CENTROID, llx, lly, urx, ury); if (stats == NULL) return -1; res = cpl_stats_get_centroid_x(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief Computes the y centroid value over an image sub-window @param image input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @return the y centroid value @see cpl_image_get_min_window() */ /*----------------------------------------------------------------------------*/ double cpl_image_get_centroid_y_window( const cpl_image * image, int llx, int lly, int urx, int ury) { cpl_stats * stats ; double res ; stats = cpl_stats_new_from_image_window(image, CPL_STATS_CENTROID, llx, lly, urx, ury); if (stats == NULL) return -1; res = cpl_stats_get_centroid_y(stats) ; cpl_stats_delete(stats) ; return res ; } /*----------------------------------------------------------------------------*/ /** @brief Computes minimum pixel value and position over an image sub window. @param image Input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @param px ptr to the x coordinate of the minimum position @param py ptr to the y coordinate of the minimum position @return the #_cpl_error_code_ or CPL_ERROR_NONE The specified bounds are included in the specified region. Images can be CPL_TYPE_FLOAT, CPL_TYPE_INT or CPL_TYPE_DOUBLE. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_get_minpos_window( const cpl_image * image, int llx, int lly, int urx, int ury, int * px, int * py) { cpl_stats * stats ; cpl_ensure_code(image, CPL_ERROR_NULL_INPUT); cpl_ensure_code(px, CPL_ERROR_NULL_INPUT); cpl_ensure_code(py, CPL_ERROR_NULL_INPUT); stats = cpl_stats_new_from_image_window(image, CPL_STATS_MINPOS, llx, lly, urx, ury); if (stats == NULL) return cpl_error_get_code(); *px = cpl_stats_get_min_x(stats) ; *py = cpl_stats_get_min_y(stats) ; cpl_stats_delete(stats) ; return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Computes minimum pixel value and position over an image. @param image Input image. @param px ptr to the x coordinate of the minimum position @param py ptr to the y coordinate of the minimum position @return the #_cpl_error_code_ or CPL_ERROR_NONE Images can be CPL_TYPE_FLOAT, CPL_TYPE_INT or CPL_TYPE_DOUBLE. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_get_minpos( const cpl_image * image, int * px, int * py) { cpl_ensure_code(image, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(px, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(py, CPL_ERROR_NULL_INPUT) ; return cpl_image_get_minpos_window(image, 1, 1, image->nx, image->ny, px, py) ; } /*----------------------------------------------------------------------------*/ /** @brief Computes maximum pixel value and position over an image sub window. @param image Input image. @param llx Lower left x position (FITS convention) @param lly Lower left y position (FITS convention) @param urx Upper right x position (FITS convention) @param ury Upper right y position (FITS convention) @param px ptr to the x coordinate of the maximum position @param py ptr to the y coordinate of the maximum position @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_image_get_minpos_window() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_get_maxpos_window( const cpl_image * image, int llx, int lly, int urx, int ury, int * px, int * py) { cpl_stats * stats ; cpl_ensure_code(image, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(px, CPL_ERROR_NULL_INPUT); cpl_ensure_code(py, CPL_ERROR_NULL_INPUT); stats = cpl_stats_new_from_image_window(image, CPL_STATS_MAXPOS, llx, lly, urx, ury); if (stats == NULL) return cpl_error_get_code(); *px = cpl_stats_get_max_x(stats) ; *py = cpl_stats_get_max_y(stats) ; cpl_stats_delete(stats) ; return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Computes maximum pixel value and position over an image. @param image Input image. @param px ptr to the x coordinate of the maximum position @param py ptr to the y coordinate of the maximum position @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_image_get_minpos() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_get_maxpos( const cpl_image * image, int * px, int * py) { cpl_ensure_code(image, CPL_ERROR_NULL_INPUT) ; return cpl_image_get_maxpos_window(image, 1, 1, image->nx, image->ny, px, py) ; } #define CPL_OPERATION CPL_IMAGE_STATS_MEDIAN_STAT /*----------------------------------------------------------------------------*/ /** @brief Computes first and second order image statistics using median. @param image Input image. @param sigma Output computed sigma value. @return the median @see cpl_image_get_median_window() This function takes an image in input. It tries to estimate the average and standard deviation of the image by approximating them by resp. the median and the average absolute distance to the median. The median is the returned double. The average absolute distance to the median is written into sigma. Images can be CPL_TYPE_FLOAT, CPL_TYPE_INT or CPL_TYPE_DOUBLE. In case of error, the #_cpl_error_code_ code is set, and the returned double is undefined. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_ILLEGAL_INPUT if the image bad pixels map cannot be retrieved - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported - CPL_ERROR_DATA_NOT_FOUND if all the pixels are bad in image */ /*----------------------------------------------------------------------------*/ double cpl_image_get_median_dev( const cpl_image * image, double * sigma) { double median_val ; int npix ; int i ; /* The number of bad pixels inside the subwindow */ int nbadpix ; /* A map of the the bad pixels in the input */ cpl_mask * badmask ; cpl_binary * badmap ; /* Test entries */ cpl_ensure(image, CPL_ERROR_NULL_INPUT, 0.0); cpl_ensure(sigma, CPL_ERROR_NULL_INPUT, 0.0); if (image->bpm == NULL) badmask = cpl_mask_new(image->nx, image->ny) ; else badmask = cpl_mask_duplicate(image->bpm) ; badmap = cpl_mask_get_data(badmask) ; /* Get the first good pixel and the number of bad pixels */ nbadpix = cpl_mask_count(badmask) ; npix = image->nx * image->ny ; median_val = (double)cpl_image_get_median(image) ; *sigma=0.0; /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_stats_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_stats_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_stats_body.h" #undef CPL_CLASS default: cpl_mask_delete(badmask) ; cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, 0.0); } cpl_mask_delete(badmask) ; npix -= nbadpix; *sigma/= (double) npix ; return median_val; } #undef CPL_OPERATION /**@}*/