;+ ; NAME: ; FIND_BOUNDARY ; ; PURPOSE: ; ; This program finds the boundary points about a region of interest (ROI) ; represented by pixel indices. ; ; AUTHOR: ; ; FANNING SOFTWARE CONSULTING ; David Fanning, Ph.D. ; 1645 Sheely Drive ; Fort Collins, CO 80526 USA ; Phone: 970-221-0438 ; E-mail: davidf@dfanning.com ; Coyote's Guide to IDL Programming: http://www.dfanning.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-npoints array of the X and Y points that describe the ; boundary. The points are in the device coodinate system. ; ; INPUT KEYWORDS: ; ; SCALE - A one-element or two-element array of the pixel scale factors, [xscale, yscale], ; used to calculate the perimeter length of the boundary. The SCALE keyword is ; used only if the PERIMETER keyword is provided as an output variable. If ; one-element, the same scale factor is used in both the X and Y direction. ; Set to [1.0,1.0] by default. ; ; 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 area enclosed by the boundary points. The ; area is scaled with the SCALE scaling factors. Note that this area can be larger ; than the area of the ROI if the ROI has interior "holes". ; ; PERIMETER - A named variable that will contain the perimeter length of the boundary ; upon returning from the function. ; ; 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=FSC_Color('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. ;- ; ;########################################################################### ; ; LICENSE ; ; This software is OSI Certified Open Source Software. ; OSI Certified is a certification mark of the Open Source Initiative. ; ; Copyright © 1997-2002 Fanning Software Consulting. ; ; This software is provided "as-is", without any express or ; implied warranty. In no event will the authors be held liable ; for any damages arising from the use of this software. ; ; Permission is granted to anyone to use this software for any ; purpose, including commercial applications, and to alter it and ; redistribute it freely, subject to the following restrictions: ; ; 1. The origin of this software must not be misrepresented; you must ; not claim you wrote the original software. If you use this software ; in a product, an acknowledgment in the product documentation ; would be appreciated, but is not required. ; ; 2. Altered source versions must be plainly marked as such, and must ; not be misrepresented as being the original software. ; ; 3. This notice may not be removed or altered from any source distribution. ; ; For more information on Open Source Software, visit the Open Source ; web site: http://www.opensource.org. ; ;########################################################################### FUNCTION Find_Boundary_ERROR_MESSAGE, theMessage, Traceback=traceback, $ NoName=noName, _Extra=extra On_Error, 2 ; Check for presence and type of message. IF N_Elements(theMessage) EQ 0 THEN theMessage = !Error_State.Msg s = Size(theMessage) messageType = s[s[0]+1] IF messageType NE 7 THEN BEGIN Message, "The message parameter must be a string.", _Extra=extra ENDIF ; Get the call stack and the calling routine's name. Help, Calls=callStack callingRoutine = (Str_Sep(StrCompress(callStack[1])," "))[0] ; Are widgets supported? Doesn't matter in IDL 5.3 and higher. widgetsSupported = ((!D.Flags AND 65536L) NE 0) OR Float(!Version.Release) GE 5.3 IF widgetsSupported THEN BEGIN IF Keyword_Set(noName) THEN answer = Dialog_Message(theMessage, _Extra=extra) ELSE BEGIN IF StrUpCase(callingRoutine) EQ "$MAIN$" THEN answer = Dialog_Message(theMessage, _Extra=extra) ELSE $ answer = Dialog_Message(StrUpCase(callingRoutine) + ": " + theMessage, _Extra=extra) ENDELSE ENDIF ELSE BEGIN Message, theMessage, /Continue, /NoPrint, /NoName, /NoPrefix, _Extra=extra Print, '%' + callingRoutine + ': ' + theMessage answer = 'OK' ENDELSE ; Provide traceback information if requested. IF Keyword_Set(traceback) THEN BEGIN Help, /Last_Message, Output=traceback Print,'' Print, 'Traceback Report from ' + StrUpCase(callingRoutine) + ':' Print, '' FOR j=0,N_Elements(traceback)-1 DO Print, " " + traceback[j] ENDIF RETURN, answer END ; ------------------------------------------------------------------------------------------ FUNCTION Find_Boundary_Outline, mask, darray, boundaryPts, ptIndex, $ xsize, ysize, from_direction 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, XSize=xsize, YSize=ysize, $ Perimeter=perimeter, Scale=scale, Area=area Catch, theError IF theError NE 0 THEN BEGIN Catch, /Cancel ok = Find_Boundary_Error_Message(/Traceback) 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 scale = [1.0,1.0] IF N_Elements(scale) EQ 1 THEN scale = [Double(scale), Double(scale)] diagonal = SQRT(scale[0]^2 + scale[1]^2) IF N_Elements(xsize) EQ 0 THEN xsize = !D.X_Size IF N_Elements(ysize) EQ 0 THEN ysize = !D.Y_Size 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. boundaryPts = IntArr(2, Long(xsize) * ysize / 4L) 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 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 ENDIF ENDREP UNTIL (boundaryPts[0,ptIndex] EQ firstPt[0] AND $ boundaryPts[1,ptIndex] EQ firstPt[1]) boundaryPts = boundaryPts[*,0:ptIndex-1] ; Need an area? IF Arg_Present(area) THEN BEGIN ; Must have at least three points for a polygon. IF N_Elements(boundaryPts) GT 6 THEN BEGIN ; Still may not have a polygon (may have a line). ; Turn POLYFILLV warning messages off. Catch the ; problem yourself. quiet = !Quiet !Quiet = 1 i = PolyFillV(boundaryPts[0,*], boundaryPts[1,*], xsize, ysize) !Quiet = quiet IF i[0] EQ -1 THEN $ area = (N_Elements(boundaryPts) / 2 + 2) * scale[0] * scale[1] ELSE $ ; A line. area = N_Elements(i) * scale[0] * scale[1] ; Normal polygon. ENDIF ELSE BEGIN area = (N_Elements(boundaryPts) / 2 + 2) * scale[0] * scale[1] ; A line. ENDELSE ENDIF RETURN, boundaryPts END ; ------------------------------------------------------------------------------------------