File: find_boundary.pro

package info (click to toggle)
coyote 2019.01.29-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster
  • size: 6,316 kB
  • sloc: python: 184; makefile: 14; sh: 13
file content (327 lines) | stat: -rw-r--r-- 13,477 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
;+
; NAME:
;       FIND_BOUNDARY
;
; PURPOSE:
;
;       This program finds the boundary points about a region of interest (ROI)
;       represented by pixel indices. It uses a "chain-code" algorithm for finding
;       the boundary pixels.
;
; AUTHOR:
;
;       FANNING SOFTWARE CONSULTING
;       David Fanning, Ph.D.
;       1645 Sheely Drive
;       Fort Collins, CO 80526 USA
;       Phone: 970-221-0438
;       E-mail: david@idlcoyote.com
;       Coyote's Guide to IDL Programming: http://www.idlcoyote.com
;
; CATEGORY:
;
;       Graphics, math.
;
; CALLING SEQUENCE:
;
;       boundaryPts = Find_Boundary(indices, XSize=xsize, YSize=ysize)
;
; OPTIONAL INPUTS:
;
;       indices - A 1D vector of pixel indices that describe the ROI. For example,
;            the indices may be returned as a result of the WHERE function.
;
; OUTPUTS:
;
;       boundaryPts - A 2-by-n points array of the X and Y points that describe the
;            boundary. The points are scaled if the SCALE keyword is used.
;
; INPUT KEYWORDS:
;
;       SCALE - A one-element or two-element array of the pixel scale factors, [xscale, yscale],
;            used to calculate the perimeter length or area of the ROI. The SCALE keyword is
;            NOT applied to the boundary points. By default, SCALE=[1,1].
;
;       XSIZE - The X size of the window or array from which the ROI indices are taken.
;            Set to !D.X_Size by default.
;
;       YSIZE - The Y size of the window or array from which the ROI indices are taken.
;            Set to !D.Y_Size by default.
;
; OUTPUT KEYWORDS:
;
;       AREA - A named variable that contains the pixel area represented by the input pixel indices,
;            scaled by the SCALE factors.
;
;       CENTER - A named variable that contains a two-element array containing the center point or
;            centroid of the ROI. The centroid is the position in the ROI that the ROI would
;            balance on if all the index pixels were equally weighted. The output is a two-element
;            floating-point array in device coordinate system, unless the SCALE keyword is used,
;            in which case the values will be in the scaled coordinate system.
;
;       PERIM_AREA - A named variable that contains the (scaled) area represented by the perimeter
;            points, as indicated by John Russ in _The Image Processing Handbook, 2nd Edition_ on
;            page 490. This is the same "perimeter" that is returned by IDLanROI in its
;            ComputeGeometry method, for example. In general, the perimeter area will be
;            smaller than the pixel area.
;
;       PERIMETER - A named variable that will contain the perimeter length of the boundary
;            upon returning from the function, scaled by the SCALE factors.
;
;  EXAMPLE:
;
;       LoadCT, 0, /Silent
;       image = BytArr(400, 300)+125
;       image[125:175, 180:245] = 255B
;       indices = Where(image EQ 255)
;       Window, XSize=400, YSize=300
;       TV, image
;       PLOTS, Find_Boundary(indices, XSize=400, YSize=300, Perimeter=length), $
;           /Device, Color=cgColor('red')
;       Print, length
;           230.0
;
; MODIFICATION HISTORY:
;
;       Written by David W. Fanning, April 2002. Based on an algorithm written by Guy
;       Blanchard and provided by Richard Adams.
;       Fixed a problem with distinction between solitary points and
;          isolated points (a single point connected on a diagonal to
;          the rest of the mask) in which the program can't get back to
;          the starting pixel. 2 Nov 2002. DWF
;       Added the ability to return the perimeter length with PERIMETER and
;           SCALE keywords. 2 Nov 2002. DWF.
;       Added AREA keyword to return area enclosed by boundary. 2 Nov 2002. DWF.
;       Fixed a problem with POLYFILLV under-reporting the area by removing
;           POLYFILLV and using a pixel counting method. 10 Dec 2002. DWF.
;       Added the PERIM_AREA and CENTER keywords. 15 December 2002. DWF.
;       Replaced the cgErrorMsg routine with the latest version. 15 December 2002. DWF.
;       Fixed a problem in which XSIZE and YSIZE have to be specified as integers to work. 6 March 2006. DWF.
;       Fixed a small problem with very small ROIs that caused the program to crash. 1 October 2008. DWF.
;       Modified the algorithm that determines the number of boundary points for small ROIs. 28 Sept 2010. DWF.
;       Modified the algorithm that determines if we are "home" to correct a problem discovered by Peter Vogt in
;          which small regions connect to larger regions by pixel corners only was cut off. 12 Mar 2014. Peter Vogt.
;-
;******************************************************************************************;
;  Copyright (c) 2008-2014, by Fanning Software Consulting, Inc.                           ;
;  All rights reserved.                                                                    ;
;                                                                                          ;
;  Redistribution and use in source and binary forms, with or without                      ;
;  modification, are permitted provided that the following conditions are met:             ;
;                                                                                          ;
;      * Redistributions of source code must retain the above copyright                    ;
;        notice, this list of conditions and the following disclaimer.                     ;
;      * Redistributions in binary form must reproduce the above copyright                 ;
;        notice, this list of conditions and the following disclaimer in the               ;
;        documentation and/or other materials provided with the distribution.              ;
;      * Neither the name of Fanning Software Consulting, Inc. nor the names of its        ;
;        contributors may be used to endorse or promote products derived from this         ;
;        software without specific prior written permission.                               ;
;                                                                                          ;
;  THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY        ;
;  EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES    ;
;  OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT     ;
;  SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT,             ;
;  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED    ;
;  TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;         ;
;  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND             ;
;  ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT              ;
;  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS           ;
;  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                            ;
;******************************************************************************************;
FUNCTION Find_Boundary_Outline, mask, darray, boundaryPts, ptIndex, $
   xsize, ysize, from_direction

