File: detect.i

package info (click to toggle)
yorick-yutils 1.3.0-2
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 604 kB
  • ctags: 17
  • sloc: makefile: 105; python: 12
file content (482 lines) | stat: -rw-r--r-- 17,599 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
/*
 * detect.i --
 *	Detection of local minima/maxima for Yorick.
 *
 * Copyright (c) 2003, Eric THIEBAUT.
 *
 * 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 (to receive a  copy  of  the GNU
 * General Public License, write to the Free Software Foundation, Inc., 675
 * Mass Ave, Cambridge, MA 02139, USA).
 *
 * Routines:
 *      find_1d_minmax - find local minima/maxima in 1-D array.
 *      plot_1d_minmax - plot local minima/maxima in 1-D array.
 *	find_2d_max - find isolated local maxima in 2-D array.
 *
 * History:
 *	$Id: detect.i,v 1.1.1.1 2007/12/11 23:55:13 frigaut Exp $
 *	$Log: detect.i,v $
 *	Revision 1.1.1.1  2007/12/11 23:55:13  frigaut
 *	Initial Import - yorick-yutils
 *	
 *
 *-----------------------------------------------------------------------------
 */

func find_1d_minmax(a, what, inf=, sup=, alev=, rlev=, hysteresis=)
/* DOCUMENT find_1d_minmax(a)
 *     -or- find_1d_minmax(a, what)
 *   Find local minima/maxima in 1-D array  A.  If WHAT is nil or zero, the
 *   function returns an array of integers  with same shape as A and set to
 *   +1 where A is a local maximum, to -1 where A is a local minimum and to
 *   0  elsewhere.  Otherwise, the  function returns  the indices  of local
 *   maxima or local minima depending  whether WHAT is positive or negative
 *   (the result may be empty).  WHAT may also be a string: "any", "min" or
 *   "max".
 *
 *   Contiguous extrema,  say a  local maximum LOCMAX  and a  local minimum
 *   LOCMIN, are separated by a strict hysteresis (or a gap) such that:
 *
 *      LOCMAX - LOCMIN > HYSTERESIS*(max(A) - min(A)) >= 0
 *
 *   The default  is HYSTERESIS=0, i.e. all strict  local minima/maxima are
 *   detected.  However,  in order  to avoid being  too sensitive  to local
 *   extrema  (for  instance  because  of  noise), the  hysteresis  of  the
 *   algorithm  can  be  adjusted  by  keywords INF,  SUP,  ALEV,  RLEV  or
 *   HYSTERESIS.   Note  that  the   easiest  keywords  to  play  with  are
 *   HYSTERESIS or ATOL and RTOL (SUP and INF are more tricky to use).
 *
 *   HYSTERESIS = level of hysteresis  relative to peak-to-peak value of A.
 *         HYSTERESIS   must   be  conformable   with   A  and   everywhere
 *         non-negative (this  is not checked).  HYSTERESIS  is another way
 *         to  specify the absolute  tolerance and  is only  significant if
 *         ALEV is not specified.   Specifying HYSTERESIS gives an absolute
 *         tolerance:
 *            ALEV = (max(A) - min(A))*HYSTERESIS.
 *         Hence the lower  is the hysteresis, the more  local extrema will
 *         be  detected.   As  a  rule  of  thumb, a  good  value  for  the
 *         hysteresis is the 3-4 divided by the signal-to-noise-ratio.
 *
 *   ALEV = absolute level of  hysteresis.  ALEV must be conformable with A
 *         and such that  ALEV >= 0 everywhere (this  is not checked).  The
 *         default is  the same as with ALEV=0  (unless  keyword HYSTERESIS
 *         is specified).
 *
 *   RLEV = relative level of hysteresis.  RLEV must be conformable with A
 *         and such  that 0 <= RLEV  < 1 everywhere (this  is not checked).
 *         The default is the same as with RLEV=0.
 *
 *   INF = inferior  bound with respect to  a maximum: A(i) may  be a local
 *         maximum with respect to A(j) if  and only if A(j) < INF(i).  INF
 *         must  have the  same number  of elements  as A.   If INF  is not
 *         specified,  it is  computed from  the  value of  ALEV and  RLEV.
 *         Given ALEV and RLEV, the inferior bound is:
 *             INF = A - (ALEV + RLEV*abs(A))
 *
 *   SUP = superior  bound with respect to  a minimum: A(i) may  be a local
 *         minimum with respect to A(j) if  and only if A(j) > SUP(i).  SUP
 *         must  have the  same number  of elements  as A.   If SUP  is not
 *         specified,  it is  computed from  the  value of  ALEV and  RLEV.
 *         Given ALEV and RLEV, the superior bound is such that:
 *             A = SUP - (ALEV + RLEV*abs(SUP))
 *         or:
 *             SUP = (A + ALEV)/(1 - sign(SUP)*RLEV)
 *         since 0 <= RLEV < 1  then SUP has the same  sign as A + ALEV and
 *         finally:
 *             SUP = (A + ALEV)/(1 - sign(A + ALEV)*RLEV)
 *
 * SEE ALSO: plot_1d_minmax.
 */
{
  if (structof(what) == string) {
    if (what == "any") what = 0;
    else if (what == "max") what = 1;
    else if (what == "min") what = -1;
    else error, "bad value for WHAT (must be \"any\", \"min\" or \"max\")";
  }
  if (! is_array(a) || (dimsof(a)(1)) != 1) error, "expecting 1-D array";
  if ((s = structof(a)) != double) {
    if (s == complex) error, "illegal complex array";
    a = double(a);
  }
  n = numberof(a);

  /* compute the inferior/superior bounds */
  if (is_void(inf) || is_void(sup)) {
    if (is_void(alev)) {
      if (is_void(hysteresis)) alev = 0.0;
      else alev = (max(a) - min(a))*hysteresis;
    }
    if (noneof(rlev)) {
      if (noneof(alev)) {
        if (is_void(inf)) eq_nocopy, inf, a;
        if (is_void(sup)) eq_nocopy, sup, a;
      } else {
        if (is_void(inf)) inf = a - alev;
        if (is_void(sup)) sup = a + alev;
      }
    } else {
      if (is_void(inf)) inf = a - rlev*abs(a) - alev;
      if (is_void(sup)) {
        sup = a + alev; /* temporary */
        sup = sup/(1.0 - sign(sup)*rlev);
      }
    }
  }
  type = array(long, n);

  /* start with leftmost element */
  imin = imax = i = 1;
  vmin = vmax = a(1);
  vinf = inf(1);
  vsup = sup(1);
  state = 0;
  for (;;) {
    val = a(++i);
    if (i == n) {
      /* end-point */
      if (state > 0) {
        type((val > vmax ? i : imax)) = 1;
      } else if (state < 0) {
        type((val < vmin ? i : imin)) = -1;
      }
      break;
    }
    if (state >= 0) {
      /* seeking for a local maximum */
      if (val < vinf) {
        /* accept maximum value found so far as a local maximum and setup
           to start seeking for next local minimum */
        type(imax) = 1;
        imin = i;
        vmin = val;
        vsup = sup(i);
        state = -1;
      } else if (val > vmax) {
        /* higher maximum found */
        imax = i;
        vmax = val;
        vinf = inf(i);
      }
    }
    if (state <= 0) {
      /* seeking for a local minimum */
      if (val > vsup) {
        /* accept maximum value found so far as a local maximum and setup
           to start seeking for a local minimum */
        type(imin) = -1;
        imax = i;
        vmax = val;
        vinf = inf(i);
        state = 1;
      } else if (val < vmin) {
        /* lower minimum found */
        imin = i;
        vmin = val;
        vsup = sup(i);
      }
    }
  }
  return (what ? (what > 0 ? where(type > 0) : where(type > 0)) : type);
}

