;+ ; 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. ;- ;******************************************************************************************; ; Copyright (c) 2008, 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 ; 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 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 ENDREP UNTIL (boundaryPts[0,ptIndex] EQ firstPt[0] AND $ boundaryPts[1,ptIndex] EQ firstPt[1]) boundaryPts = boundaryPts[*,0:ptIndex-1] ; 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) * Indgen(xsize) ) / totalMass ycm = Total( Total(mask, 1) * Indgen(ysize) ) / totalMass center = [xcm, ycm] RETURN, boundaryPts END ; ------------------------------------------------------------------------------------------