On_Error, 2

FOR j=1,7 DO BEGIN

   to_direction = (from_direction + j) MOD 8
   newPt = boundaryPts[*,ptIndex-1] + darray[*,to_direction]

      ; If this is the edge, assume it is a background point.

   IF (newpt[0] LT 0 OR newpt[0] GE xsize OR newpt[1] LT 0 OR $
       newpt[1] GE ysize) THEN CONTINUE
   IF mask[newPt[0], newPt[1]] NE 0 THEN BEGIN
      boundaryPts[*,ptIndex] = newPt

     ; Return the "from" direction.
      RETURN, (to_direction + 4) MOD 8
   ENDIF

ENDFOR

   ; If we get this far, this is either a solitary point or an isolated point.

IF TOTAL(mask GT 0) GT 1 THEN BEGIN ; Isolated point.
   newPt = boundaryPts[*,ptIndex-1] + darray[*,from_direction]
   boundaryPts[*,ptIndex] = newPt
   RETURN, (from_direction + 4) MOD 8
ENDIF ELSE BEGIN ; Solitary point.
   boundaryPts[*,ptIndex] = boundaryPts[*,ptIndex-1]
   RETURN, -1
ENDELSE
END ; ------------------------------------------------------------------------------------------



FUNCTION Find_Boundary, indices, $
   AREA=area, $
   CENTER=center, $
   PERIM_AREA=perim_area, $
   PERIMETER=perimeter, $
   SCALE=scale, $
   XSIZE=xsize, $
   YSIZE=ysize


Catch, theError
IF theError NE 0 THEN BEGIN
   Catch, /Cancel
   ok = cgErrorMsg()
   RETURN, -1
ENDIF


IF N_Elements(indices) EQ 0 THEN Message, 'Indices of boundary region are required. Returning...'
IF N_Elements(scale) EQ 0 THEN BEGIN
   diagonal = SQRT(2.0D)
ENDIF ELSE BEGIN
   scale = Double(scale)
   diagonal = SQRT(scale[0]^2 + scale[1]^2)
ENDELSE
IF N_Elements(xsize) EQ 0 THEN xsize = !D.X_Size ELSE xsize = Long(xsize)
IF N_Elements(ysize) EQ 0 THEN ysize = !D.Y_Size ELSE ysize = Long(ysize)
IF Arg_Present(perimeter) THEN perimeter = 0.0

   ; Create a mask with boundary region inside.

indices = indices[Uniq(indices)]
mask = BytArr(xsize, ysize)
mask[indices] = 255B

   ; Set up a direction array.

darray = [[1,0],[1,1],[0,1],[-1,1],[-1,0],[-1,-1],[0,-1],[1,-1]]

   ; Find a starting point. The pixel to the left of
   ; this point is background

i = Where(mask GT 0)
firstPt = [i[0] MOD xsize, i[0] / xsize]
from_direction = 4

   ; Set up output points array. For narrow ROIs, we need to construct
   ; a different sort of algoritm for the number of boundary points.
IF (xsize LT 4) OR (ysize LT 4) THEN BEGIN
    boundaryPts = LonArr(2, (2*Max([xsize,ysize]) + 2*Min([xsize,ysize])))