func plot_1d_minmax(y, x, list, nocurve=, type=, width=, color=,
                    symbol=, size=, fill=,
                    what=, inf=, sup=, alev=, rlev=, hysteresis=)
/* DOCUMENT plot_1d_minmax, y;
       -or- plot_1d_minmax, y, x;
       -or- plot_1d_minmax, y, x, list;

     Plot (X,Y) curve with local minima/maxima.  LIST is the list of
     extrema as returned by find_1d_minmax; if LIST is nil, find_1d_minmax
     is used to find the extrema (with argument Y, and values of keywords
     WHAT, INF, SUP, ALEV, RLEV and/or HYSTERESIS).

     Unless keyword NOCURVE is true, the curve (X,Y) is plotted as well
     (with values of keywords TYPE, WIDTH and/or COLOR).

     Keywords SYMBOL, SIZE, FILL, and COLOR can be used to customize the
     plotting of local minima/maxima.  If SYMBOL is unspecified and both
     minima and maxima are to be plotted, triangles pointing to the top (to
     the bottom) will be used to display maxima (minima).

     When called as a function, actual LIST is returned.

   SEE ALSO: find_1d_minmax, plp, plg. */
{
  if (is_void(list)) list = find_1d_minmax(y, what, inf=inf, sup=sup,
                                           alev=alev, rlev=rlev,
                                           hysteresis=hysteresis);
  if (is_void(x)) x = double(indgen(numberof(y)));
  if (! nocurve) plg, y, x, color=color, type=type, width=width;
  if (numberof(list) == numberof(y)) {
    /* plot both minima and maxima */
    if (is_array((i = where(list < 0)))) {
      plp, y(i), x(i), symbol=(is_void(symbol) ? 7 : symbol),
        size=size, color=color, fill=fill;
    }
    if (is_array((i = where(list > 0)))) {
      plp, y(i), x(i), symbol=(is_void(symbol) ? 3 : symbol),
        size=size, color=color, fill=fill;
    }
  } else if (! is_void(list)) {
    plp, y(list), x(list), symbol=symbol, size=size, color=color, fill=fill;
  }
  return list;
}

