/* $Id: cpl_bivector.c,v 1.21 2007/07/17 13:50:25 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/07/17 13:50:25 $ * $Revision: 1.21 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include #include "cpl_memory.h" #include "cpl_bivector.h" #include "cpl_tools.h" #include "cpl_error.h" /*----------------------------------------------------------------------------- Define -----------------------------------------------------------------------------*/ #define HALF_CENTROID_DOMAIN 5 /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_bivector Bi-vector object * * This module provides functions to handle @em cpl_bivector. * * A @em cpl_bivector is composed of two vectors of the same size. It can * be used to store 1d functions, with the x and y positions of the * samples, offsets in x and y or simply positions in an image. * * This module provides the possibility to compute a local maximum, to * cross-correlate two signals, etc... * * @par Synopsis: * @code * #include "cpl_bivector.h" * @endcode */ /*----------------------------------------------------------------------------*/ /**@{*/ /*----------------------------------------------------------------------------- Type definition -----------------------------------------------------------------------------*/ struct _cpl_bivector_ { cpl_vector * x; cpl_vector * y; }; /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /** @brief Create a new cpl_bivector @param n Positive number of points @return 1 newly allocated cpl_bivector or NULL on error The returned object must be deallocated using cpl_bivector_delete() or cpl_bivector_unwrap_vectors(), provided the two cpl_vectors are deallocated separately. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_ILLEGAL_INPUT if n is < 1. */ /*----------------------------------------------------------------------------*/ cpl_bivector * cpl_bivector_new(int n) { cpl_bivector * f; /* Test input */ cpl_ensure(n > 0, CPL_ERROR_ILLEGAL_INPUT, NULL); /* Allocate memory */ f = cpl_malloc(sizeof(cpl_bivector)); f->x = cpl_vector_new(n); f->y = cpl_vector_new(n); return f; } /*----------------------------------------------------------------------------*/ /** @brief Create a new cpl_bivector from two cpl_vectors @param x the x cpl_vector @param y the y cpl_vector @return 1 cpl_bivector or NULL on error @note The input cpl_vectors must have identical sizes. Afterwards one of those two vectors may be resized, which will corrupt the bivector. Such a corrupted bivector should not be used any more, but rather deallocated, using cpl_bivector_unwrap_vectors() or cpl_bivector_delete(). The returned object must be deallocated using cpl_bivector_delete() or with cpl_bivector_unwrap_vectors(), provided the two cpl_vectors are deallocated separately. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if the input vectors have different sizes */ /*----------------------------------------------------------------------------*/ cpl_bivector * cpl_bivector_wrap_vectors( cpl_vector * x, cpl_vector * y) { cpl_bivector * f; /* Test input */ cpl_ensure(x, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(y, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(cpl_vector_get_size(x)==cpl_vector_get_size(y), CPL_ERROR_ILLEGAL_INPUT, NULL); /* Allocate memory */ f = cpl_malloc(sizeof(cpl_bivector)); f->x = x; f->y = y; return f; } /*----------------------------------------------------------------------------*/ /** @brief Duplicate a cpl_bivector @param in cpl_bivector to duplicate @return 1 newly allocated cpl_bivector or NULL on error The returned object must be deallocated using cpl_bivector_delete() Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if the input bivector contains vectors of different sizes */ /*----------------------------------------------------------------------------*/ cpl_bivector * cpl_bivector_duplicate(const cpl_bivector * in) { cpl_bivector * out; cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(cpl_vector_get_size(in->x)==cpl_vector_get_size(in->y), CPL_ERROR_ILLEGAL_INPUT, NULL); out = cpl_malloc(sizeof(cpl_bivector)); out->x = cpl_vector_duplicate(in->x); out->y = cpl_vector_duplicate(in->y); return out; } /*----------------------------------------------------------------------------*/ /** @brief Delete a cpl_bivector @param f cpl_bivector to delete @return void */ /*----------------------------------------------------------------------------*/ void cpl_bivector_delete(cpl_bivector * f) { if (f == NULL) return; cpl_vector_delete(f->x); cpl_vector_delete(f->y); cpl_free(f); return; } /*----------------------------------------------------------------------------*/ /** @brief Free memory associated to a cpl_bivector, excluding the two vectors. @param f cpl_bivector to delete @return void @see cpl_bivector_wrap_vectors */ /*----------------------------------------------------------------------------*/ void cpl_bivector_unwrap_vectors(cpl_bivector * f) { cpl_free(f); return; } /*----------------------------------------------------------------------------*/ /** @brief Dump a cpl_bivector as ASCII to a stream @param f Input cpl_bivector to dump or NULL @param stream Output stream, accepts @c stdout or @c stderr or NULL @return void Comment lines start with the hash character. stream may be NULL in which case @c stdout is used. @note In principle a cpl_bivector can be saved using cpl_bivector_dump() and re-read using cpl_bivector_read(). This will however introduce significant precision loss due to the limited accuracy of the ASCII representation. */ /*----------------------------------------------------------------------------*/ void cpl_bivector_dump(const cpl_bivector * f, FILE * stream) { const double * data_x; const double * data_y; int i; if (stream == NULL) stream = stdout; if (f == NULL) { fprintf(stream, "#NULL bivector\n"); return; } if (cpl_vector_get_size(f->x) != cpl_vector_get_size(f->y)) return; data_x = cpl_vector_get_data_const(f->x); data_y = cpl_vector_get_data_const(f->y); fprintf(stream, "#--- Bi-vector ---\n"); fprintf(stream, "# X Y \n"); for (i=0 ; i < cpl_bivector_get_size(f) ; i++) fprintf(stream, "%g\t%g\n", data_x[i], data_y[i]); fprintf(stream, "#-----------------------\n"); return; } /*----------------------------------------------------------------------------*/ /** @brief Read a list of values from an ASCII file and create a cpl_bivector @param filename Name of the input ASCII file @return 1 newly allocated cpl_bivector or NULL on error @see cpl_vector_load The input ASCII file must contain two values per line. The returned object must be deallocated using cpl_bivector_delete() Two columns of numbers are expected in the input file. In addition to normal files, FIFO (see man mknod) are also supported. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_FILE_IO if the file cannot be read - CPL_ERROR_BAD_FILE_FORMAT if the file contains no valid lines */ /*----------------------------------------------------------------------------*/ cpl_bivector * cpl_bivector_read(const char * filename) { FILE * in; cpl_vector * v1; cpl_vector * v2; int np = 0; int size = 1000; /* Default size */ char line[1024]; double x, y; cpl_ensure(filename, CPL_ERROR_NULL_INPUT, NULL); /* Open the file */ in = fopen(filename, "r"); cpl_ensure(in, CPL_ERROR_FILE_IO, NULL); /* Create and fill the vectors */ v1 = cpl_vector_new(size); v2 = cpl_vector_new(size); while (fgets(line, 1024, in) != NULL) { if (line[0] != '#' && sscanf(line, "%lg %lg", &x, &y) == 2) { /* Found new element-pair - increase vector sizes if necessary, - insert element at end and - increment size counter */ if (np == size) { size *= 2; cpl_vector_set_size(v1, size); cpl_vector_set_size(v2, size); } cpl_vector_set(v1, np, x); cpl_vector_set(v2, np, y); np++; } } /* Check that the loop ended due to eof and not an error */ if (ferror(in)) { fclose(in); cpl_vector_delete(v1); cpl_vector_delete(v2); cpl_ensure(0, CPL_ERROR_FILE_IO, NULL); } fclose(in); /* Check that the file was not empty and set the size to its true value */ if (np == 0 || cpl_vector_set_size(v1, np) || cpl_vector_set_size(v2, np)) { cpl_vector_delete(v1); cpl_vector_delete(v2); cpl_ensure(0, CPL_ERROR_BAD_FILE_FORMAT, NULL); } /* Cannot fail at this point */ return cpl_bivector_wrap_vectors(v1, v2); } /*----------------------------------------------------------------------------*/ /** @brief Get the size of the cpl_bivector @param in the input bivector @return The size or a negative number on error Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if the input bivector contains vectors of different sizes */ /*----------------------------------------------------------------------------*/ int cpl_bivector_get_size(const cpl_bivector * in) { int size; cpl_ensure(in, CPL_ERROR_NULL_INPUT, -1); assert( in->x ); assert( in->y ); /* Get sizes of x and y vectors */ size = cpl_vector_get_size(in->x); cpl_ensure(size == cpl_vector_get_size(in->y), CPL_ERROR_ILLEGAL_INPUT, -2); return size; } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the x vector of the cpl_bivector @param in a cpl_bivector @return Pointer to the x vector or NULL on error The returned pointer refers to an already created cpl_vector. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_bivector_get_x(cpl_bivector * in) { cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL); assert( in->x ); return in->x; } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the x vector of the cpl_bivector @param in a cpl_bivector @return Pointer to the x vector or NULL on error @see cpl_bivector_get_x */ /*----------------------------------------------------------------------------*/ const cpl_vector * cpl_bivector_get_x_const (const cpl_bivector * in) { return cpl_bivector_get_x((cpl_bivector*)in); } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the y vector of the cpl_bivector @param in a cpl_bivector @return Pointer to the y vector or NULL on error The returned pointer refers to an already created cpl_vector. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_bivector_get_y(cpl_bivector * in) { cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL); assert( in->y ); return in->y; } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the y vector of the cpl_bivector @param in a cpl_bivector @return Pointer to the y vector or NULL on error */ /*----------------------------------------------------------------------------*/ const cpl_vector * cpl_bivector_get_y_const(const cpl_bivector * in) { return cpl_bivector_get_y((cpl_bivector*)in); } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the x data part of the cpl_bivector @param in a cpl_bivector @return Pointer to the double x array or NULL on error @see cpl_vector_get_data The returned pointer refers to already allocated data. @note Use at your own risk: direct manipulation of vector data rules out any check performed by the vector object interface, and may introduce inconsistencies between the information maintained internally, and the actual vector data and structure. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ double * cpl_bivector_get_x_data(cpl_bivector * in) { cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL); assert( in->x ); return cpl_vector_get_data(in->x); } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the x data part of the cpl_bivector @param in a cpl_bivector @return Pointer to the double x array or NULL on error @see cpl_bivector_get_x_data */ /*----------------------------------------------------------------------------*/ const double * cpl_bivector_get_x_data_const(const cpl_bivector * in) { return cpl_bivector_get_x_data((cpl_bivector*)in); } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the y data part of the cpl_bivector @param in a cpl_bivector @return Pointer to the double y array or NULL on error @see cpl_vector_get_x_data Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL */ /*----------------------------------------------------------------------------*/ double * cpl_bivector_get_y_data(cpl_bivector * in) { cpl_ensure(in, CPL_ERROR_NULL_INPUT, NULL); assert( in->y ); return cpl_vector_get_data(in->y); } /*----------------------------------------------------------------------------*/ /** @brief Get a pointer to the y data part of the cpl_bivector @param in a cpl_bivector @return Pointer to the double y array or NULL on error @see cpl_bivector_get_y_data */ /*----------------------------------------------------------------------------*/ const double * cpl_bivector_get_y_data_const(const cpl_bivector * in) { return cpl_bivector_get_y_data((cpl_bivector*)in); } /*----------------------------------------------------------------------------*/ /** @brief Linear interpolation of a 1d-function @param fref Reference 1d-function @param fout Preallocated 1d-function to contain result @return CPL_ERROR_NONE or the relevant #_cpl_error_code_ fref must have both its abscissa and ordinate defined. fout must have its abscissa defined and its ordinate allocated. The linear interpolation will be done from the values in fref to the abscissa points in fout. The abscissa points of both fref and fout must be growing, x_i < x_i+1. The abscissa points of fout must be in range of those of fref (i.e. extrapolation is not allowed). fref must be of at least length 2. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_DATA_NOT_FOUND if fout is too short - CPL_ERROR_ILLEGAL_INPUT if one of the requirements on the 2 input bivectors is not respected. */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_bivector_interpolate_linear(cpl_bivector * fout, const cpl_bivector * fref) { const cpl_error_code err = CPL_ERROR_ILLEGAL_INPUT; const int m = cpl_bivector_get_size(fref); const int n = cpl_bivector_get_size(fout); const double * xref; const double * yref; double * xout; double * yout; /* Initialize to avoid unjustified compiler warning */ double grad = 0; double y0 = 0; /* Start interpolation from below */ int iabove = 0; int ibelow = iabove - 1; /* Avoid (false) uninit warning */ int i; cpl_ensure_code(n > 0, cpl_error_get_code()); cpl_ensure_code(m > 0, cpl_error_get_code()); cpl_ensure_code(m > 1, CPL_ERROR_DATA_NOT_FOUND); xref = cpl_bivector_get_x_data_const(fref); yref = cpl_bivector_get_y_data_const(fref); xout = cpl_bivector_get_x_data(fout); yout = cpl_bivector_get_y_data(fout); assert( xref != NULL ); assert( yref != NULL); assert( xout != NULL); assert( yout != NULL); /* Verify that extrapolation is not necessary */ cpl_ensure_code(xref[0 ] <= xout[0 ], err); cpl_ensure_code(xout[0 ] < xout[n-1], err); cpl_ensure_code(xout[n-1] <= xref[m-1], err); for (i = 0; i < n; i++) { /* When possible reuse reference function abscissa points */ if (xout[i] > xref[iabove] || i == 0) { /* No, need new points */ while (xout[i] > xref[++iabove]); ibelow = iabove - 1; /* Verify that reference abscissa points are valid */ cpl_ensure_code(xref[iabove] > xref[ibelow], err); grad = (yref[iabove] - yref[ibelow]) / (xref[iabove] - xref[ibelow]); y0 = yref[ibelow] - grad * xref[ibelow]; } else /* Interpolation point may not be smaller than the lower reference point */ cpl_ensure_code(xout[i] >= xref[ibelow], err); yout[i] = y0 + grad * xout[i]; } return CPL_ERROR_NONE; } /**@}*/