ENDIF ELSE BEGIN
    boundaryPts = LonArr(2, (Long(xsize) * ysize / 4L) + 1)
ENDELSE
boundaryPts[0] = firstPt
ptIndex = 0L
secPt = LonArr(2)

   ;   We shall not cease from exploration
   ;   And the end of all our exploring
   ;   Will be to arrive where we started
   ;   And know the place for the first time
   ;
   ;                     T.S. Eliot
REPEAT BEGIN

   ; The break condition is normally to get back to the first point, but this is not enough to
   ; always find the perimeter. In particualr, small areas that are connected only by corner
   ; pixels touching can get missed or cut off by using this algorithm. To get around this, we also
   ; check the second pixel. If the first and second pixel are the same, then we know we have
   ; come around to the same spot and are retracing our steps.
   AGAIN:
   ptIndex = ptIndex + 1L
   from_direction = Find_Boundary_Outline(mask, darray, $
      boundaryPts, ptIndex, xsize, ysize, from_direction)

   IF N_Elements(perimeter) NE 0 THEN BEGIN
      IF N_Elements(scale) EQ 0 THEN BEGIN
         CASE from_direction OF
            0: perimeter = perimeter + 1.0D
            1: perimeter = perimeter + diagonal
            2: perimeter = perimeter + 1.0D
            3: perimeter = perimeter + diagonal
            4: perimeter = perimeter + 1.0D
            5: perimeter = perimeter + diagonal
            6: perimeter = perimeter + 1.0D
            7: perimeter = perimeter + diagonal
            ELSE: perimeter = 4
         ENDCASE
       ENDIF ELSE BEGIN
         CASE from_direction OF
            0: perimeter = perimeter + scale[0]
            1: perimeter = perimeter + diagonal
            2: perimeter = perimeter + scale[1]
            3: perimeter = perimeter + diagonal
            4: perimeter = perimeter + scale[0]
            5: perimeter = perimeter + diagonal
            6: perimeter = perimeter + scale[1]
            7: perimeter = perimeter + diagonal
            ELSE: perimeter = (2*scale[0]) + (2*scale[1])
         ENDCASE
      ENDELSE
   ENDIF
   IF ptIndex EQ 1 THEN BEGIN
       secPt[0] = boundaryPts[0, ptIndex]
       secPt[1] = boundaryPts[1, ptIndex]
       GOTO, AGAIN
   ENDIF
ENDREP UNTIL (boundaryPts[0,ptIndex] EQ secPt[0] AND boundaryPts[1,ptIndex] EQ secPt[1] AND $
    boundaryPts[0,ptIndex-1] EQ firstPt[0] AND boundaryPts[1,ptIndex-1] EQ firstPt[1])
boundaryPts = boundaryPts[*, 0:ptIndex-2]

   ; Calculate area.

IF N_Elements(scale) EQ 0 THEN BEGIN

   area = N_Elements(i)

      ; Calculate area from the perimeter.
      ; The first point must be the same as the last point. Method
      ; of Russ, p.490, _Image Processing Handbook, 2nd Edition_.

   bx = Double(Reform(boundaryPts[0,*]))
   by = Double(Reform(boundaryPts[1,*]))
   bx = [bx, bx[0]]
   by = [by, by[0]]
   n = N_Elements(bx)
   perim_area = Total( (bx[1:n-1] + bx[0:n-2]) * (by[1:n-1] - by[0:n-2]) ) / 2.


ENDIF ELSE BEGIN

   area = N_Elements(i) * scale[0] * scale[1]

      ; Calculate area from the perimeter.
      ; The first point must be the same as the last point. Method
      ; of Russ, p.490, _Image Processing Handbook, 2nd Edition_.

   bx = Double(Reform(boundaryPts[0,*])) * scale[0]
   by = Double(Reform(boundaryPts[1,*])) * scale[1]
   bx = [bx, bx[0]]
   by = [by, by[0]]
   n = N_Elements(bx)
   perim_area = Total( (bx[1:n-1] + bx[0:n-2]) * (by[1:n-1] - by[0:n-2]) ) / 2.

   boundaryPts = Double(Temporary(boundaryPts))
   boundaryPts[0,*] = boundaryPts[0,*] * scale[0]
   boundaryPts[1,*] = boundaryPts[1,*] * scale[1]
ENDELSE

   ; Calculate the centroid

mask = mask GT 0
totalMass = Total(mask)
xcm = Total( Total(mask, 2) * lindgen(xsize) ) / totalMass
ycm = Total( Total(mask, 1) * lindgen(ysize) ) / totalMass
center = [xcm, ycm]

RETURN, boundaryPts
END ; ------------------------------------------------------------------------------------------