/*---------------------------------------------------------------------------*/

func find_2d_max(img, alev=, rlev=, cmin=, cmax=, bad=, debug=)
/* DOCUMENT map = find_2d_max(img)
     Find disjoint  local maxima in  2-D array IMG  and return an  array of
     integers MAP with same shape as IMG and set as follow:
         MAP(x,y) = -1       for bad pixels (or strictly above CMAX)
         MAP(x,y) =  0       for pixels not assigned to any local maximum
         MAP(x,y) =  N (N>0) for pixels assigned to N-th local maximum
     The  maxima  are  labelled  from   the  higher  to  the  lower.  Hence
     where(MAP==1) gives the indices of pixels around the stronger maximum.

     The selection  works as follows.   The algorithm starts with  the next
     (unassigned) higher maximum and  then marks all connected pixels which
     are greater or equal to a given  threshold.  If a bad pixel or a pixel
     already assigned to another  maximum is encountered during this stage,
     the algorithm  gives up this maximum  and proceeds with  the next one.
     Otherwise,  the marked  region  is  assigned to  the  maximum and  the
     algorithm proceeds with the  next one.  This algorithm guarantees that
     the marked regions are all separated (by at least one pixel) from each
     other and from any bad pixel.  The threshold reads:
         THRESHOLD = PEAK - RLEV*abs(PEAK) - ALEV
     where PEAK is the current  maximum, ALEV (ALEV>=0 everywhere) and RLEV
     (0<=RLEV<1 everywhere) are the absolute and relative threshold levels.
     ALEV and RLEV are given by keyword. By default, ALEV=0 and RLEV=0.
     
     Keyword CMIN can be used to  set the minimum value of a local maximum.
     Since searching all  maxima may be prohibitively long,  it is strongly
     recommended to limit the depth of the search by the keyword CMIN.

     Keywords BAD  and/or CMAX can be used  to mark as bad  pixels the ones
     for which BAD is non-zero and/or which are (strictly) above CMAX.

     If  specified, keywords ALEV,  RLEV, CMIN,  CMAX and  BAD must  all be
     conformable with  IMG; you can  therefore setup things on  a per-pixel
     basis, or columnwise, or rowwise...

     
  EXAMPLE:  
     For instance, if SIGMA is the standard deviation of noise in the image
     and  BACKGROUND is  its background  level (both  could  be pixelwise),
     then:
         find_2d_max(IMG, cmin=BACKGROUND+3*SIGMA, alev=4*SIGMA)
     will find all the maxima in  IMG which are above the background with a
     3 SIGMA  confidence level  and mark the  regions around  every maximum
     (with value PEAK) where connected pixels are such that:
         IMG >= PEAK - 4*SIGMA

   HINTS:
     1. Use  keyword CMIN  to limit  the  search (possibly  on a  per-pixel
        basis).

     2. The search  necessitates to sort  the pixels elligible to  be local
        maxima (all which are above CMIN,  below CMAX and not in BAD), this
        sorting can be very long for large images (again use CMIN) but also
        for integer valued images  (a drawback of the quicksort algorithm?)
        to overcome this  it is sufficient to add a  small amount of random
        noise in the image, for instance:
            find_2d_max(IMG + 1e-9*(random(dimsof(IMG)) - 0.5), ...)
        but beware that this can make the result (slightly) inpredictible.
      
   SEE ALSO: sort, find_1d_minmax. */
{
  if (! is_array(img) ||
      ((dims = dimsof(img))(1)) != 2) error, "expecting 2-D array";
#if 0
  if ((s = structof(img)) != double) {
    if (s == complex) error, "illegal complex array";
    img = double(img);
  }
#endif

  /* mark bad points */
  region = array(long, dims); /* needed to make BAD conformable with IMG */
  if (is_void(bad)) {
    if (! is_void(cmax)) bad = (img > cmax);
  } else if (is_void(cmax)) {
    bad |= region; /* make BAD conformable with IMG */
  } else {
    bad |= (img > cmax);
  }
  if (anyof(bad)) region(where(bad)) = -1;

  /* sort pixels eligible for being local maxima */
  if (debug) write, format="%s...", "sorting";
  if (is_void(cmin) && is_void(bad)) {
    index = sort(img(*));
  } else {
    if (is_void(cmin))     index = where(! bad);
    else if (is_void(bad)) index = where(img >= cmin);
    else                   index = where((! bad) & (img >= cmin));
    if (! is_array(index)) {
      write, "warning no pixel are eligible for being local maxima";
      return region;
    }
    index = index(sort(img(index)));
  }
  if (debug) write, format="%s\n", "done";
  bad = cmin = cmax = []; /* free some memory */

  /* compute threshold */
  if (is_void(rlev)) {
    if (is_void(alev)) threshold = img(index);
    else               threshold = (img - alev)(index);
  } else {
    if (is_void(alev)) threshold = (img - rlev*abs(img))(index);
    else               threshold = (img - alev - rlev*abs(img))(index);
  }
  alev = rlev = []; /* free some memory */

  /* serach local maxima */
  number = numberof(img);
  width = dims(2);
  height = dims(3);
  list = array(long, number); /* indices of pixels in current region */
  state = array(long, dims);
  mark = 1;
  for (i=numberof(index) ; i>=1 ; --i) {
    j = index(i);
    if (region(j)) continue;
    level = threshold(i);
#if 0
    if (debug) {
      write, format="search max around (%d,%d) %g >= %g\n",
        1 + (j - 1)%width, 1 + (j - 1)/width,double(img(j)),
        double(level);
    }
#endif
    
    /* Use a kind of non-recursive flood-fill algorithm.
     *
     * The 3 following bits are used to indicate the directions to
     * investigate (so that we limit the number of checks undergone by a
     * pixel):
     *
     *     +---+
     *     | 4 |
     * +---+---+---+
     * | 2 | x | 1 |
     * +---+---+---+
     *     | 8 |
     *     +---+
     *
     * region(x,y) =  0  if unused
     * region(x,y) = -1  if invalid
     * region(x,y) =  n  if inside n-th blob
     *
     * TO DO: use same array for REGION and STATE
     *
     */
    region(j) = mark; /* mark current maximum */
    state(j) = 15;    /* will check all neighbors of current maximum */
    count = 1;        /* number of pixels in the current region */
    list(1) = j;      /* current maximum belongs to current region */
    discard = 0n;     /* no error yet */
    for (j=1 ; j<=count ; ++j) {
      k = list(j);
      x = 1 + (k - 1)%width;
      y = 1 + (k - 1)/width;
      s = state(k);
      //if (debug) write,format="state(%d,%d)=\n";
      if (s & 1) {
        if (x < width) {
          l = k + 1;
          if (! (r = region(l))) {
            if (img(l) >= level) {
              region(l) = mark;
              state(l)  = 13; /* (1|2|4|8) & ~2 */
              list(++count) = l;
            }
          } else if (r == mark) {
            state(l) = (s & 13);
          } else {
            discard = 1n;
            break;
          }
        }
      }
      if (s & 2) {
        if (x > 1) {
          l = k - 1;
          if (! (r = region(l))) {
            if (img(l) >= level) {
              region(l) = mark;
              state(l)  = 14; /* (1|2|4|8) & ~1 */
              list(++count) = l;
            }
          } else if (r == mark) {
            state(l) = (s & 14);
          } else {
            discard = 1n;
            break;
          }
        }
      }
      if (s & 4) {
        if (y < height) {
          l = k + width;
          if (! (r = region(l))) {
            if (img(l) >= level) {
              region(l) = mark;
              state(l)  = 7; /* (1|2|4|8) & ~8 */
              list(++count) = l;
            }
          } else if (r == mark) {
            state(l) = (s & 7);
          } else {
            discard = 1n;
            break;
          }
        }
      }
      if (s & 8) {
        if (y > 1) {
          l = k - width;
          if (! (r = region(l))) {
            if (img(l) >= level) {
              region(l) = mark;
              state(l)  = 11; /* (1|2|4|8) & ~4 */
              list(++count) = l;
            }
          } else if (r == mark) {
            state(l) = (s & 11);
          } else {
            discard = 1n;
            break;
          }
        }
      }
    }
    l = list(1:count);
    if (discard) {
      region(l) = 0;
    } else {
      ++mark;
      if (debug) {
        fma;
        //pli, region > 0;
        pli, region;
        pause, 1;
      }
    }
  }
  return region;
}

/*---------------------------------------------------------------------------*/