/* $Id: cpl_image_basic.c,v 1.143 2008/02/07 16:30:20 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: 2008/02/07 16:30:20 $ * $Revision: 1.143 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include #include #include #include #include "cpl_memory.h" #include "cpl_vector.h" #include "cpl_image_basic.h" #include "cpl_image_stats.h" #include "cpl_stats.h" #include "cpl_image_bpm.h" #include "cpl_image_iqe.h" #include "cpl_mask.h" #include "cpl_tools.h" #include "cpl_errorstate.h" #include "cpl_math_const.h" #include "cpl_image_defs.h" /*----------------------------------------------------------------------------- Defines -----------------------------------------------------------------------------*/ #define CPL_IMAGE_BASIC_ASSIGN 1 #define CPL_IMAGE_BASIC_ASSIGN_LOCAL 2 #define CPL_IMAGE_BASIC_THRESHOLD 3 #define CPL_IMAGE_BASIC_ABS 4 #define CPL_IMAGE_BASIC_AVERAGE 5 #define CPL_IMAGE_BASIC_SUBMIN 6 #define CPL_IMAGE_BASIC_EXTRACT 7 #define CPL_IMAGE_BASIC_EXTRACTROW 8 #define CPL_IMAGE_BASIC_EXTRACTCOL 9 #define CPL_IMAGE_BASIC_COLLAPSE 10 #define CPL_IMAGE_BASIC_COLLAPSE_MEDIAN 11 #define CPL_IMAGE_BASIC_ROTATE_INT_LOCAL 12 #define CPL_IMAGE_BASIC_SHIFT_INT_LOCAL 13 #define CPL_IMAGE_BASIC_COPY 14 #define CPL_IMAGE_BASIC_FLIP_LOCAL 15 #define CPL_IMAGE_BASIC_MOVE_PIXELS 16 #define CPL_IMAGE_CHK_POSITIVE 17 #define CPL_IMAGE_CHK_INTEGER 18 #define CPL_IMAGE_CHK_NON_ZERO 19 #define CPL_IMAGE_CHK_NON_NEGATIVE 20 #define CPL_IMAGE_BASIC_OP_SCALAR 21 #define CPL_IMAGE_BASIC_SQRT 22 #define CPL_IMAGE_BASIC_DECLARE 23 #define CPL_IMAGE_BASIC_DEFINE 24 #define CPL_IMAGE_ADDITION(a,b,c) a = (b) + (c) #define CPL_IMAGE_ADDITIONASSIGN(a,b) a += (b) #define CPL_IMAGE_SUBTRACTION(a,b,c) a = (b) - (c) #define CPL_IMAGE_SUBTRACTIONASSIGN(a,b) a -= (b) #define CPL_IMAGE_MULTIPLICATION(a,b,c) a = (b) * (c) #define CPL_IMAGE_MULTIPLICATIONASSIGN(a,b) a *= (b) #define CPL_IMAGE_DIVISION(a,b,c) a = (b) / (c) #define CPL_IMAGE_DIVISIONASSIGN(a,b) a /= (b) /* Caveat: The result is evaluated twice here */ #define CPL_IMAGE_POWASSIGN(a,b) a = pow(a, b) #define CPL_IMAGE_EXPASSIGN(a,b) a = pow(b, a) #define CPL_IMAGE_LOGASSIGN(a,b) a = log(a)/log(b) /*----------------------------------------------------------------------------- Private Function prototypes -----------------------------------------------------------------------------*/ static double cpl_vector_get_noise(const cpl_vector *, int) ; static double cpl_vector_get_fwhm(const cpl_vector *, int, double) ; static cpl_error_code cpl_fft(double *, double *, const unsigned *, int, int); /* Declare the C-type dependent functions */ #define CPL_OPERATION CPL_IMAGE_BASIC_DECLARE #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_image_basic Images basic operation routines * * This module provides functions to handle basic image processing functions * like simple image operations, thresholding, image collapsing, * extraction, etc... * * @par Synopsis: * @code * #include "cpl_image_basic.h" * @endcode */ /*----------------------------------------------------------------------------*/ /**@{*/ /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ #define CPL_OPERATION CPL_IMAGE_BASIC_ASSIGN /*----------------------------------------------------------------------------*/ /** @brief Add two images. @param image1 first operand @param image2 second operand @return 1 newly allocated image or NULL on error Creates a new image, being the result of the operation, and returns it to the caller. The returned image must be deallocated using cpl_image_delete(). The function supports images with different types among CPL_TYPE_INT, CPL_TYPE_FLOAT and CPL_TYPE_DOUBLE. The returned image type is the one of the first passed image. The bad pixels map of the result is the union of the bad pixels maps of the input images. 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 input images have different sizes - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_add_create( const cpl_image * image1, const cpl_image * image2) { #define CPL_OPERATOR CPL_IMAGE_ADDITION #include "cpl_image_basic_body.h" #undef CPL_OPERATOR } /*----------------------------------------------------------------------------*/ /** @brief Subtract two images. @param image1 first operand @param image2 second operand @return 1 newly allocated image or NULL on error @see cpl_image_add_create() */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_subtract_create( const cpl_image * image1, const cpl_image * image2) { #define CPL_OPERATOR CPL_IMAGE_SUBTRACTION #include "cpl_image_basic_body.h" #undef CPL_OPERATOR } /*----------------------------------------------------------------------------*/ /** @brief Multiply two images. @param image1 first operand @param image2 second operand @return 1 newly allocated image or NULL on error @see cpl_image_add_create() */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_multiply_create( const cpl_image * image1, const cpl_image * image2) { #define CPL_OPERATOR CPL_IMAGE_MULTIPLICATION #include "cpl_image_basic_body.h" #undef CPL_OPERATOR } #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Divide two images. @param image1 first operand @param image2 second operand @return 1 newly allocated image or NULL on error @see cpl_image_add_create() @see cpl_image_divide() The result of division with a zero-valued pixel is marked as a bad pixel. */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_divide_create( const cpl_image * image1, const cpl_image * image2) { cpl_image * image_out ; register float * pf1, * pf2, * outpf ; register int * pi1, * pi2, * outpi ; register double * pd1, * pd2, * outpd ; cpl_mask * zeros ; cpl_binary * pzeros ; int i ; cpl_ensure(image1 && image2, CPL_ERROR_NULL_INPUT, NULL); /* Input data images shall have the same sizes */ cpl_ensure(image1->nx == image2->nx && image1->ny == image2->ny, CPL_ERROR_ILLEGAL_INPUT, NULL); /* Create the zeros map */ zeros = cpl_mask_new(image2->nx, image2->ny) ; pzeros = cpl_mask_get_data(zeros) ; /* Switch on the first passed image type */ switch (image1->type) { case CPL_TYPE_INT: image_out = cpl_image_new(image1->nx, image1->ny, CPL_TYPE_INT); pi1 = (int*)image1->pixels ; outpi = (int*)image_out->pixels ; /* Switch on the second passed image type */ switch (image2->type) { case CPL_TYPE_INT: pi2 = (int*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pi2 == 0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpi++ ; pi1++ ; pi2++ ; } else { CPL_IMAGE_DIVISION(*outpi++, *pi1++, *pi2++); pzeros++ ; } } break ; case CPL_TYPE_FLOAT: pf2 = (float*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pf2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpi++ ; pi1++ ; pf2++ ; } else { CPL_IMAGE_DIVISION(*outpi++, *pi1++, *pf2++); pzeros++ ; } } break ; case CPL_TYPE_DOUBLE: pd2 = (double*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pd2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpi++ ; pi1++ ; pd2++ ; } else { CPL_IMAGE_DIVISION(*outpi++, *pi1++, *pd2++); pzeros++ ; } } break ; default: cpl_image_delete(image_out) ; cpl_mask_delete(zeros) ; cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } break ; case CPL_TYPE_FLOAT: image_out = cpl_image_new(image1->nx, image1->ny, CPL_TYPE_FLOAT); pf1 = (float*)image1->pixels ; outpf = (float*)image_out->pixels ; /* Switch on the second passed image type */ switch (image2->type) { case CPL_TYPE_INT: pi2 = (int*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pi2 == 0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpf++ ; pf1++ ; pi2++ ; } else { CPL_IMAGE_DIVISION(*outpf++, *pf1++, *pi2++); pzeros++ ; } } break ; case CPL_TYPE_FLOAT: pf2 = (float*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pf2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpf++ ; pf1++ ; pf2++ ; } else { CPL_IMAGE_DIVISION(*outpf++, *pf1++, *pf2++); pzeros++ ; } } break ; case CPL_TYPE_DOUBLE: pd2 = (double*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pd2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpf++ ; pf1++ ; pd2++ ; } else { CPL_IMAGE_DIVISION(*outpf++, *pf1++, *pd2++); pzeros++ ; } } break ; default: cpl_image_delete(image_out) ; cpl_mask_delete(zeros) ; cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } cpl_tools_add_flops( image_out->nx * image_out->ny ); break ; case CPL_TYPE_DOUBLE: image_out = cpl_image_new(image1->nx, image1->ny, CPL_TYPE_DOUBLE); pd1 = (double*)image1->pixels ; outpd = (double*)image_out->pixels ; /* Switch on the second passed image type */ switch (image2->type) { case CPL_TYPE_INT: pi2 = (int*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pi2 == 0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpd++ ; pd1++ ; pi2++ ; } else { CPL_IMAGE_DIVISION(*outpd++, *pd1++, *pi2++); pzeros++ ; } } break ; case CPL_TYPE_FLOAT: pf2 = (float*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pf2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpd++ ; pd1++ ; pf2++ ; } else { CPL_IMAGE_DIVISION(*outpd++, *pd1++, *pf2++); pzeros++ ; } } break ; case CPL_TYPE_DOUBLE: pd2 = (double*)image2->pixels ; for (i=0 ; i<(image_out->nx * image_out->ny) ; i++) { if (*pd2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; outpd++ ; pd1++ ; pd2++ ; } else { CPL_IMAGE_DIVISION(*outpd++, *pd1++, *pd2++); pzeros++ ; } } break ; default: cpl_image_delete(image_out) ; cpl_mask_delete(zeros) ; cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } cpl_tools_add_flops( image_out->nx * image_out->ny ); break ; default: cpl_mask_delete(zeros) ; cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } /* Handle bad pixels map */ if (image1->bpm == NULL && image2->bpm == NULL) { image_out->bpm = NULL ; } else if (image1->bpm == NULL) { image_out->bpm = cpl_mask_duplicate(image2->bpm) ; } else if (image2->bpm == NULL) { image_out->bpm = cpl_mask_duplicate(image1->bpm) ; } else { image_out->bpm = cpl_mask_duplicate(image1->bpm) ; cpl_mask_or(image_out->bpm, image2->bpm) ; } /* Handle division by zero in the BPM */ if (cpl_mask_count(zeros) > 0) { if (image_out->bpm == NULL) { image_out->bpm = zeros ; } else { cpl_mask_or(image_out->bpm, zeros) ; cpl_mask_delete(zeros) ; } } else cpl_mask_delete(zeros) ; return image_out ; } #define CPL_OPERATION CPL_IMAGE_BASIC_ASSIGN_LOCAL /*----------------------------------------------------------------------------*/ /** @brief Add two images, store the result in the first image. @param im1 first operand. @param im2 second operand. @return the #_cpl_error_code_ or CPL_ERROR_NONE The first input image is modified to contain the results of the operation. The function supports images with different types among CPL_TYPE_INT, CPL_TYPE_FLOAT and CPL_TYPE_DOUBLE. The bad pixels map of the first image becomes the union of the bad pixels maps of the input images. 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 input images have different sizes - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_add( cpl_image * im1, const cpl_image * im2) { #define CPL_OPERATOR CPL_IMAGE_ADDITIONASSIGN #include "cpl_image_basic_body.h" #undef CPL_OPERATOR } /*----------------------------------------------------------------------------*/ /** @brief Subtract two images, store the result in the first image. @param im1 first operand. @param im2 second operand. @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_image_add() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_subtract( cpl_image * im1, const cpl_image * im2) { #define CPL_OPERATOR CPL_IMAGE_SUBTRACTIONASSIGN #include "cpl_image_basic_body.h" #undef CPL_OPERATOR } /*----------------------------------------------------------------------------*/ /** @brief Multiply two images, store the result in the first image. @param im1 first operand. @param im2 second operand. @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_image_add() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_multiply( cpl_image * im1, const cpl_image * im2) { #define CPL_OPERATOR CPL_IMAGE_MULTIPLICATIONASSIGN #include "cpl_image_basic_body.h" #undef CPL_OPERATOR } #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Divide two images, store the result in the first image. @param im1 first operand. @param im2 second operand. @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_image_add() @note The result of division with a zero-valued pixel is marked as a bad pixel. */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_divide( cpl_image * im1, const cpl_image * im2) { register float * pf1, * pf2 ; register int * pi1, * pi2 ; register double * pd1, * pd2 ; cpl_mask * zeros ; cpl_binary * pzeros ; register int i ; cpl_ensure_code(im1, CPL_ERROR_NULL_INPUT); cpl_ensure_code(im2, CPL_ERROR_NULL_INPUT); /* Input data images shall have the same sizes */ cpl_ensure_code(im1->nx == im2->nx, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(im1->ny == im2->ny, CPL_ERROR_ILLEGAL_INPUT); assert( im1->pixels ); assert( im2->pixels ); /* Create the zeros map */ zeros = cpl_mask_new(im2->nx, im2->ny) ; pzeros = cpl_mask_get_data(zeros) ; /* Switch on the first passed image type */ switch (im1->type) { case CPL_TYPE_INT: pi1 = (int*)im1->pixels ; /* Switch on the second passed image type */ switch (im2->type) { case CPL_TYPE_INT: pi2 = (int*)im2->pixels ; for (i=0 ; i<(im1->nx * im1->ny) ; i++) { if (*pi2 == 0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pi1++ ; pi2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pi1++, *pi2++); pzeros++ ; } } break ; case CPL_TYPE_FLOAT: pf2 = (float*)im2->pixels ; for (i=0 ; i<(im1->nx * im1->ny) ; i++) { if (*pf2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pi1++ ; pf2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pi1++, *pf2++); pzeros++ ; } } break ; case CPL_TYPE_DOUBLE: pd2 = (double*)im2->pixels ; for (i=0 ; i<(im1->nx * im1->ny) ; i++) { if (*pd2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pi1++ ; pd2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pi1++, *pd2++); pzeros++ ; } } break ; default: cpl_mask_delete(zeros) ; cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } break ; case CPL_TYPE_FLOAT: pf1 = (float*)im1->pixels ; /* Switch on the second passed image type */ switch (im2->type) { case CPL_TYPE_INT: pi2 = (int*)im2->pixels ; for (i=0 ; i<(im1->nx * im1->ny) ; i++) { if (*pi2 == 0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pf1++ ; pi2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pf1++, *pi2++); pzeros++ ; } } break ; case CPL_TYPE_FLOAT: pf2 = (float*)im2->pixels ; for (i=0 ; i<(im1->nx * im1->ny) ; i++) { if (*pf2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pf1++ ; pf2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pf1++, *pf2++); pzeros++ ; } } break ; case CPL_TYPE_DOUBLE: pd2 = (double*)im2->pixels ; for (i=0 ; i<(im1->nx * im1->ny) ; i++) { if (*pd2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pf1++ ; pd2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pf1++, *pd2++); pzeros++ ; } } break ; default: cpl_mask_delete(zeros) ; cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } cpl_tools_add_flops( im1->nx * im1->ny ); break ; case CPL_TYPE_DOUBLE: pd1 = (double*)im1->pixels ; /* Switch on the second passed image type */ switch (im2->type) { case CPL_TYPE_INT: pi2 = (int*)im2->pixels ; for (i=0 ; i<(im1->nx * im1->ny) ; i++) { if (*pi2 == 0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pd1++ ; pi2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pd1++, *pi2++); pzeros++ ; } } break ; case CPL_TYPE_FLOAT: pf2 = (float*)im2->pixels ; for (i=0 ; i<(im1->nx * im2->ny) ; i++) { if (*pf2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pd1++ ; pf2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pd1++, *pf2++); pzeros++ ; } } break ; case CPL_TYPE_DOUBLE: pd2 = (double*)im2->pixels ; for (i=0 ; i<(im1->nx * im2->ny) ; i++) { if (*pd2 == 0.0) { *pzeros++ = (cpl_binary)CPL_BINARY_1 ; pd1++ ; pd2++ ; } else { CPL_IMAGE_DIVISIONASSIGN(*pd1++, *pd2++); pzeros++ ; } } break ; default: cpl_mask_delete(zeros) ; cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } cpl_tools_add_flops( im1->nx * im1->ny ); break ; default: cpl_mask_delete(zeros) ; cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH) ; } /* Handle bad pixels map */ if (im2->bpm != NULL) { if (im1->bpm == NULL) { im1->bpm = cpl_mask_duplicate(im2->bpm) ; } else { cpl_mask_or(im1->bpm, im2->bpm) ; } } /* Handle division by zero in the BPM */ if (cpl_mask_count(zeros) > 0) { if (im1->bpm == NULL) { im1->bpm = zeros ; } else { cpl_mask_or(im1->bpm, zeros) ; cpl_mask_delete(zeros) ; } } else cpl_mask_delete(zeros) ; return CPL_ERROR_NONE; } #define CPL_OPERATION CPL_IMAGE_BASIC_OP_SCALAR #define CPL_OPERATOR CPL_IMAGE_ADDITIONASSIGN /*----------------------------------------------------------------------------*/ /** @brief Elementwise addition of a scalar to an image @param image Image to be modified in place. @param addend Number to add @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error Modifies the image by adding a number to each of its pixels. The operation is always performed in double precision, with a final cast of the result to the target image type. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT, 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_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_add_scalar( cpl_image * image, double addend) { const double scalar = addend; int i ; /* Error handling */ cpl_ensure_code(image != NULL, CPL_ERROR_NULL_INPUT); assert( image->pixels ); /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE ; } #undef CPL_OPERATOR #define CPL_OPERATOR CPL_IMAGE_SUBTRACTIONASSIGN /*----------------------------------------------------------------------------*/ /** @brief Elementwise subtraction of a scalar from an image @param image Image to be modified in place. @param subtrahend Number to subtract @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error @see cpl_image_add_scalar() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_subtract_scalar( cpl_image * image, double subtrahend) { const double scalar = subtrahend; int i ; /* Error handling */ cpl_ensure_code(image != NULL, CPL_ERROR_NULL_INPUT); assert( image->pixels ); /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE ; } #undef CPL_OPERATOR #define CPL_OPERATOR CPL_IMAGE_MULTIPLICATIONASSIGN /*----------------------------------------------------------------------------*/ /** @brief Elementwise multiplication of an image with a scalar @param image Image to be modified in place. @param factor Number to multiply with @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error @see cpl_image_add_scalar() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_multiply_scalar( cpl_image * image, double factor) { const double scalar = factor; int i ; /* Error handling */ cpl_ensure_code(image != NULL, CPL_ERROR_NULL_INPUT); assert( image->pixels ); /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE ; } #undef CPL_OPERATOR #define CPL_OPERATOR CPL_IMAGE_DIVISIONASSIGN /*----------------------------------------------------------------------------*/ /** @brief Elementwise division of an image with a scalar @param image Image to be modified in place. @param divisor Non-zero number to divide with @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error Modifies the image by dividing each of its pixels with a number. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT, CPL_TYPE_DOUBLE. If divisor is zero, the image is not modified (and an error is returned). Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported - CPL_ERROR_DIVISION_BY_ZERO a division by 0 occurs */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_divide_scalar( cpl_image * image, double divisor) { const double scalar = divisor; int i ; /* Error handling */ cpl_ensure_code(image != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(divisor != 0.0, CPL_ERROR_DIVISION_BY_ZERO); assert( image->pixels ); /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE ; } #undef CPL_OPERATOR #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Compute the elementwise logarithm of the image. @param image Image to be modified in place. @param base Base of the logarithm. @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error Modifies the image by computing the base-scalar logarithm of each of its pixels. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. The base and all the pixels must be positive and the base must be different from 1, or a cpl_error_code will be returned and the image will be left unmodified. 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 base is negative or if the image contains any negative value. - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported - CPL_ERROR_DIVISION_BY_ZERO a division by 0 occurs */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_logarithm( cpl_image * image, double base) { /* Assume that the compiler exploits that log(base) is constant */ const double scalar = base; int i ; unsigned char ok = 1; /* Error handling */ cpl_ensure_code(image != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(base > 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(base != 1, CPL_ERROR_DIVISION_BY_ZERO); assert( image->pixels ); #define CPL_OPERATION CPL_IMAGE_CHK_POSITIVE /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS #undef CPL_OPERATION default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } cpl_ensure_code(ok, CPL_ERROR_ILLEGAL_INPUT); #define CPL_OPERATION CPL_IMAGE_BASIC_OP_SCALAR #define CPL_OPERATOR CPL_IMAGE_LOGASSIGN /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE ; } #undef CPL_OPERATOR #define CPL_OPERATOR CPL_IMAGE_EXPASSIGN /*----------------------------------------------------------------------------*/ /** @brief Compute the elementwise exponential of the image. @param image Image to be modified in place. @param base Base of the exponential. @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error Modifies the image by computing the base-scalar exponential of each of its pixels. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. If the base is zero all pixels must be positive and if the base is negative all pixels must be integer, otherwise a cpl_error_code is returned and the image is unmodified. 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 requirements on the image are not respected - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported - CPL_ERROR_DIVISION_BY_ZERO a division by 0 occurs */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_exponential( cpl_image * image, double base) { const double scalar = base; int i ; /* Error handling */ cpl_ensure_code(image != NULL, CPL_ERROR_NULL_INPUT); assert( image->pixels ); if (base == 0) { unsigned char ok = 1; #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_CHK_POSITIVE /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } #undef CPL_OPERATION cpl_ensure_code(ok, CPL_ERROR_DIVISION_BY_ZERO); } else if (base < 0) { unsigned char ok = 1; #define CPL_OPERATION CPL_IMAGE_CHK_INTEGER /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_OP_SCALAR cpl_ensure_code(ok, CPL_ERROR_ILLEGAL_INPUT); } /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE ; } #undef CPL_OPERATOR #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Compute the elementwise power of the image. @param image Image to be modified in place. @param exponent Scalar exponent. @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error Modifies the image by lifting each of its pixels to exponent. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. If the exponent is non-positive all pixels must be non-zero and if the exponent is non-integer all pixels must be non-negative, otherwise a cpl_error_code is returned and the image is unmodified. 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 requirements on the image are not respected - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported - CPL_ERROR_DIVISION_BY_ZERO a division by 0 occurs */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_power( cpl_image * image, double exponent) { const double scalar = exponent; int i ; /* Error handling */ cpl_ensure_code(image != NULL, CPL_ERROR_NULL_INPUT); if (exponent <= 0) { unsigned char ok = 1; #define CPL_OPERATION CPL_IMAGE_CHK_NON_ZERO /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } #undef CPL_OPERATION cpl_ensure_code(ok, CPL_ERROR_DIVISION_BY_ZERO); } else if (exponent != ceil(exponent)) { unsigned char ok = 1; #define CPL_OPERATION CPL_IMAGE_CHK_NON_NEGATIVE /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } #undef CPL_OPERATION cpl_ensure_code(ok, CPL_ERROR_ILLEGAL_INPUT); } /* Special case with exponent == 0.5 -> use sqrt() for efficiency */ if (scalar == 0.5) { #define CPL_OPERATION CPL_IMAGE_BASIC_SQRT /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } #undef CPL_OPERATION } else { /* General case */ #define CPL_OPERATION CPL_IMAGE_BASIC_OP_SCALAR #define CPL_OPERATOR CPL_IMAGE_POWASSIGN /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } #undef CPL_OPERATOR #undef CPL_OPERATION } return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief Normalise pixels in an image. @param image Image operand. @param mode Normalisation mode. @return CPL_ERROR_NONE, or the relevant #_cpl_error_code_ on error. Normalises an image according to a given criterion. Possible normalisations are: - CPL_NORM_SCALE sets the pixel interval to [0,1]. - CPL_NORM_MEAN sets the mean value to 1. - CPL_NORM_FLUX sets the flux to 1. - CPL_NORM_ABSFLUX sets the absolute flux to 1. 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_normalise( cpl_image * image, cpl_norm mode) { double scale; cpl_ensure_code(image, CPL_ERROR_NULL_INPUT); switch (mode) { case CPL_NORM_SCALE: { cpl_stats * stats = cpl_stats_new_from_image(image, CPL_STATS_MIN | CPL_STATS_MAX); cpl_ensure_code( stats != NULL, cpl_error_get_code()); scale = cpl_stats_get_max(stats) - cpl_stats_get_min(stats) ; if (scale > 0 && cpl_image_subtract_scalar(image, cpl_stats_get_min(stats))) { cpl_stats_delete(stats); cpl_ensure_code(0, cpl_error_get_code()) ; } cpl_stats_delete(stats); break ; } case CPL_NORM_MEAN: { scale = cpl_image_get_mean(image); break ; } case CPL_NORM_FLUX: { scale = cpl_image_get_flux(image); break ; } case CPL_NORM_ABSFLUX: { scale = cpl_image_get_absflux(image); break ; } default: /* This case can only be reached if cpl_norm is extended in error */ assert( 0 ) ; } cpl_ensure_code( !cpl_image_divide_scalar(image, scale), cpl_error_get_code()); return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief Create a new normalised image from an existing image. @param image_in Image operand. @param mode Normalisation mode. @return 1 newly allocated image or NULL on error @see cpl_image_normalise Stores the result in a newly allocated image and returns it. The returned image must be deallocated using cpl_image_delete(). - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_normalise_create( const cpl_image * image_in, cpl_norm mode) { cpl_image * image_out; cpl_ensure(image_in, CPL_ERROR_NULL_INPUT, NULL); image_out = cpl_image_duplicate(image_in); if (cpl_image_normalise(image_out, mode)) { cpl_image_delete(image_out); cpl_ensure(0, cpl_error_get_code(), NULL); } return image_out ; } /*----------------------------------------------------------------------------*/ /** @brief Create a new image by elementwise addition of a scalar to an image @param image Image to add @param addend Number to add @return 1 newly allocated image or NULL in case of an error @see cpl_image_add_scalar Creates a new image, being the result of the operation, and returns it to the caller. The returned image must be deallocated using cpl_image_delete(). The function supports images with different types among CPL_TYPE_INT, CPL_TYPE_FLOAT and CPL_TYPE_DOUBLE. The type of the created image is that of the passed image. */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_add_scalar_create( const cpl_image * image, double addend) { cpl_image * result = cpl_image_duplicate(image); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_add_scalar(result, addend)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } /*----------------------------------------------------------------------------*/ /** @brief Create an image by elementwise subtraction of a scalar from an image @param image Image to be subtracted from @param subtrahend Number to subtract @return 1 newly allocated image or NULL in case of an error @see cpl_image_subtract_scalar @see cpl_image_add_scalar_create */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_subtract_scalar_create( const cpl_image * image, double subtrahend) { cpl_image * result = cpl_image_duplicate(image); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_subtract_scalar(result, subtrahend)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } /*----------------------------------------------------------------------------*/ /** @brief Create a new image by multiplication of a scalar and an image @param image Image to be multiplied @param factor Number to multiply with @return 1 newly allocated image or NULL in case of an error @see cpl_image_multiply_scalar @see cpl_image_add_scalar_create */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_multiply_scalar_create( const cpl_image * image, double factor) { cpl_image * result = cpl_image_duplicate(image); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_multiply_scalar(result, factor)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } /*----------------------------------------------------------------------------*/ /** @brief Create a new image by elementwise division of an image with a scalar @param image Image to divide @param divisor Non-zero number to divide with @return 1 newly allocated image or NULL in case of an error @see cpl_image_divide_scalar @see cpl_image_add_scalar_create */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_divide_scalar_create( const cpl_image * image, double divisor) { cpl_image * result = cpl_image_duplicate(image); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_divide_scalar(result, divisor)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } /*----------------------------------------------------------------------------*/ /** @brief Create a new image by taking the elementwise logarithm of an image @param image Image to take logarithm of @param base Base of the logarithm. @return 1 newly allocated image or NULL in case of an error @see cpl_image_logarithm @see cpl_image_add_scalar_create */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_logarithm_create( const cpl_image * image, double base) { cpl_image * result = cpl_image_duplicate(image); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_logarithm(result, base)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } /*----------------------------------------------------------------------------*/ /** @brief Create a new image by elementwise exponentiation of an image @param image Image to exponentiate @param base Base of the exponential @return 1 newly allocated image or NULL in case of an error @see cpl_image_logarithm @see cpl_image_add_scalar_create */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_exponential_create( const cpl_image * image, double base) { cpl_image * result = cpl_image_duplicate(image); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_exponential(result, base)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } /*----------------------------------------------------------------------------*/ /** @brief Create a new image by elementwise raising of an image to a power @param image Image to raise to a power @param exponent scalar exponent @return 1 newly allocated image or NULL in case of an error @see cpl_image_exponent @see cpl_image_add_scalar_create */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_power_create( const cpl_image * image, double exponent) { cpl_image * result = cpl_image_duplicate(image); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_power(result, exponent)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } #define CPL_OPERATION CPL_IMAGE_BASIC_THRESHOLD /*----------------------------------------------------------------------------*/ /** @brief Threshold an image to a given interval. @param image_in Image to threshold. @param lo_cut Lower bound. @param hi_cut Higher bound. @param assign_lo_cut Value to assign to pixels below low bound. @param assign_hi_cut Value to assign to pixels above high bound. @return the #_cpl_error_code_ or CPL_ERROR_NONE Pixels outside of the provided interval are assigned the given values. Use FLT_MIN and FLT_MAX for floating point images and DBL_MIN and DBL_MAX for double images for the lo_cut and hi_cut to avoid any pixel replacement. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. lo_cut must be smaller than or equal to hi_cut. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported - CPL_ERROR_ILLEGAL_INPUT if lo_cut is greater than hi_cut */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_threshold( cpl_image * image_in, double lo_cut, double hi_cut, double assign_lo_cut, double assign_hi_cut) { int i ; cpl_ensure_code(image_in != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(lo_cut <= hi_cut, CPL_ERROR_ILLEGAL_INPUT); /* Switch on image type */ switch (image_in->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE ; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_ABS /*----------------------------------------------------------------------------*/ /** @brief Take the absolute value of an image. @param image Image to be modified in place @return CPL_ERROR_NONE or the relevant the #_cpl_error_code_ on error Set each pixel to its absolute value. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_abs(cpl_image * image) { int i ; /* Check entries */ cpl_ensure_code(image, CPL_ERROR_NULL_INPUT); /* Switch on image type */ switch (image->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE; } #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Take the absolute value of an image. @param image_in Image operand. @return 1 newly allocated image or NULL on error @see cpl_image_abs For each pixel, out = abs(in). The returned image must be deallocated using cpl_image_delete(). */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_abs_create(const cpl_image * image_in) { cpl_image * result = cpl_image_duplicate(image_in); cpl_ensure(result, cpl_error_get_code(), NULL); if (cpl_image_abs(result)) { cpl_image_delete(result); cpl_ensure(0, cpl_error_get_code(), NULL); } return result; } #define CPL_OPERATION CPL_IMAGE_BASIC_AVERAGE /*----------------------------------------------------------------------------*/ /** @brief Build the average of two images. @param image_1 First image operand. @param image_2 Second image operand. @return 1 newly allocated image or NULL on error Builds the average of two images and returns a newly allocated image, to be deallocated using cpl_image_delete(). The average is arithmetic, i.e. outpix=(pix1+pix2)/2 Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT 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_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_average_create( const cpl_image * image_1, const cpl_image * image_2) { cpl_image * image_out ; int * pii2; float * pfi2; double * pdi2; int i ; /* Check entries */ cpl_ensure(image_1 && image_2, CPL_ERROR_NULL_INPUT, NULL); /* Switch on first passed image type */ switch (image_1->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } /* Handle bad pixels map */ if (image_1->bpm == NULL && image_2->bpm == NULL) { image_out->bpm = NULL ; } else if (image_1->bpm == NULL) { image_out->bpm = cpl_mask_duplicate(image_2->bpm) ; } else if (image_2->bpm == NULL) { image_out->bpm = cpl_mask_duplicate(image_1->bpm) ; } else { image_out->bpm = cpl_mask_duplicate(image_1->bpm) ; cpl_mask_or(image_out->bpm, image_2->bpm) ; } return image_out ; } #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Collapse an image region along its rows or columns. @param self Image to collapse. @param llx lower left x coord. @param lly lower left y coord @param urx upper right x coord @param ury upper right y coord @param direction Collapsing direction. @return a newly allocated image or NULL on error @see cpl_image_collapse_create() llx, lly, urx, ury are the image region coordinates in FITS convention. Those specified bounds are included in the collapsed region. 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 specified window is not valid */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_collapse_window_create(const cpl_image * self, int llx, int lly, int urx, int ury, int direction) { cpl_image * other; /* Switch on image type */ switch (cpl_image_get_type(self)) { case CPL_TYPE_DOUBLE: other = cpl_image_collapse_window_create_double(self, llx, lly, urx, ury, direction); break; case CPL_TYPE_FLOAT: other = cpl_image_collapse_window_create_float(self, llx, lly, urx, ury, direction); break; case CPL_TYPE_INT: other = cpl_image_collapse_window_create_int(self, llx, lly, urx, ury, direction); break; default: /* NULL input will be caught here */ cpl_ensure(0, CPL_ERROR_NULL_INPUT, NULL); } /* Propagate error, if any */ cpl_ensure(other != NULL, cpl_error_get_code(), NULL); return other; } /*----------------------------------------------------------------------------*/ /** @brief Collapse an image along its rows or columns. @param self Image to collapse. @param direction Collapsing direction. @return 1 newly allocated image having 1 row or 1 column or NULL on error Collapsing an image means building up a 1d signal by adding up all pixels on the same row or column. @verbatim Collapse along y: p7 p8 p9 Input image is a 3x3 image containing 9 pixels. p4 p5 p6 The output is an image containing one row with p1 p2 p3 3 pixels A, B, C, where: ---------- A B C A = p1+p4+p7 B = p2+p5+p8 C = p3+p6+p9 If p7 is a bad pixel, A = (p1+p4)*3/2. If p1, p4, p7 are bad, A is flagged as bad. @endverbatim Provide the collapsing direction as an int. Give 0 to collapse along y (sum of rows) and get an image with a single row in output, or give 1 to collapse along x (sym of columns) to get an image with a single column in output. Only the good pixels are collapsed. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. The returned image must be deallocated using cpl_image_delete(). Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_collapse_create(const cpl_image * self, int direction) { cpl_image * other = cpl_image_collapse_window_create(self, 1, 1, cpl_image_get_size_x(self), cpl_image_get_size_y(self), direction); /* Propagate error, if any */ cpl_ensure(other != NULL, cpl_error_get_code(), NULL); return other; } #define CPL_OPERATION CPL_IMAGE_BASIC_COLLAPSE_MEDIAN /*----------------------------------------------------------------------------*/ /** @brief Collapse an image along its rows and columns, with filtering. @param in Input image. @param direction Collapsing direction. @param discard_lo Low rejection parameter. @param discard_hi High rejection parameter. @return 1 newly allocated image having 1 row or 1 column or NULL for error @see cpl_image_collapse_create() Collapsing is done as for cpl_image_collapse_create(). The difference is that pixels are not just summed along rows or columns but sorted, the highest and lowest pixel values are removed, and the median is computed on the remaining values. All pixels are collapsed, even the bad pixels. The bad pixel map of the collapsed image is empty whatever the input is, the bad pixel maps are not used in this function. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. The returned image must be deallocated using cpl_image_delete(). 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 sum of discarded pixels is bigger than the image size in the collapsing direction - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_collapse_median_create( const cpl_image * in, int direction, int discard_lo, int discard_hi) { cpl_image * collapsed ; double * line ; int width ; int i, j ; /* Test entries */ cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(!((direction == 0) && (discard_lo+discard_hi>=in->ny)), CPL_ERROR_ILLEGAL_INPUT, NULL) ; cpl_ensure(!((direction == 1) && (discard_lo+discard_hi>=in->nx)), CPL_ERROR_ILLEGAL_INPUT, NULL) ; /* Switch on image type */ switch (in->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } collapsed->bpm = NULL ; return collapsed ; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_EXTRACT /*----------------------------------------------------------------------------*/ /** @brief Extract a rectangular zone from an image into another image. @param in Input image @param llx Lower left X coordinate @param lly Lower left Y coordinate @param urx Upper right X coordinate @param ury Upper right Y coordinate @return 1 newly allocated image or NULL on error @note The returned image must be deallocated using cpl_image_delete() The input coordinates define the extracted region by giving the coordinates of the lower left and upper right corners (inclusive). Coordinates must be provided in the FITS convention: lower left corner of the image is at (1,1), x increasing from left to right, y increasing from bottom to top. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT 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_ILLEGAL_INPUT if the window coordinates are not valid - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_image * cpl_image_extract(const cpl_image * in, int llx, int lly, int urx, int ury) { cpl_image * self = NULL; cpl_ensure(in != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(llx >= 1, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(llx <= urx, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(urx <= in->nx, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(lly >= 1, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(lly <= ury, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(ury <= in->ny, CPL_ERROR_ILLEGAL_INPUT, NULL); /* Switch on image type */ switch (in->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: /* It is an error in CPL to enter here */ cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL); } cpl_ensure(self != NULL, cpl_error_get_code(), NULL); /* Bad pixels handling */ self->bpm = in->bpm == NULL ? NULL : cpl_mask_extract(in->bpm, llx, lly, urx, ury); return self; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_EXTRACTROW /*----------------------------------------------------------------------------*/ /** @brief Extract a row from an image @param image_in Input image @param pos Position of the row (1 for the bottom one) @return 1 newly allocated cpl_vector or NULL on error Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. The returned vector must be deallocated using cpl_vector_delete(). The bad pixels map is not taken into account in this function. 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 pos is not valid - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_new_from_image_row( const cpl_image * image_in, int pos) { cpl_vector * out ; double * out_data ; int i ; /* Test entries */ cpl_ensure(image_in, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(pos>=1 && pos<=image_in->ny, CPL_ERROR_ILLEGAL_INPUT,NULL); /* Allocate output vector */ out = cpl_vector_new(image_in->nx) ; out_data = cpl_vector_get_data(out) ; /* Switch on image type */ switch (image_in->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_vector_delete(out) ; cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } return out ; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_EXTRACTCOL /*----------------------------------------------------------------------------*/ /** @brief Extract a column from an image @param image_in Input image @param pos Position of the column (1 for the left one) @return 1 newly allocated cpl_vector or NULL on error Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. The returned vector must be deallocated using cpl_vector_delete(). The bad pixels map is not taken into account in this function. 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 pos is not valid - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_new_from_image_column( const cpl_image * image_in, int pos) { cpl_vector * out ; double * out_data ; int i ; /* Check entries */ cpl_ensure(image_in, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(pos>=1 && pos<=image_in->nx, CPL_ERROR_ILLEGAL_INPUT,NULL); /* Allocate output vector */ out = cpl_vector_new(image_in->ny) ; out_data = cpl_vector_get_data(out) ; /* Switch on image type */ switch (image_in->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_vector_delete(out) ; cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ; } return out ; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_ROTATE_INT_LOCAL /*----------------------------------------------------------------------------*/ /** @brief Rotate an image by a multiple of 90 degrees counterclockwise. @param self The image to rotate in place. @param rot The multiple: -1 is a rotation of 90 deg clockwise. @return CPL_ERROR_NONE on success, otherwise the relevant #_cpl_error_code_ @note The dimension of a rectangular image is changed. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. For rotations of +90 or -90 degrees on rectangular non-1D-images, the pixel buffer is temporarily duplicated. rot may be any integer value, its modulo 4 determines the rotation: - -3 to turn 270 degrees clockwise. - -2 to turn 180 degrees clockwise. - -1 to turn 90 degrees clockwise. - 0 to not turn - +1 to turn 90 degrees counterclockwise (same as -3) - +2 to turn 180 degrees counterclockwise (same as -2). - +3 to turn 270 degrees counterclockwise (same as -1). Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_turn(cpl_image * self, int rot) { cpl_ensure_code(self, CPL_ERROR_NULL_INPUT); rot %= 4; if (rot < 0) rot += 4; /* rot is 0, 1, 2 or 3. */ /* Rotate the bad pixel map */ if (rot != 0 && self->bpm != NULL) cpl_mask_turn(self->bpm, rot); if (rot == 1 && self->ny == 1) { /* The image is a single horizontal row, the order of the pixels needs to be reversed */ /* This is equivalent to a 180 degree image rotation */ rot = 2; self->ny = self->nx; self->nx = 1; } else if (rot == 3 && self->nx == 1) { /* The image is a single vertical solumn, the order of the pixels needs to be reversed */ /* This is equivalent to a 180 degree image rotation */ rot = 2; self->nx = self->ny; self->ny = 1; } /* Switch on the image type */ switch (self->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: /* It is a bug in CPL to reach this point */ cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } return CPL_ERROR_NONE; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_SHIFT_INT_LOCAL /*----------------------------------------------------------------------------*/ /** @brief Shift an image by integer offsets @param im the image to shift. @param x_shift x shift @param y_shift y shift @return the #_cpl_error_code_ or CPL_ERROR_NONE This function operates locally on the pixel buffer. The new zones (in the result image) where no new value is computed are set to 0 and flagged as bad pixels. Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT 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_ILLEGAL_INPUT if the requested shift is bigger than the image size - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_shift( cpl_image * im, int x_shift, int y_shift) { int nx, ny ; cpl_image * tmp_im ; cpl_mask * bpm ; int i, j ; /* Check entries */ cpl_ensure(im, CPL_ERROR_NULL_INPUT, CPL_ERROR_NULL_INPUT) ; /* Initialize */ nx = cpl_image_get_size_x(im) ; ny = cpl_image_get_size_y(im) ; cpl_ensure(nx>fabs(x_shift) && ny>fabs(y_shift), CPL_ERROR_ILLEGAL_INPUT, CPL_ERROR_ILLEGAL_INPUT) ; /* Duplicate the input image */ tmp_im = cpl_image_duplicate(im) ; /* Switch on the image type */ switch (im->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_image_delete(tmp_im); cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH); } cpl_image_delete(tmp_im) ; /* Shift the bad pixel map */ if (im->bpm == NULL) bpm = cpl_mask_new(im->nx, im->ny) ; else bpm = cpl_mask_duplicate(im->bpm) ; cpl_mask_shift(bpm, x_shift, y_shift) ; cpl_image_reject_from_mask(im, bpm) ; cpl_mask_delete(bpm) ; return CPL_ERROR_NONE ; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_COPY /*----------------------------------------------------------------------------*/ /** @brief Copy one image into another @param im1 the image in which im2 is inserted @param im2 the inserted image @param xpos the x pixel position in im1 where the lower left pixel of im2 should go (from 1 to the x size of im1) @param ypos the y pixel position in im1 where the lower left pixel of im2 should go (from 1 to the y size of im1) @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error. (xpos, ypos) must be a valid position in im1. If im2 is bigger than the place left in im1, the part that falls outside of im1 is simply ignored, an no error is raised. The bad pixels are inherited from im2 in the concerned im1 zone. The two input images must be of the same type, namely one of CPL_TYPE_INT, CPL_TYPE_FLOAT, CPL_TYPE_DOUBLE. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported or if the input images are of different types - CPL_ERROR_ACCESS_OUT_OF_RANGE if xpos or ypos are outside the specified range */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_copy( cpl_image * im1, const cpl_image * im2, int xpos, int ypos) { const int nx1 = cpl_image_get_size_x(im1) ; const int ny1 = cpl_image_get_size_y(im1) ; const int nx2 = cpl_image_get_size_x(im2) ; const int ny2 = cpl_image_get_size_y(im2) ; cpl_mask * bpm1 ; cpl_mask * bpm2 ; int urx, ury ; int im1_pos, im2_pos ; int i, j ; /* Check entries */ cpl_ensure_code(nx1 > 0, cpl_error_get_code()); cpl_ensure_code(ny1 > 0, cpl_error_get_code()); cpl_ensure_code(nx2 > 0, cpl_error_get_code()); cpl_ensure_code(ny2 > 0, cpl_error_get_code()); cpl_ensure_code(xpos >= 1, CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_ensure_code(xpos <= nx1, CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_ensure_code(ypos >= 1, CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_ensure_code(ypos <= ny1, CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_ensure_code(im1->type == im2->type, CPL_ERROR_TYPE_MISMATCH); /* Define the zone to modify in im1: xpos, ypos, urx, ury */ if (xpos+nx2-1 >= nx1) urx = nx1; else urx = xpos+nx2-1; if (ypos+ny2-1 >= ny1) ury = ny1; else ury = ypos+ny2-1; /* Switch on the image type */ switch (im1->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH) ; } /* Handle the bad pixels */ if (im1->bpm != NULL || im2->bpm != NULL) { if (im1->bpm == NULL) bpm1 = cpl_mask_new(im1->nx, im1->ny) ; else bpm1 = cpl_mask_duplicate(im1->bpm) ; if (im2->bpm == NULL) bpm2 = cpl_mask_new(im2->nx, im2->ny) ; else bpm2 = cpl_mask_duplicate(im2->bpm) ; cpl_mask_copy(bpm1, bpm2, xpos, ypos) ; cpl_mask_delete(bpm2) ; cpl_image_reject_from_mask(im1, bpm1) ; cpl_mask_delete(bpm1) ; } return CPL_ERROR_NONE ; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_FLIP_LOCAL /*----------------------------------------------------------------------------*/ /** @brief Flip an image on a given miror line. @param im the image to flip. @param angle mirror line in polar coord. is theta = (PI/4) * angle @return the #_cpl_error_code_ or CPL_ERROR_NONE This function operates locally on the pixel buffer. angle can take one of the following values: - 0 (theta=0) to flip the image around the horizontal - 1 (theta=pi/4) to flip the image around y=x - 2 (theta=pi/2) to flip the image around the vertical - 3 (theta=3pi/4) to flip the image around y=-x Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT 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_ILLEGAL_INPUT if the angle is different from the allowed values - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_flip( cpl_image * im, int angle) { int nx, ny ; cpl_image * tmp_im ; int i, j ; /* Check entries */ cpl_ensure_code(im, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(angle==0 || angle==1 || angle==2 || angle==3, CPL_ERROR_ILLEGAL_INPUT) ; /* Initialise */ nx = cpl_image_get_size_x(im) ; ny = cpl_image_get_size_y(im) ; /* Duplicate the input image */ /* FIXME: Avoid this by swapping pairs of pixels */ tmp_im = cpl_image_duplicate(im) ; /* Switch on the image type */ switch (im->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, CPL_ERROR_TYPE_MISMATCH) ; } cpl_image_delete(tmp_im) ; /* Flip the bad pixel map */ if (im->bpm != NULL) cpl_mask_flip(im->bpm, angle) ; return CPL_ERROR_NONE ; } #undef CPL_OPERATION #define CPL_OPERATION CPL_IMAGE_BASIC_MOVE_PIXELS /*----------------------------------------------------------------------------*/ /** @brief Reorganize the pixels in an image @param im the image to reorganize @param nb_cut the number of cut in x and y @param new_pos array with the nb_cut^2 new positions @return the #_cpl_error_code_ or CPL_ERROR_NONE nb_cut^2 defines in how many tiles the images will be moved. Each tile will then be moved to an other place defined in new_pos. 1 will leave the image unchanged, 2 is used to move the quadrants, etc.. new_pos contains nb_cut^2 values between 1 and nb_cut^2. The zones positions are counted from the lower left part of the image. It is not allowed to move two tiles to the same position (the relation between th new tiles positions and the initial position is bijective !). Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. The image x and y sizes have to be multiples of nb_cut. @verbatim Example: 16 17 18 6 5 4 13 14 15 3 2 1 10 11 12 ----> 12 11 10 7 8 9 9 8 7 4 5 6 18 17 16 1 2 3 15 14 13 image 3x6 cpl_image_move(image, 3, new_pos) ; with new_pos = {9,8,7,6,5,4,3,2,1} ; @endverbatim The bad pixels are moved accordingly. 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 nb_cut is not strictly positive or cannot divide one of the image sizes or if the new_pos array specifies to move two tiles to the same position. - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_move( cpl_image * im, int nb_cut, const int * new_pos) { int test_sum ; int tile_sz_x, tile_sz_y ; int tile_x, tile_y ; int npos, opos ; int i, j, k, l ; /* Check entries */ cpl_ensure_code(im, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(new_pos, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(nb_cut > 0, CPL_ERROR_ILLEGAL_INPUT) ; cpl_ensure_code(im->nx % nb_cut == 0, CPL_ERROR_ILLEGAL_INPUT) ; cpl_ensure_code(im->ny % nb_cut == 0, CPL_ERROR_ILLEGAL_INPUT) ; /* Test that new_pos takes all values between 1 and nb_cut*nb_cut */ /* The test here is not strict, but should be sufficient */ test_sum = 0 ; for (i=0 ; inx / nb_cut ; tile_sz_y = im->ny / nb_cut ; /* Switch on the image type */ switch (im->type) { #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS default: cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH) ; } /* Handle the bad pixels */ if (im->bpm != NULL) cpl_mask_move(im->bpm, nb_cut, new_pos) ; return CPL_ERROR_NONE ; } #undef CPL_OPERATION /*----------------------------------------------------------------------------*/ /** @brief Apply a gaussian fit on an image sub window @param im the input image @param xpos the x position of the center (1 for the first pixel) @param ypos the y position of the center (1 for the first pixel) @param size the window size in pixels, at least 4 @param norm the norm of the gaussian or NULL @param xcen the x center of the gaussian or NULL @param ycen the y center of the gaussian or NULL @param sig_x the sigma in x of the gaussian or NULL @param sig_y the sigma in y of the gaussian or NULL @param fwhm_x the FHHM in x or NULL @param fwhm_y the FHHM in y or NULL @return the #_cpl_error_code_ or CPL_ERROR_NONE @see cpl_image_iqe() The computed norm, xcen, ycen, sig_x, sig_y coefficients are defining the gaussian: f(x, y) = (norm/(2*pi*sig_x*sig_y)) * exp(-(x-xcen)^2/(2*sig_x^2)) * exp(-(y-ycen)^2/(2*sig_y^2)) sig_x and sig_y are derived from fwhm_x and fwhm_y with: fwhm = 2 * sqrt(2*ln(2)) * sigma Images can be CPL_TYPE_INT, CPL_TYPE_FLOAT or CPL_TYPE_DOUBLE. The cpl_image_iqe() function is directly called by this function. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if the input image pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if xpos, ypos or size or illegal - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_fit_gaussian( const cpl_image * im, int xpos, int ypos, int size, double * norm, double * xcen, double * ycen, double * sig_x, double * sig_y, double * fwhm_x, double * fwhm_y) { const cpl_image * im_local; int llx, lly, urx, ury ; cpl_bivector * stats ; double * pstats ; /* Check entries */ cpl_ensure_code(im != NULL, CPL_ERROR_NULL_INPUT) ; cpl_ensure_code(xpos >= 1, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(ypos >= 1, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(xpos <= im->nx, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(ypos <= im->ny, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(size >= 1, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(size < im->nx, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(size < im->ny, CPL_ERROR_ILLEGAL_INPUT); /* Extraction zone */ llx = xpos - (int)(size/2) ; lly = ypos - (int)(size/2) ; urx = xpos + (int)(size/2) ; ury = ypos + (int)(size/2) ; if (llx < 1) llx = 1 ; if (lly < 1) lly = 1 ; if (urx > im->nx) urx = im->nx ; if (ury > im->ny) ury = im->ny ; /* Convert the image to FLOAT, if needed */ im_local = cpl_image_get_type(im) == CPL_TYPE_FLOAT ? im : cpl_image_cast(im, CPL_TYPE_FLOAT); /* Call cpl_image_iqe */ stats = cpl_image_iqe(im_local, llx, lly, urx, ury); if (im_local != im) cpl_image_delete((cpl_image*)im_local); cpl_ensure_code(stats != NULL, cpl_error_get_code()); /* Write the results */ pstats = cpl_bivector_get_x_data(stats) ; if (xcen) *xcen = pstats[0] ; if (ycen) *ycen = pstats[1] ; if (fwhm_x) *fwhm_x = pstats[2] ; if (fwhm_y) *fwhm_y = pstats[3] ; if (sig_x) *sig_x = pstats[2] / CPL_MATH_FWHM_SIG ; if (sig_y) *sig_y = pstats[3] / CPL_MATH_FWHM_SIG ; if (norm) *norm = pstats[5] * 2 * CPL_MATH_PI * (pstats[2]*pstats[3]) / (CPL_MATH_FWHM_SIG*CPL_MATH_FWHM_SIG) ; cpl_bivector_delete(stats) ; return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @brief Compute FWHM values in x and y for an object @param in the input image @param xpos the x position of the object (1 for the first pixel) @param ypos the y position of the object (1 for the first pixel) @param fwhm_x the computed FWHM in x or -1 on error @param fwhm_y the computed FWHM in y or -1 on error @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ This function uses a basic method: start from the center of the object and go away until the half maximum value is reached in x and y. For the FWHM in x (resp. y) to be computed, the image size in the x (resp. y) direction should be at least of 5 pixels. If for any reason, one of the FHWMs cannot be computed, its returned value is -1.0, but an error is not necessarily raised. For example, if a 4 column image is passed, the fwhm_x would be -1.0, the fwhm_y would be correctly computed, and no error would be raised. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL - CPL_ERROR_DATA_NOT_FOUND if (xpos, ypos) specifies a rejected pixel or a pixel with a non-positive value - CPL_ERROR_ACCESS_OUT_OF_RANGE if xpos or ypos is outside the image size range */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_get_fwhm( const cpl_image * in, int xpos, int ypos, double * fwhm_x, double * fwhm_y) { double half_max; double thres; const int minimum_size = 5; int is_rejected; /* Check entries - and initialize *fwhm_{x,y} */ if (fwhm_y != NULL) *fwhm_y = -1; cpl_ensure_code(fwhm_x, CPL_ERROR_NULL_INPUT); *fwhm_x = -1; cpl_ensure_code(fwhm_y, CPL_ERROR_NULL_INPUT); /* This call will check the validity of image, xpos and ypos */ half_max = 0.5 * cpl_image_get(in, xpos, ypos, &is_rejected); cpl_ensure_code(is_rejected >= 0, cpl_error_get_code()); cpl_ensure_code(!is_rejected, CPL_ERROR_DATA_NOT_FOUND); cpl_ensure_code(half_max > 0, CPL_ERROR_DATA_NOT_FOUND); /* FWHM in x */ if (in->nx >= minimum_size) { cpl_errorstate pstate; /* Extract the vector centered on the maximum */ cpl_vector * row = cpl_vector_new_from_image_row(in, ypos); /* If an error happened, update its location */ cpl_ensure_code(row, cpl_error_get_code()); pstate = cpl_errorstate_get(); /* Find out threshold */ thres = cpl_vector_get_noise(row, xpos); /* Compute the FWHM */ if (cpl_errorstate_is_equal(pstate)) *fwhm_x = cpl_vector_get_fwhm(row, xpos, half_max + thres * 0.5); cpl_vector_delete(row); /* Propagate the error, if any */ cpl_ensure_code(cpl_errorstate_is_equal(pstate), cpl_error_get_code()); } /* FWHM in y */ if (in->ny >= minimum_size) { cpl_errorstate pstate; /* Extract the vector centered on the maximum */ cpl_vector * col = cpl_vector_new_from_image_column(in, xpos); /* If an error happened, update its location */ cpl_ensure_code(col, cpl_error_get_code()); pstate = cpl_errorstate_get(); /* Find out threshold */ thres = cpl_vector_get_noise(col, ypos); /* Compute the FWHM */ if (cpl_errorstate_is_equal(pstate)) *fwhm_y = cpl_vector_get_fwhm(col, ypos, half_max + thres * 0.5); cpl_vector_delete(col); /* Propagate the error, if any */ cpl_ensure_code(cpl_errorstate_is_equal(pstate), cpl_error_get_code()); } return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Fast Fourier Transform an image. @param real The image real part to be transformed in place @param imaginary The image imaginary part to be transformed in place @param mode The desired FFT options (combined with bitwise or) @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ on error The input images must be of double type. If the second passed image is NULL, the resulting imaginary part cannot be returned. This can be useful if the input is real and the output is known to also be real. But if the output has a significant imaginary part, you might want to pass a 0-valued image as the second parameter. Any rejected pixel us used as if it were a good pixel. The image must be square with a size that is a power of two. These are the supported FFT modes: CPL_FFT_DEFAULT: Default, forward FFT transform CPL_FFT_INVERSE: Inverse FFT transform CPL_FFT_UNNORMALIZED: Do not normalize (with N*N for N-by-N image) on inverse. Has no effect on forward transform. CPL_FFT_SWAP_HALVES: Swap the four quadrants of the result image. 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 is not square or if the image size is not a power of 2. - CPL_ERROR_INVALID_TYPE if mode is 1, e.g. due to a logical or (||) of the allowed FFT options. - CPL_ERROR_UNSUPPORTED_MODE if mode is otherwise different from the allowed FFT options. - CPL_ERROR_TYPE_MISMATCH if the passed image type is not supported */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_image_fft( cpl_image * real, cpl_image * imaginary, unsigned mode) { unsigned dim[2]; double * imag_part ; /* Check entries */ cpl_ensure_code(real, CPL_ERROR_NULL_INPUT); cpl_ensure_code(mode != 1, CPL_ERROR_INVALID_TYPE); cpl_ensure_code( mode <= (CPL_FFT_INVERSE|CPL_FFT_UNNORMALIZED|CPL_FFT_SWAP_HALVES), CPL_ERROR_UNSUPPORTED_MODE); cpl_ensure_code(real->nx==real->ny, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(real->type==CPL_TYPE_DOUBLE,CPL_ERROR_TYPE_MISMATCH); cpl_ensure_code(cpl_tools_is_power_of_2(real->nx)>=0, CPL_ERROR_ILLEGAL_INPUT); if (imaginary != NULL) { cpl_ensure_code(real->nx==imaginary->nx,CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(real->ny==imaginary->ny,CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(imaginary->type==CPL_TYPE_DOUBLE, CPL_ERROR_TYPE_MISMATCH); } /* Initialize */ dim[0] = dim[1] = real->nx; if (imaginary == NULL) { /* Create the imaginary part and set it to 0 */ imag_part = cpl_calloc(dim[0]*dim[1], sizeof(double)) ; } else { /* Put the input imaginary part in a local object */ imag_part = imaginary->pixels ; } /* APPLY THE FFT HERE */ cpl_ensure_code(!cpl_fft(real->pixels, imag_part, dim, 2, (mode & CPL_FFT_INVERSE) ? -1 : 1),cpl_error_get_code()); /* Free the imaginary part result in the input image */ if (imaginary == NULL) { cpl_free(imag_part) ; } /* Normalize on the inverse transform */ if (!(mode & CPL_FFT_UNNORMALIZED) && (mode & CPL_FFT_INVERSE)) { cpl_ensure_code(!cpl_image_divide_scalar(real, dim[0]*dim[1]), cpl_error_get_code()); if (imaginary != NULL) { cpl_ensure_code(!cpl_image_divide_scalar(imaginary, dim[0]*dim[1]), cpl_error_get_code()); } } /* Swap halves in both dimensions */ if (mode & CPL_FFT_SWAP_HALVES) { const int new_pos[] = {4, 3, 2, 1}; cpl_ensure_code(!cpl_image_move(real, 2, new_pos), cpl_error_get_code()); if (imaginary != NULL) { cpl_ensure_code(!cpl_image_move(imaginary, 2, new_pos), cpl_error_get_code()); } } return CPL_ERROR_NONE; } /**@}*/ /*----------------------------------------------------------------------------*/ /** @brief N-dimensional FFT. @param real N-dimensional data set stored in 1d (real part). @param imag N-dimensional data set stored in 1d (imaginary part). @param nn Dimensions of the set. @param ndim How many dimensions this set has. @param isign Transform direction. @return the #_cpl_error_code_ or CPL_ERROR_NONE This routine is a public domain FFT. See extract of Usenet article below. Found on http://www.tu-chemnitz.de/~arndt/joerg.html @verbatim From: alee@tybalt.caltech.edu (Andrew Lee) Newsgroups: comp.sources.misc Subject: N-dimensional, Radix 2 FFT Routine Date: 17 Jul 87 22:26:29 GMT Approved: allbery@ncoast.UUCP X-Archive: comp.sources.misc/8707/48 [..] Now for the usage (finally): real[] and imag[] are the array of complex numbers to be transformed, nn[] is the array giving the dimensions (I mean size) of the array, ndim is the number of dimensions of the array, and isign is +1 for a forward transform, and -1 for an inverse transform. real[], imag[] and nn[] are stored in the "natural" order for C: nn[0] gives the number of elements along the leftmost index, nn[ndim - 1] gives the number of elements along the rightmost index. Additional notes: The routine does NO NORMALIZATION, so if you do a forward, and then an inverse transform on an array, the result will be identical to the original array MULTIPLIED BY THE NUMBER OF ELEMENTS IN THE ARRAY. Also, of course, the dimensions of real[] and imag[] must all be powers of 2. @endverbatim Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if (one of) the input pointer(s) is NULL */ /*----------------------------------------------------------------------------*/ static cpl_error_code cpl_fft( double * real, double * imag, const unsigned * nn, int ndim, int isign) { int idim; unsigned i1, i2rev, i3rev, ibit; unsigned ip2, ifp1, ifp2, k2, n; unsigned nprev = 1, ntot = 1; register unsigned i2, i3; double theta; double w_r, w_i, wp_r, wp_i; double wtemp; double temp_r, temp_i, wt_r, wt_i; double t1, t2; /* Check entries */ cpl_ensure_code(real, CPL_ERROR_NULL_INPUT); cpl_ensure_code(imag, CPL_ERROR_NULL_INPUT); cpl_ensure_code(nn, CPL_ERROR_NULL_INPUT); for (idim = 0; idim < ndim; ++idim) cpl_ensure_code(cpl_tools_is_power_of_2(nn[idim])>=0, CPL_ERROR_ILLEGAL_INPUT); /* Compute total number of complex values */ for (idim = 0; idim < ndim; ++idim) ntot *= nn[idim]; for (idim = ndim - 1; idim >= 0; --idim) { n = nn[idim]; ip2 = nprev * n; /* Unit step for next dimension */ i2rev = 0; /* Bit reversed i2 */ /* This is the bit reversal section of the routine */ /* Loop over current dimension */ for (i2 = 0; i2 < ip2; i2 += nprev) { if (i2 < i2rev) /* Loop over lower dimensions */ for (i1 = i2; i1 < i2 + nprev; ++i1) /* Loop over higher dimensions */ for (i3 = i1; i3 < ntot; i3 += ip2) { i3rev = i3 + i2rev - i2; temp_r = real[i3]; temp_i = imag[i3]; real[i3] = real[i3rev]; imag[i3] = imag[i3rev]; real[i3rev] = temp_r; imag[i3rev] = temp_i; } ibit = ip2; /* Increment from high end of i2rev to low */ do { ibit >>= 1; i2rev ^= ibit; } while (ibit >= nprev && !(ibit & i2rev)); } /* Here begins the Danielson-Lanczos section of the routine */ /* Loop over step sizes */ for (ifp1 = nprev; ifp1 < ip2; ifp1 <<= 1) { ifp2 = ifp1 << 1; /* Initialize for the trig. recurrence */ theta = isign * CPL_MATH_2PI / (ifp2 / nprev); wp_r = sin(0.5 * theta); wp_r *= -2.0 * wp_r; wp_i = sin(theta); w_r = 1.0; w_i = 0.0; /* Loop by unit step in current dimension */ for (i3 = 0; i3 < ifp1; i3 += nprev) { /* Loop over lower dimensions */ for (i1 = i3; i1 < i3 + nprev; ++i1) /* Loop over higher dimensions */ for (i2 = i1; i2 < ntot; i2 += ifp2) { /* Danielson-Lanczos formula */ k2 = i2 + ifp1; wt_r = real[k2]; wt_i = imag[k2]; /* Complex multiply using 3 real multiplies. */ real[k2] = real[i2] - (temp_r = (t1 = w_r * wt_r) - (t2 = w_i * wt_i)); imag[k2] = imag[i2] - (temp_i = (w_r + w_i) * (wt_r + wt_i) - t1 - t2); real[i2] += temp_r; imag[i2] += temp_i; } /* Trigonometric recurrence */ wtemp = w_r; /* Complex multiply using 3 real multiplies. */ w_r += (t1 = w_r * wp_r) - (t2 = w_i * wp_i); w_i += (wtemp + w_i) * (wp_r + wp_i) - t1 - t2; } } nprev *= n; } return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Compute the noise around a peak in a vector @param vec the input cpl_vector @param pos the peak position (from 1 to the vector size) @return the noise value. The passed cpl_vector object must have at least two elements. 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 */ /*----------------------------------------------------------------------------*/ static double cpl_vector_get_noise( const cpl_vector * vec, int pos) { int nelem ; cpl_vector * smooth_vec ; const double * vector_data ; double noise_left, noise_right ; int i ; /* Check entries */ cpl_ensure(vec, CPL_ERROR_NULL_INPUT, -1.0) ; nelem = cpl_vector_get_size(vec) ; cpl_ensure(pos >= 1, CPL_ERROR_ILLEGAL_INPUT, -1.0) ; cpl_ensure(pos <= nelem, CPL_ERROR_ILLEGAL_INPUT, -1.0) ; cpl_ensure(nelem > 1, CPL_ERROR_ILLEGAL_INPUT, -1.0) ; /* Smooth out the array to be less sensitive to noise */ smooth_vec = cpl_vector_filter_lowpass_create(vec, CPL_LOWPASS_LINEAR, 1) ; /* If that fails, use the unsmoothed vector */ vector_data = cpl_vector_get_data_const(smooth_vec ? smooth_vec : vec); /* Find noise level on the left side of the peak. */ i = pos - 1 ; while (i > 0) { if (vector_data[i] > vector_data[i-1]) i--; else break; } noise_left = vector_data[i] ; /* Find noise level on the right side of the peak */ i = pos - 1 ; while (i < nelem-1) { if (vector_data[i] > vector_data[i+1]) i++ ; else break ; } noise_right = vector_data[i] ; cpl_vector_delete(smooth_vec) ; return 0.5 * (noise_left + noise_right); } /*----------------------------------------------------------------------------*/ /** @brief Compute the FWHM of an object in a cpl_vector @param vec the input cpl_vector @param pos the object position (from 1 to the vector size) @param half_max the half maximum value @return the FWHM or a negative value on error @note The return value may be -1 with no error condition Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if pos is less than 1 or greater than the vec-length */ /*----------------------------------------------------------------------------*/ static double cpl_vector_get_fwhm( const cpl_vector * vec, int pos, double half_max) { const double * vec_data; int nelem; double x_left, x_right; double y_1, y_2; int i; /* Check entries */ cpl_ensure(vec, CPL_ERROR_NULL_INPUT, -1.0) ; nelem = cpl_vector_get_size(vec) ; cpl_ensure(pos >= 1, CPL_ERROR_ILLEGAL_INPUT, -1.0) ; cpl_ensure(pos <= nelem, CPL_ERROR_ILLEGAL_INPUT, -1.0) ; vec_data = cpl_vector_get_data_const(vec) ; /* Object may be too noisy - or strange in some other way */ if (vec_data[pos - 1] <= half_max) return -1.0; /* Find first pair of values, y(i) <= half_max < y(i+1) on the left of the maximum */ i = pos - 1; while ((vec_data[i] > half_max) && (i > 0)) i--; if (vec_data[i] > half_max) return -1.0; /* y_1 could not be found */ y_1 = vec_data[i] ; y_2 = vec_data[i+1] ; /* assert ( y_1 <= half_max && half_max < y_2 ); */ /* Assume linearity between y_1 and y_2 */ x_left = i + (half_max-y_1) / (y_2-y_1) ; /* assert( x_left >= i ); */ /* Find first pair of values, y(i-1) > half_max >= y(i) on the right of the maximum */ i = pos - 1 ; while ((vec_data[i] > half_max) && (i < nelem-1)) i++; if (vec_data[i] > half_max) return -1.0; /* y_2 could not be found */ y_1 = vec_data[i-1] ; y_2 = vec_data[i] ; /* assert( y_1 > half_max && half_max >= y_2 ); */ /* Assume linearity between y_1 and y_2 */ x_right = i + (half_max-y_2) / (y_2-y_1) ; /* assert( x_right < i ); */ if (x_right < x_left || x_right - x_left > FLT_MAX) return -1; return x_right - x_left ; } #define CPL_OPERATION CPL_IMAGE_BASIC_DEFINE /* Define the C-type dependent functions */ #define CPL_CLASS CPL_CLASS_DOUBLE #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_FLOAT #include "cpl_image_basic_body.h" #undef CPL_CLASS #define CPL_CLASS CPL_CLASS_INT #include "cpl_image_basic_body.h" #undef CPL_CLASS #undef CPL_OPERATION