detect_peaks.h

00001 
00002 /*----------------------------------------------------------------------------
00003  *                                  E.S.O.
00004  *----------------------------------------------------------------------------
00005  * File name    :   detect_peaks.h
00006  * Author       :   Nicolas Devillard
00007  * Created on   :   March 17, 1997
00008  * Hardware     :   Sun Sparc 20
00009  * Software     :   ANSI C under Solaris Unix
00010  *                  Part of ECLIPSE library for Adonis
00011  * Description  :   peak detection routines
00012  *--------------------------------------------------------------------------*/
00013 
00014 /*
00015 
00016  $Id: detect_peaks.h,v 1.1 2003/09/03 12:50:47 amodigli Exp $
00017  $Author: amodigli $
00018  $Date: 2003/09/03 12:50:47 $
00019  $Revision: 1.1 $
00020 
00021  */
00022 
00023 #ifndef _DETECT_PEAKS_H_
00024 #define _DETECT_PEAKS_H_
00025 
00026 
00027 /*----------------------------------------------------------------------------
00028  *                              Includes
00029  *--------------------------------------------------------------------------*/
00030 
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 
00034 #include "memory.h"
00035 #include "pixelmaps.h"
00036 #include "binary_filters.h"
00037 
00038 #define DEFAULT_SIGMA_THRESHOLD     2.0
00039 
00040 /*
00041  * Printing flags
00042  * to set a flag, do for example:
00043  *
00044  *      int pflag = 0 ;
00045  *      pflag |= (1<<PPOS_X) ;
00046  *      pflag |= (1<<PPOS_Y) ;
00047  *
00048  * to read a flag, do for example:
00049  *
00050  *      boolean output_x, output_y ;
00051  *      output_x = (boolean)(pflag & (1<<PPOS_X)) ;
00052  *      output_y = (boolean)(pflag & (1<<PPOS_Y)) ;
00053  *      if (output_x == TRUE) { ... }
00054  */
00055 
00056 #define PPOS_NO_OUTPUT              0
00057 #define PPOS_SAOIMAGE               1
00058 #define PPOS_X                      2
00059 #define PPOS_Y                      3
00060 #define PPOS_MAG                    4
00061 
00062 
00063 /*----------------------------------------------------------------------------
00064  *                              New Types
00065  *--------------------------------------------------------------------------*/
00066 
00067 /*
00068  * The following structure contains all information necessary for
00069  * recursive search of consistent zones of pixels.  see the functions
00070  * below to have an idea of how it is used. floodfill_from_pixel() is
00071  * recursive get_pixelrank_if_valid() is the checking function is for
00072  * recursive evals.
00073  */
00074 
00075 struct pixel_accumulator {
00076     /*
00077      * This contains accumulated pixel coordinates, in X and Y
00078      */
00079     double   accXpos,
00080              accYpos ;
00081 
00082     /*
00083      * This stores the accumulated pixelvalues
00084      */
00085 
00086     pixelvalue  accIntensity ;
00087             
00088     /*
00089      * PixelsInCurrentBlock
00090      * Contains the number of pixels accumulated so far in the current
00091      * consistent block of pixels.
00092      */
00093     int     PixelsInCurrentBlock ;
00094 
00095     /*
00096      * Contains the image width, necessary to compute positions in the
00097      * image
00098      * ... the image height is not requested, because no test is made to
00099      *     check positions outside the image plane.
00100      */
00101     int     imWidth;
00102 
00103     /*
00104      * contains the rank of the current pixel being examined
00105      * mostly used to optimize speed.
00106      */
00107     int     CurrentRank ;
00108 
00109     /*
00110      * Contains the total number of pixels set to 1 in the currently
00111      * examined pixel_map
00112      */
00113     int     Total1Pixels ;
00114 
00115     /*
00116      * Contains the number of pixels which still have to be examined
00117      */
00118     int     RemainingPixels ;
00119 
00120     /*
00121      * This table contains a flag per pixel set to 1 in the pixel_map.
00122      * a flag set to TRUE means the pixel has not been associated to
00123      * a zone which has been found. If it is set to FALSE, the pixel
00124      * has already been associated with a white zone.
00125      */
00126     boolean *ValidPixels ;
00127 
00128     /*
00129      * This is the most important stuff: this array contains all
00130      * the positions of the white pixels in the input image.
00131      * Positions are not referenced by a couple (x,y) but by
00132      * a single number POS = x + ImWidth * y
00133      * Knowing Imwidth, it is possible to retrieve x and y values,
00134      * but in many parts of the algorithm, the single pixel position
00135      * in the image array is enough.
00136      */
00137     ulong32 *pos ;
00138 
00139     /*
00140      * This is the attached input image
00141      * it should be just a pointer copy
00142      */
00143     OneImage * refImage ;
00144 } ;
00145 
00146 
00147  
00148 /*----------------------------------------------------------------------------
00149  *                      Function ANSI C prototypes
00150  *--------------------------------------------------------------------------*/
00151 
00152 
00153 /*----------------------------------------------------------------------------
00154  * Function :   find_bright_objects()
00155  * In       :   an image, a sigma rejection threshold (floa), and
00156  *              intermediate printing flags.
00157  *              To have the function return a list of pixel positions,
00158  *              set print_flags to 0. n_objects is then not updated and
00159  *              NULL is returned.
00160  *              Other flags are described above.
00161  * Out      :   a list of pixel positions, and a number of pixels
00162  *              or: a list of positions in requested format on stdout.
00163  * Job      :   find out the peaks in an astronomical image
00164  * Notice   :
00165  * Algorithm:
00166  *
00167  *  The idea is first to threshold the image to get all pixels
00168  * which are (quite a lot) above the average image value. This threshold
00169  * is made at average + k * sigma , where:
00170  *
00171  *      average is the average image pixel value
00172  *      sigma is the standard deviation
00173  *      k is a user-defined parameter (actually, the 'sigma' parameter
00174  *      in this function), it is usually defaulted to 2.0
00175  *      Note that for noisy signals, this estimation is biased and it is
00176  *      not possible to use mean and stdev directly. To increase
00177  *      robustness, we use the plane's median value as an estimate of
00178  *      the mean, and the average absolute distance to this median as an
00179  *      estimate of the stdev. This proves quite reliable.
00180  *
00181  * This first operation yields a binary map of all high pixel-valued
00182  * positions in the image. To avoid detecting bad pixels which have an
00183  * extremely high value but are located to a single pixel, a morphological
00184  * closing is applied to the image, which removes all regions smaller than
00185  * the structuring element, in this case a 3x3 square.
00186  *
00187  * NB:
00188  * This also removes very small objects, and is the major cause for
00189  * non-detections of this algorithm.
00190  *
00191  * Then a FloodFill algorithm is applied to gather pixel positions
00192  * in groups belonging to the same white blob. This algorithm is
00193  * recursively called, information is stored in a pixel_accumulator
00194  * structure for recursivity.
00195  * For each found white blob, the algorithm returns:
00196  *      the average x position
00197  *      the average y position
00198  *
00199  * Possible improvements to this algorithm:
00200  *
00201  *      The threshold criterion (avg + k* sigma) is certainly not
00202  *      the best. It could be possible to determine a better one
00203  *      from the histogram, for example.
00204  *
00205  *      By changing the structuring element in the morphological closing,
00206  *      it should be possible to detect smaller objects without detecting
00207  *      bad pixels.
00208  *
00209  *      N. Devillard, Wed Mar 19th, 1997
00210  *
00211  *--------------------------------------------------------------------------*/
00212 
00213 pixel_position *
00214 find_bright_objects(
00215     OneImage    *   image1, 
00216     double          sigma,
00217     int         *   n_objects,
00218     int             print_flags
00219 ) ;
00220 
00221 
00222 /*----------------------------------------------------------------------------
00223  * Function :   find_centers_on_pixelmap()
00224  * In       :   binary pixel map
00225  *              flag: printIntermediate set to TRUE if results are printed
00226  *              as they are computed
00227  * Out      :   list of pixel positions, number of objects
00228  * Job      :   find out the positions for the centers of white blobs in
00229  *              a binary pixel map.
00230  * Notice   :   if the image contains numerous white blobs, may be long.
00231  *              print_flags are described below in fprint_1pix()
00232  *--------------------------------------------------------------------------*/
00233 
00234 pixel_position *
00235 find_centers_on_pixelmap(
00236     OneImage    *   refImage,
00237     pixel_map   *   map,
00238     int         *   n_objects,
00239     int             print_flags
00240 ) ;
00241 
00242 
00243 /*
00244  * --------- Pixel position output to file:
00245  * the following functions are divided into a general function to print
00246  * a list of pixel positions and a function to print a single position.
00247  * This allows in 'detpeak' to output pixel positions at the time they
00248  * are found, and allows other functions to print a list of pixels once
00249  * they are established.
00250  */
00251 
00252 
00253 /*----------------------------------------------------------------------------
00254  * Function :   fprint_pixel_positions()
00255  * In       :   FILE pointer to print out to, list of pixel positions,
00256  *              number of pixels
00257  * Out      :   formatted list of pixel positions in requested file
00258  * Job      :   print out to txt file a list of pixel positions
00259  * Notice   :   it is up to the caller to provide an allocated and valid
00260  *              file pointer and to close it afterwards.
00261  *--------------------------------------------------------------------------*/
00262 
00263 void
00264 fprint_pixel_positions(
00265     FILE            *   f_out,
00266     pixel_position  *   p_list,
00267     int                 npix,
00268     int                 pflags
00269 ) ;
00270 
00271 
00272 /*----------------------------------------------------------------------------
00273  * Function :   fprint_1pix()
00274  * In       :   1 pixel_position, print flags, a file pointer
00275  * Out      :   output to file pointer
00276  * Job      :   prints out to the given file pointer some information
00277  *              about the provided pixel position
00278  * Notice   :   printing flags definition: 
00279  *
00280  *              PPOS_NO_OUTPUT      No output [*]
00281  *
00282  *              PPOS_SAOIMAGE       SAOimage type [*]
00283  *
00284  *              PPOS_X              Output X position
00285  *              PPOS_Y              Output Y position
00286  *              PPOS_MAG            Output magnitude
00287  *
00288  *              [*] means: excludes all other flags
00289  *--------------------------------------------------------------------------*/
00290 
00291 void fprint_1pix(
00292     FILE            *   f_out,
00293     pixel_position  *   pix1,
00294     int                 pflags
00295 ) ;
00296 
00297 
00298 
00299 /*----------------------------------------------------------------------------
00300  * Function :   fine_position_centers()
00301  * In       :   an image, a list of pixel positions, a number of pixels
00302  *              corresponding to that list, and 3 doubles
00303  * Out      :   void
00304  * Job      :   refine the peak centers, by computing a barycenter
00305  * Notice   :   the 3 doubles (r1, r2, r3) define 3 circle radiuses:
00306  *
00307  *              r2 and r3 define a ring centered on each pixel in the list,
00308  *              around which the background will be estimated.
00309  *
00310  *              r1 defines a disk in which the barycenter will be computed,
00311  *              i.e. the centroid pixel, weighted by pixel values, from
00312  *              which the background value has been subtracted:
00313  *
00314  *              x_center = (1/sum(pixel[i]-bg))*(sum(x[i] * (pixel[i]-bg)))
00315  *              y_center = (1/sum(pixel[i]-bg))*(sum(y[i] * (pixel[i]-bg)))
00316  *
00317  *--------------------------------------------------------------------------*/
00318 
00319 void
00320 fine_position_centers(
00321     OneImage        *   in,
00322     pixel_position  *   p_list,
00323     int                 npix,
00324     double              r1,
00325     double              r2,
00326     double              r3
00327 ) ;
00328 
00329 
00330 /*----------------------------------------------------------------------------
00331  * Function :   sort_ppos_by_decmag()
00332  * In       :   2 pixel positions
00333  * Out      :   integer: 1 or -1
00334  * Job      :   compares two pixel positions
00335  * Notice   :   to be used with qsort() (from libc)
00336  *              to sort out a list of pixel positions by decreasing
00337  *              magnitude.
00338  *--------------------------------------------------------------------------*/
00339 
00340 int sort_ppos_by_decmag(
00341     const void * p1,
00342     const void * p2
00343 ) ;
00344 
00345 
00346 /*----------------------------------------------------------------------------
00347  * Function :   sort_ppos_by_incmag()
00348  * In       :   2 pixel positions
00349  * Out      :   integer: 1 or -1
00350  * Job      :   compares two pixel positions
00351  * Notice   :   to be used with qsort() (from libc)
00352  *              to sort out a list of pixel positions by increasing
00353  *              magnitude.
00354  *--------------------------------------------------------------------------*/
00355 
00356 int sort_ppos_by_incmag(
00357     const void * p1,
00358     const void * p2
00359 ) ;
00360 
00361 /*---------------------------------------------------------------------------
00362    Function :   sort_out_ppos_by_radius()
00363    In       :   1 list of pixel_position, number of positions in the list,
00364                 pivot position in 2d.
00365    Out      :   void
00366    Job      :   sort out a given list according to a criterion of
00367                 increasing distance to a given point at (cx, cy).
00368    Notice   :
00369  ---------------------------------------------------------------------------*/
00370 
00371 
00372 void sort_out_ppos_by_radius(
00373     pixel_position  *   peaks,
00374     int                 npeaks,
00375     double              cx,
00376     double              cy
00377 ) ;   
00378 
00379 
00380 
00381 /*---------------------------------------------------------------------------
00382    Function :   locate_peak_in_window()
00383    In       :   1 image, position to look for a peak, window size
00384    Out      :   refined position (pixel accuracy) will be written in
00385                 input offset table (2 ints)
00386                 returns 0 if Ok, -1 otherwise
00387    Job      :   precise the location of a local maximum in a window
00388    Notice   :
00389    Algorithm:
00390 
00391                 * Extract the sub-window as a vignette
00392                 * Filter it out with a median filter
00393                 * Locate the maximum pixel in this filtered version
00394                 * Return coordinates of this local maximum
00395 
00396  ---------------------------------------------------------------------------*/
00397 
00398 int
00399 locate_peak_in_window(
00400     OneImage    *   img,
00401     int             px,
00402     int             py,
00403     int             search_hx,
00404     int             search_hy,
00405     int         *   refpos
00406 );
00407 
00408 #endif
00409 /*--------------------------------------------------------------------------*/

Generated on Wed Oct 26 13:08:52 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001