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 /*--------------------------------------------------------------------------*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001