/* $Id: cpl_apertures.c,v 1.6 2007/08/09 16:05:15 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/08/09 16:05:15 $ * $Revision: 1.6 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include /* Needed by memcpy() */ #include #include #include #include #include #include #include #include #include "cpl_apertures.h" /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_apertures High level functions to handle apertures * * The aperture object contains a list of zones in an image. It is * typically used to contain the results of an objects detection, or if * one wants to work on a very specific zone in an image. * * This module provides functions to handle @em cpl_apertures. */ /*----------------------------------------------------------------------------*/ /**@{*/ /*----------------------------------------------------------------------------- Type definition -----------------------------------------------------------------------------*/ struct _cpl_apertures_ { /* Number of apertures */ int naperts ; /* All positions are in FITS conventions : (1,1) is the lower left) */ /* Aperture center: x=sum(x*1)/sum(1) */ double * x ; double * y ; /* Aperture weighted center : xcentroid=sum(x*f(x)) / sum(f(x)) */ double * xcentroid ; double * ycentroid ; int * npix ; int * left_x ; int * left_y ; int * right_x ; int * right_y ; int * top_x ; int * top_y ; int * bottom_x ; int * bottom_y ; double * max_val ; double * min_val ; double * mean ; double * median ; double * stdev ; double * flux ; } ; /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /** @brief Constructor for cpl_apertures @param naperts Number of apertures in the structure @return 1 newly allocated cpl_apertures or NULL in error case The returned object must be deleted using cpl_apertures_delete(). Possible #_cpl_error_code_ set in this function: - CPL_ERROR_ILLEGAL_INPUT if naperts is not strictly positive */ /*----------------------------------------------------------------------------*/ cpl_apertures * cpl_apertures_new(int naperts) { cpl_apertures * aperts ; /* Check entries */ cpl_ensure(naperts>0, CPL_ERROR_ILLEGAL_INPUT, NULL) ; /* Create a cpl_apertures */ aperts = cpl_calloc(1, sizeof(cpl_apertures)) ; aperts->naperts = naperts ; /* Allocate data holders */ aperts->x = cpl_calloc(naperts, sizeof(double)) ; aperts->y = cpl_calloc(naperts, sizeof(double)) ; aperts->xcentroid = cpl_calloc(naperts, sizeof(double)) ; aperts->ycentroid = cpl_calloc(naperts, sizeof(double)) ; aperts->npix = cpl_calloc(naperts, sizeof(int)) ; aperts->left_x = cpl_calloc(naperts, sizeof(int)) ; aperts->left_y = cpl_calloc(naperts, sizeof(int)) ; aperts->right_x = cpl_calloc(naperts, sizeof(int)) ; aperts->right_y = cpl_calloc(naperts, sizeof(int)) ; aperts->top_x = cpl_calloc(naperts, sizeof(int)) ; aperts->top_y = cpl_calloc(naperts, sizeof(int)) ; aperts->bottom_x = cpl_calloc(naperts, sizeof(int)) ; aperts->bottom_y = cpl_calloc(naperts, sizeof(int)) ; aperts->max_val = cpl_calloc(naperts, sizeof(double)) ; aperts->min_val = cpl_calloc(naperts, sizeof(double)) ; aperts->mean = cpl_calloc(naperts, sizeof(double)) ; aperts->median = cpl_calloc(naperts, sizeof(double)) ; aperts->stdev = cpl_calloc(naperts, sizeof(double)) ; aperts->flux = cpl_calloc(naperts, sizeof(double)) ; /* Return */ return aperts ; } /*----------------------------------------------------------------------------*/ /** @brief Destructor for cpl_apertures @param apert Object to delete. @return void This function deallocates all possibly allocated arrays inside the given object, then deallocates the main pointer. */ /*----------------------------------------------------------------------------*/ void cpl_apertures_delete(cpl_apertures * apert) { if (apert==NULL) return ; if (apert->x!=NULL) cpl_free(apert->x); if (apert->y!=NULL) cpl_free(apert->y); if (apert->xcentroid!=NULL) cpl_free(apert->xcentroid); if (apert->ycentroid!=NULL) cpl_free(apert->ycentroid); if (apert->npix!=NULL) cpl_free(apert->npix); if (apert->left_x!=NULL) cpl_free(apert->left_x); if (apert->left_y!=NULL) cpl_free(apert->left_y); if (apert->right_x!=NULL) cpl_free(apert->right_x); if (apert->right_y!=NULL) cpl_free(apert->right_y); if (apert->top_x!=NULL) cpl_free(apert->top_x); if (apert->top_y!=NULL) cpl_free(apert->top_y); if (apert->bottom_x!=NULL) cpl_free(apert->bottom_x); if (apert->bottom_y!=NULL) cpl_free(apert->bottom_y); if (apert->max_val!=NULL) cpl_free(apert->max_val); if (apert->min_val!=NULL) cpl_free(apert->min_val); if (apert->mean!=NULL) cpl_free(apert->mean); if (apert->median!=NULL) cpl_free(apert->median); if (apert->stdev!=NULL) cpl_free(apert->stdev); if (apert->flux!=NULL) cpl_free(apert->flux); cpl_free(apert); return ; } /*----------------------------------------------------------------------------*/ /** @brief Dump a cpl_apertures to an opened file pointer. @param aperts cpl_apertures to dump @param fp Opened file pointer, ready to receive data @return void This function dumps all informations contained into a cpl_apertures to the passed (opened) file pointer. It is Ok to pass stdout or stderr. If the object is unallocated or contains nothing, this function does nothing. */ /*----------------------------------------------------------------------------*/ void cpl_apertures_dump( const cpl_apertures * aperts, FILE * fp) { int i ; if (aperts==NULL || fp==NULL) return ; if (aperts->naperts<1) return ; fprintf(fp, "# X Y"); fprintf(fp, " XCENTROID YCENTROID"); if (aperts->npix != NULL) fprintf(fp, " pix") ; if (aperts->max_val != NULL) fprintf(fp, " max") ; if (aperts->min_val != NULL) fprintf(fp, " min") ; if (aperts->mean != NULL) fprintf(fp, " mean") ; if (aperts->median != NULL) fprintf(fp, " med") ; if (aperts->stdev != NULL) fprintf(fp, " dev") ; if (aperts->flux != NULL) fprintf(fp, " flux") ; fprintf(fp, "\n"); for (i=0 ; inaperts ; i++) { fprintf(fp, "% 3d %6.1f %6.1f", i+1, aperts->x[i], aperts->y[i]); fprintf(fp, " %6.1f %6.1f", aperts->xcentroid[i], aperts->ycentroid[i]); if (aperts->npix != NULL) fprintf(fp, " % 6d", aperts->npix[i]) ; if (aperts->max_val != NULL) fprintf(fp, " %6.2f", aperts->max_val[i]); if (aperts->min_val != NULL) fprintf(fp, " %6.2f", aperts->min_val[i]); if (aperts->mean != NULL) fprintf(fp, " %6.2f", aperts->mean[i]); if (aperts->median != NULL) fprintf(fp, " %6.2f", aperts->median[i]); if (aperts->stdev != NULL) fprintf(fp, " %6.2f", aperts->stdev[i]); if (aperts->flux != NULL) fprintf(fp, " %8.2f", aperts->flux[i]); fprintf(fp, "\n"); } return ; } /*----------------------------------------------------------------------------*/ /** @brief Compute statistics on selected apertures @param in Reference image. @param lab labelized image (type CPL_TYPE_INT) @return An apertures statistics holder or NULL in error case The returned object must be deleted using cpl_apertures_delete(). The labelized image must contain at least one pixel for each successive values from 1 to its maximum value. For the centroiding computation of an aperture, if some pixels have values lower or equal to 0, all the values of the aperture are locally shifted such as the minimum value of the aperture has a value of epsilon. The centroid is then computed on these positive values. In principle, centroid should always be computed on positive values, this is done to avoid raising an error in case the caller of the function wants to use it on negative values images without caring about the centroid results. In such cases, the centroid result would be meaningful, but slightly depend on the hardcoded value chosen for epsilon (1e-8). 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 lab is not of type CPL_TYPE_INT or if lab and in have different sizes or if lab maximum value is lower or equal to 0 or if one of the lab values is missing. */ /*----------------------------------------------------------------------------*/ cpl_apertures * cpl_apertures_new_from_image( const cpl_image * in, const cpl_image * lab) { cpl_apertures * aperts ; double min ; cpl_image * dimage ; const cpl_mask * bpm ; cpl_image * lab_local ; int naperts ; double * pd ; int * pi ; double * sqsum ; double pix ; int npix ; int * plab ; double * med_array ; double weight; int count ; int nx, ny ; double ep ; int i, j, k; /* Initialise */ nx = cpl_image_get_size_x(in) ; ny = cpl_image_get_size_y(in) ; ep = 1e-10 ; /* Test entries */ cpl_ensure(in && lab, CPL_ERROR_NULL_INPUT, NULL) ; cpl_ensure(cpl_image_get_type(lab)==CPL_TYPE_INT, CPL_ERROR_ILLEGAL_INPUT, NULL) ; cpl_ensure(nx==cpl_image_get_size_x(lab) && ny==cpl_image_get_size_y(lab), CPL_ERROR_ILLEGAL_INPUT, NULL) ; /* Set the bad pixels of the input image to 0 in the lab image */ lab_local = cpl_image_duplicate(lab) ; bpm = cpl_image_get_bpm_const(in) ; if (bpm != NULL) { cpl_image_reject_from_mask(lab_local, bpm) ; cpl_image_fill_rejected(lab_local, 0.0) ; } /* Initialise */ pi = cpl_image_get_data_int(lab_local) ; /* Compute the minimum of the image */ min = cpl_image_get_min(in) ; /* Get number of apertures from the labelised image */ if ((naperts = (int)cpl_image_get_max(lab_local)) <= 0) { cpl_image_delete(lab_local) ; cpl_ensure(0, CPL_ERROR_ILLEGAL_INPUT, NULL) ; } /* Convert input image in a double image */ dimage = cpl_image_cast(in, CPL_TYPE_DOUBLE) ; pd = cpl_image_get_data_double(dimage) ; /* Create a cpl_apertures */ aperts = cpl_apertures_new(naperts) ; sqsum = cpl_calloc(naperts, sizeof(double)) ; /* Initialise min and max */ for (i=0 ; imax_val[i] = -DBL_MAX; aperts->min_val[i] = DBL_MAX; } /* Compute some stats */ for (j=0 ; jx[k] += (double)(i+1) ; aperts->y[k] += (double)(j+1) ; /* Increase number of pixels */ aperts->npix[k] ++ ; /* Store pixel sum and squared sum */ aperts->flux[k] += pix ; sqsum[k] += (pix * pix) ; /* Check max pos value */ if (pix > aperts->max_val[k]) aperts->max_val[k] = pix ; /* Check min pos value */ if (pix < aperts->min_val[k]) aperts->min_val[k] = pix ; /* Check object extremities */ if ((i+1 < aperts->left_x[k]) || (aperts->left_x[k] == 0)) aperts->left_x[k] = i+1 ; if ((i+1 < aperts->left_y[k]) || (aperts->left_y[k] == 0)) aperts->left_y[k] = j+1 ; if ((i+1 > aperts->right_x[k]) || (aperts->right_x[k] == 0)) aperts->right_x[k] = i+1 ; if ((i+1 > aperts->right_y[k]) || (aperts->right_y[k] == 0)) aperts->right_y[k] = j+1 ; if ((j+1 < aperts->bottom_x[k]) || (aperts->bottom_x[k] == 0)) aperts->bottom_x[k] = i+1 ; if ((j+1 < aperts->bottom_y[k]) || (aperts->bottom_y[k] == 0)) aperts->bottom_y[k] = j+1 ; if ((j+1 > aperts->top_x[k]) || (aperts->top_x[k] == 0)) aperts->top_x[k] = i+1 ; if ((j+1 > aperts->top_y[k]) || (aperts->top_y[k] == 0)) aperts->top_y[k] = j+1 ; } } /* Centroiding */ for (j=0 ; jmin_val[k] <= 0) { /* Accumulate weighted positions */ aperts->xcentroid[k]+=(double)(i+1)*(pix-aperts->min_val[k]+ep); aperts->ycentroid[k]+=(double)(j+1)*(pix-aperts->min_val[k]+ep); } else { /* Accumulate weighted positions */ aperts->xcentroid[k] += (double)(i+1) * pix ; aperts->ycentroid[k] += (double)(j+1) * pix ; } } } /* Compute average and std dev for each aperture, normalize centers */ for (k=0 ; knaperts ; k++) { npix = aperts->npix[k] ; if (npix<=0) { cpl_free(sqsum); cpl_apertures_delete(aperts) ; cpl_image_delete(dimage) ; cpl_image_delete(lab_local) ; cpl_ensure(0, CPL_ERROR_ILLEGAL_INPUT, NULL) ; } aperts->mean[k] = aperts->flux[k] / (double)npix ; if (npix>1) { aperts->stdev[k] = (sqsum[k] - ((aperts->flux[k]*aperts->flux[k])/(double)npix)) / ((double)npix-1.0); /* Can become negative due to rounding */ aperts->stdev[k] = aperts->stdev[k] >0 ? sqrt(aperts->stdev[k]) : 0; } else { aperts->stdev[k] = -1; } aperts->x[k] /= (double)npix ; aperts->y[k] /= (double)npix ; weight = aperts->flux[k]; if (aperts->min_val[k] <= 0) { weight -= (double)(npix * aperts->min_val[k]); weight += (double)(npix * ep) ; } if (aperts->xcentroid[k] < weight * aperts->left_x[k] || aperts->xcentroid[k] > weight * aperts->right_x[k] || weight == 0) { aperts->xcentroid[k] = 0; } else { aperts->xcentroid[k] /= weight; } if (aperts->ycentroid[k] < weight * aperts->bottom_y[k] || aperts->ycentroid[k] > weight * aperts->top_y[k] || weight == 0) { aperts->ycentroid[k] = 0; } else { aperts->ycentroid[k] /= weight; } } cpl_free(sqsum); /* Compute median for each aperture */ plab = cpl_image_get_data_int(lab_local) ; for (k=0 ; knpix[k] * sizeof(double)); count=0 ; for (j=0 ; jmedian[k] = cpl_tools_get_median_double(med_array, count); cpl_free(med_array) ; } /* Free and return */ cpl_image_delete(dimage) ; cpl_image_delete(lab_local) ; return aperts ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the number of apertures @param in a cpl_apertures object @return the number of apertures or -1 in error case Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_size(const cpl_apertures * in) { cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; return in->naperts ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the x position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the x position of the aperture 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_max_x( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->x[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the y position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the y position of the aperture @see cpl_apertures_get_max_x() */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_max_y( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->y[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the x centroid of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the x centroid position of the aperture 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_centroid_x( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->xcentroid[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the y centroid of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the y centroid position of the aperture 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_centroid_y( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->ycentroid[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the number of pixels of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the number of pixels of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_npix( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->npix[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the left x position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the left x position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_left( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->left_x[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the left y position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the left y position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_left_y( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->left_y[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the right x position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the right x position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_right( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->right_x[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the right y position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the right y position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_right_y( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->right_y[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the bottom x position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the bottom x position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_bottom_x( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->bottom_x[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the bottom y position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the bottom y position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_bottom( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->bottom_y[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the top x position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the top x position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_top_x( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->top_x[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the top y position of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the top y position of the aperture or -1 in error case 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ int cpl_apertures_get_top( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT,-1); return in->top_y[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the maximum value of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the maximum value of the aperture 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_max( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->max_val[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the minimum value of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the minimum value of the aperture 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_min( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->min_val[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the mean value of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the mean value of the aperture 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 ind is out of the bounds of in */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_mean( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->mean[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the median value of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the median value of the aperture @see cpl_apertures_get_mean() */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_median( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->median[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the std dev. value of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the standard deviation of the aperture @see cpl_apertures_get_mean() */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_stdev( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; /* A negative value is used as internal representation of an undefined stdev. An error is raied only when the stdev is asked for */ cpl_ensure(in->stdev[ind-1] >= 0, CPL_ERROR_DATA_NOT_FOUND, 0.00) ; return in->stdev[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Accessor to get the flux of an aperture @param in a cpl_apertures object @param ind the aperture index (1 for the first one) @return the flux of the aperture @see cpl_apertures_get_mean() */ /*----------------------------------------------------------------------------*/ double cpl_apertures_get_flux( const cpl_apertures * in, int ind) { /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, 0.00) ; cpl_ensure((ind>0) && (ind<=in->naperts), CPL_ERROR_ILLEGAL_INPUT, 0.00) ; return in->flux[ind-1] ; } /*----------------------------------------------------------------------------*/ /** @brief Sort by decreasing aperture size @param in Apertures to sort (MODIFIED) @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ 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_apertures_sort_by_npix(cpl_apertures * in) { cpl_apertures * sorted_apert ; void * swap ; int naperts ; int * sorted_ind ; int * computed ; int max_npix ; int max_ind = -1; /* Avoid (false) uninit warning */ int i, j ; /* Test entries */ cpl_ensure_code(in, CPL_ERROR_NULL_INPUT); /* Initialize */ naperts = in->naperts ; /* Compute the sorted index array - FIXME: Use better sort algorithm ? */ sorted_ind = cpl_malloc(naperts * sizeof(int)) ; computed = cpl_calloc(naperts, sizeof(int)) ; for (j=0 ; jnpix[i])) { max_npix = in->npix[i] ; max_ind = i ; } } computed[max_ind] = 1 ; sorted_ind[j] = max_ind ; } cpl_free(computed) ; /* Now sort the input apertures */ sorted_apert = cpl_apertures_new(naperts) ; for (i=0 ; inaperts ; i++) { sorted_apert->x[i] = in->x[sorted_ind[i]] ; sorted_apert->y[i] = in->y[sorted_ind[i]] ; sorted_apert->xcentroid[i] = in->xcentroid[sorted_ind[i]] ; sorted_apert->ycentroid[i] = in->ycentroid[sorted_ind[i]] ; sorted_apert->npix[i] = in->npix[sorted_ind[i]] ; sorted_apert->left_x[i] = in->left_x[sorted_ind[i]] ; sorted_apert->left_y[i] = in->left_y[sorted_ind[i]] ; sorted_apert->right_x[i] = in->right_x[sorted_ind[i]] ; sorted_apert->right_y[i] = in->right_y[sorted_ind[i]] ; sorted_apert->top_x[i] = in->top_x[sorted_ind[i]] ; sorted_apert->top_y[i] = in->top_y[sorted_ind[i]] ; sorted_apert->bottom_x[i] = in->bottom_x[sorted_ind[i]] ; sorted_apert->bottom_y[i] = in->bottom_y[sorted_ind[i]] ; sorted_apert->max_val[i] = in->max_val[sorted_ind[i]] ; sorted_apert->min_val[i] = in->min_val[sorted_ind[i]] ; sorted_apert->mean[i] = in->mean[sorted_ind[i]] ; sorted_apert->median[i] = in->median[sorted_ind[i]] ; sorted_apert->stdev[i] = in->stdev[sorted_ind[i]] ; sorted_apert->flux[i] = in->flux[sorted_ind[i]] ; } cpl_free(sorted_ind) ; /* Swap the sorted and the input apertures */ swap = cpl_malloc(sizeof(cpl_apertures)); memcpy(swap, sorted_apert, sizeof(cpl_apertures)); memcpy(sorted_apert, in, sizeof(cpl_apertures)); memcpy(in, swap, sizeof(cpl_apertures)); cpl_free(swap); cpl_apertures_delete(sorted_apert); return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief Sort by decreasing aperture peak value @param in Apertures to sort (MODIFIED) @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_apertures_sort_by_npix() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_apertures_sort_by_max(cpl_apertures * in) { cpl_apertures * sorted_apert ; void * swap ; int naperts ; int * sorted_ind ; int * computed ; double max_max = DBL_MAX; /* Avoid (false) uninit warning */ int max_ind ; int i, j ; /* Test entries */ cpl_ensure_code(in, CPL_ERROR_NULL_INPUT); /* Initialize */ naperts = in->naperts ; /* Compute the sorted index array - FIXME: Use better sort algorithm ? */ sorted_ind = cpl_malloc(naperts * sizeof(int)) ; computed = cpl_calloc(naperts, sizeof(int)) ; for (j=0 ; j < naperts ; j++) { max_ind = -1; for (i=0 ; i < naperts ; i++) { if ((computed[i]==0) && (max_ind < 0 || max_max < in->max_val[i])) { max_max = in->max_val[i] ; max_ind = i ; } } computed[max_ind] = 1 ; sorted_ind[j] = max_ind ; } cpl_free(computed) ; /* Now sort the input apertures */ sorted_apert = cpl_apertures_new(naperts) ; for (i=0 ; inaperts ; i++) { sorted_apert->x[i] = in->x[sorted_ind[i]] ; sorted_apert->y[i] = in->y[sorted_ind[i]] ; sorted_apert->xcentroid[i] = in->xcentroid[sorted_ind[i]] ; sorted_apert->ycentroid[i] = in->ycentroid[sorted_ind[i]] ; sorted_apert->npix[i] = in->npix[sorted_ind[i]] ; sorted_apert->left_x[i] = in->left_x[sorted_ind[i]] ; sorted_apert->left_y[i] = in->left_y[sorted_ind[i]] ; sorted_apert->right_x[i] = in->right_x[sorted_ind[i]] ; sorted_apert->right_y[i] = in->right_y[sorted_ind[i]] ; sorted_apert->top_x[i] = in->top_x[sorted_ind[i]] ; sorted_apert->top_y[i] = in->top_y[sorted_ind[i]] ; sorted_apert->bottom_x[i] = in->bottom_x[sorted_ind[i]] ; sorted_apert->bottom_y[i] = in->bottom_y[sorted_ind[i]] ; sorted_apert->max_val[i] = in->max_val[sorted_ind[i]] ; sorted_apert->min_val[i] = in->min_val[sorted_ind[i]] ; sorted_apert->mean[i] = in->mean[sorted_ind[i]] ; sorted_apert->median[i] = in->median[sorted_ind[i]] ; sorted_apert->stdev[i] = in->stdev[sorted_ind[i]] ; sorted_apert->flux[i] = in->flux[sorted_ind[i]] ; } cpl_free(sorted_ind) ; /* Swap the sorted and the input apertures */ swap = cpl_malloc(sizeof(cpl_apertures)); memcpy(swap, sorted_apert, sizeof(cpl_apertures)); memcpy(sorted_apert, in, sizeof(cpl_apertures)); memcpy(in, swap, sizeof(cpl_apertures)); cpl_free(swap); cpl_apertures_delete(sorted_apert); return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief Sort by decreasing aperture flux @param in Apertures to sort (MODIFIED) @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_apertures_sort_by_npix() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_apertures_sort_by_flux(cpl_apertures * in) { cpl_apertures * sorted_apert ; void * swap ; int naperts ; int * sorted_ind ; int * computed ; double max_flux = DBL_MAX; /* Avoid (false) uninit warning */ int max_ind ; int i, j ; /* Test entries */ cpl_ensure_code(in, CPL_ERROR_NULL_INPUT); /* Initialize */ naperts = in->naperts ; /* Compute the sorted index array - FIXME: Use better sort algorithm ? */ sorted_ind = cpl_malloc(naperts * sizeof(int)) ; computed = cpl_calloc(naperts, sizeof(int)) ; for (j=0 ; jflux[i])) { max_flux = in->flux[i] ; max_ind = i ; } } computed[max_ind] = 1 ; sorted_ind[j] = max_ind ; } cpl_free(computed) ; /* Now sort the input apertures */ sorted_apert = cpl_apertures_new(naperts) ; for (i=0 ; inaperts ; i++) { sorted_apert->x[i] = in->x[sorted_ind[i]] ; sorted_apert->y[i] = in->y[sorted_ind[i]] ; sorted_apert->xcentroid[i] = in->xcentroid[sorted_ind[i]] ; sorted_apert->ycentroid[i] = in->ycentroid[sorted_ind[i]] ; sorted_apert->npix[i] = in->npix[sorted_ind[i]] ; sorted_apert->left_x[i] = in->left_x[sorted_ind[i]] ; sorted_apert->left_y[i] = in->left_y[sorted_ind[i]] ; sorted_apert->right_x[i] = in->right_x[sorted_ind[i]] ; sorted_apert->right_y[i] = in->right_y[sorted_ind[i]] ; sorted_apert->top_x[i] = in->top_x[sorted_ind[i]] ; sorted_apert->top_y[i] = in->top_y[sorted_ind[i]] ; sorted_apert->bottom_x[i] = in->bottom_x[sorted_ind[i]] ; sorted_apert->bottom_y[i] = in->bottom_y[sorted_ind[i]] ; sorted_apert->max_val[i] = in->max_val[sorted_ind[i]] ; sorted_apert->min_val[i] = in->min_val[sorted_ind[i]] ; sorted_apert->mean[i] = in->mean[sorted_ind[i]] ; sorted_apert->median[i] = in->median[sorted_ind[i]] ; sorted_apert->stdev[i] = in->stdev[sorted_ind[i]] ; sorted_apert->flux[i] = in->flux[sorted_ind[i]] ; } cpl_free(sorted_ind) ; /* Swap the sorted and the input apertures */ swap = cpl_malloc(sizeof(cpl_apertures)); memcpy(swap, sorted_apert, sizeof(cpl_apertures)); memcpy(sorted_apert, in, sizeof(cpl_apertures)); memcpy(in, swap, sizeof(cpl_apertures)); cpl_free(swap); cpl_apertures_delete(sorted_apert); return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief Simple detection of apertures in an image @param in Input image @param sigmas Positive, decreasing sigmas to apply @param pisigma Index of the sigma that was used or unchanged on error @return The detected apertures or NULL on error @see cpl_apertures_extract_sigma() pisigma may be NULL. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if in or sigmas is NULL - CPL_ERROR_DATA_NOT_FOUND if the apertures cannot be detected */ /*----------------------------------------------------------------------------*/ cpl_apertures * cpl_apertures_extract(const cpl_image * in, const cpl_vector * sigmas, int * pisigma) { cpl_errorstate prestate = cpl_errorstate_get(); cpl_apertures * self = NULL; int i; cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(sigmas, CPL_ERROR_NULL_INPUT, NULL); for (i = 0; i < cpl_vector_get_size(sigmas); i++) { const double sigma = cpl_vector_get(sigmas, i); if (sigma <= 0.0) break; /* Apply the detection */ self = cpl_apertures_extract_sigma(in, sigma); if (self != NULL) break; } cpl_ensure(self != NULL, CPL_ERROR_DATA_NOT_FOUND, NULL); /* Recover from any errors set in cpl_apertures_extract_sigma() */ cpl_errorstate_set(prestate); if (pisigma) *pisigma = i; return self; } /*----------------------------------------------------------------------------*/ /** @brief Simple detection of apertures in an image window @param in Input image @param sigmas Positive, decreasing sigmas to apply @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 pisigma Index of the sigma that was used or undefined on error @return The list of detected apertures or NULL on error @see cpl_apertures_extract() @see cpl_image_extract() */ /*----------------------------------------------------------------------------*/ cpl_apertures * cpl_apertures_extract_window( const cpl_image * in, const cpl_vector * sigmas, int llx, int lly, int urx, int ury, int * pisigma) { cpl_apertures * aperts; cpl_image * extracted; int i ; /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL) ; cpl_ensure(sigmas, CPL_ERROR_NULL_INPUT, NULL) ; /* Extract the subwindow - and check entries */ extracted = cpl_image_extract(in, llx, lly, urx, ury); cpl_ensure(extracted, cpl_error_get_code(), NULL) ; /* Get the apertures */ aperts = cpl_apertures_extract(extracted, sigmas, pisigma); cpl_image_delete(extracted) ; cpl_ensure(aperts, cpl_error_get_code(), NULL) ; /* Update the detected positions */ for (i=0 ; inaperts ; i++) { aperts->x[i] += llx - 1 ; aperts->y[i] += lly - 1 ; aperts->xcentroid[i] += llx - 1 ; aperts->ycentroid[i] += lly - 1 ; aperts->left_x[i] += llx - 1 ; aperts->left_y[i] += lly - 1 ; aperts->right_x[i] += llx - 1 ; aperts->right_y[i] += lly - 1 ; aperts->top_x[i] += llx - 1 ; aperts->top_y[i] += lly - 1 ; aperts->bottom_x[i] += llx - 1 ; aperts->bottom_y[i] += lly - 1 ; } return aperts ; } /*----------------------------------------------------------------------------*/ /** @brief Simple apertures detection in an image using a provided sigma @param in Input image @param sigma detection level @return The list of detected apertures or NULL if nothing detected or in error case. The threshold used for the detection is the median plus the average distance to the median times sigma. The input image type can be CPL_TYPE_DOUBLE, CPL_TYPE_FLOAT or CPL_TYPE_INT. 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 sigma is negative or 0.0 - CPL_ERROR_ILLEGAL_OUTPUT if the image cannot be binarised or labelised */ /*----------------------------------------------------------------------------*/ cpl_apertures * cpl_apertures_extract_sigma( const cpl_image * in, double sigma) { double median, med_dist ; double threshold ; cpl_mask * selection ; cpl_matrix * kernel ; cpl_image * labels ; cpl_apertures * aperts ; int nlabels ; /* Check entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL) ; cpl_ensure(sigma>0.0, CPL_ERROR_ILLEGAL_INPUT, NULL) ; /* Compute the threshold */ median = cpl_image_get_median_dev(in, &med_dist) ; threshold = median + sigma * med_dist ; /* Binarise the image */ selection=cpl_mask_threshold_image_create(in, threshold, DBL_MAX); cpl_ensure(selection, CPL_ERROR_ILLEGAL_OUTPUT, NULL) ; /* Apply a morphological opening to remove the single pixel detections */ kernel = cpl_matrix_new(3, 3) ; cpl_matrix_fill(kernel, 1.00) ; if (cpl_mask_opening(selection, kernel) != CPL_ERROR_NONE) { cpl_mask_delete(selection) ; cpl_matrix_delete(kernel) ; cpl_ensure(0, CPL_ERROR_ILLEGAL_OUTPUT, NULL) ; } cpl_matrix_delete(kernel) ; /* Labelise the different detected apertures */ if ((labels = cpl_image_labelise_mask_create(selection, &nlabels))==NULL) { cpl_mask_delete(selection) ; cpl_ensure(0, CPL_ERROR_ILLEGAL_OUTPUT, NULL) ; } cpl_mask_delete(selection) ; /* Nothing detected */ if (nlabels == 0) { cpl_image_delete(labels) ; return NULL ; } /* Create the detected apertures list */ if ((aperts = cpl_apertures_new_from_image(in, labels)) == NULL) { cpl_image_delete(labels) ; cpl_ensure(0, CPL_ERROR_ILLEGAL_OUTPUT, NULL) ; } /* Free and return */ cpl_image_delete(labels) ; return aperts ; } /**@}*/