| Re: Chain-Link Algorithm for Perimeter [message #30322] |
Sat, 20 April 2002 17:01  |
Ted Cary
Messages: 53
Registered: October 2001
|
Member |
|
|
David Fanning wrote
>
> I think what I want to use is the Chain-Link Coding
> algorithm of R.L.T. Cederberg as described in The
> Image Processing Handbook by John Russ. Has anyone
> coded this up in IDL by any chance?
>
> I've been known to pay for code that saves me a ton
> of time. :-)
Hi David,
I coded up a routine that should order the boundary points of your blobs
for you. I don't know if it uses the chain-link algorithm. It's based
on C++ code by Dave Eberly of Magic Software. It will handle
concavities and folds but not all bow ties, at least in my tests.
It took me a good part of the day, but Eberly and IDL did most of the
work. I've been meaning to code something like this for myself anyway
and your request just got me started. It's still not the UNMASK routine
I was looking for in an earlier post, but it gets the job done.
The routine is called simply "BOUNDARY" and is included in this post as
both an attachment and as text below. See the header for instructions.
Enjoy,
TC
;+
; NAME:
; BOUNDARY
;
; PURPOSE:
;
; Returns ordered lists of x and y coordinates of ALL points on the
boundary of a region.
; Will work on regions with concavities and regions whose boundaries
have folds.
;
; AUTHOR:
;
; Ted Cary
; adapted from C++ code by David Eberly
; (see www.magic_software.com/ImageAnalysis.html)
;
; CATEGORY:
;
; Image Analysis/Contouring.
;
; CALLING SEQUENCE:
;
; BOUNDARY, mask, X, Y
;
; MASK : A (binary) 2D array of the type put out by
IDLanROI::ComputeMask.
; The interior points of the region to be contoured should be
set
; to 1, while background points should be set to 0.
;
; X : An output vector that will contain the x coordinates of the
boundary points.
;
; Y : An output vector that will contain the y coordinates of the
boundary points.
;
; EXAMPLE:
;
; Find boundary of blob made of overlapping disks::
;
; ; Create disk of radius r
;
; r = 50
; disk = Shift(Dist(2*r+1), r, r) LT r
;
; ; Create "blob" of overlapping disks
;
; mask1 = BytArr(275, 200)
; mask2 = BytArr(275, 200)
; mask1[49, 49] = disk
; mask2[124, 49] = disk
; mask = mask1 OR mask2 ; overlap two disks
;
; ; Show blob
;
; TVSCL, mask
;
; ; Get boundary
;
; Boundary, mask, x, y
;
; ; Plot Boundary in red (using 24-bit color)
;
; PlotS, x, y, Color=255, /Device
;
; LICENSE:
;
; Please give me and Dave Eberly of Magic Software credit for the code
and
; note any changes you make to it. Don't remove this notice, etc...
;
; "AS IS", no warranty, express or implied. Authors are not liable.
;
; Enjoy.
;
;########################################################### ################
PRO Boundary, mask, X, Y
; Create distance map of mask, padded by two on each side.
dims = Size(mask, /Dimensions)
distMap = BytArr(dims[0] + 4 , dims[1] + 4)
distMap[2,2] = Morph_Distance(mask, Neighbor_Sampling=0)
; Set initial coordinate of starting boundary point.
startIndex = Max(Where(distMap EQ 1))
x0 = startIndex MOD (dims[0] + 4)
y0 = startIndex/(dims[0] + 4)
; Set up x, y directional tables.
dX = [-1, 0, +1, +1, +1, 0, -1, -1]
dY = [-1, -1, -1, 0, +1, +1, +1, 0]
; Determine direction from background to start point.
iCx = x0
iCy = y0
FOR idir = 0, 7 DO BEGIN
iNx = iCx + dX[idir]
iNy = iCy + dY[idir]
IF distMap[iNx, iNy] NE 0 THEN BEGIN
idir = (idir + 1) MOD 8
BREAK
ENDIF
ENDFOR
; Initialize arrays to hold found boundary vertex coordinates.
xList = [x0 - 2]
yList = [y0 - 2]
; Traverse boundary in CW order.
; Mark traversed points as -1 in distance map.
distMap[x0, y0] = -1 ; Already visited first point.
WHILE(1) DO BEGIN
iNbr = iDir
FOR i=0, 7 DO BEGIN
iNx = iCx + dx[iNbr]
iNy = iCy + dy[iNbr]
iNbr = (iNbr + 1) MOD 8
IF distMap[iNx, iNy] EQ 1 THEN BREAK
ENDFOR
IF i EQ 8 THEN BREAK
IF iNx EQ x0 AND iNy EQ y0 THEN BREAK
; Update vertex coords lists.
xList = [xList, iNx-2]
yList = [yList, iNy-2]
; Mark visited pixels.
distMap[iNx, iNy] = -1
; Set new point coords, new direction.
iCx = iNx
iCy = iNy
iDir = (i + 5 + iDir) MOD 8
ENDWHILE
; Output vectors.
x = xList
y = yList
END;--------------------------------------------------------
;+
; NAME:
; BOUNDARY
;
; PURPOSE:
;
; Returns *ordered* lists of x and y coordinates of ALL points on the boundary a region.
; Will work on regions with concavities and regions whose boundaries have folds.
;
; AUTHOR:
;
; Ted Cary
; adapted from C++ code by David Eberly
; (see www.magic_software.com/ImageAnalysis.html)
;
; CATEGORY:
;
; Image Analysis/Contouring.
;
; CALLING SEQUENCE:
;
; BOUNDARY, mask, X, Y
;
; MASK : A (binary) 2D array of the type put out by IDLanROI::ComputeMask.
; The interior points of the region to be contoured should be set
; to 1, while background points should be set to 0.
;
; X : An output vector that will contain the x coordinates of the boundary points.
;
; Y : An output vector that will contain the y coordinates of the boundary points.
;
; EXAMPLE:
;
; Find boundary of blob made of overlapping disks::
;
; ; Create disk of radius r
;
; r = 50
; disk = Shift(Dist(2*r+1), r, r) LT r
;
; ; Create "blob" of overlapping disks
;
; mask1 = BytArr(275, 200)
; mask2 = BytArr(275, 200)
; mask1[49, 49] = disk
; mask2[124, 49] = disk
; mask = mask1 OR mask2 ; overlap two disks
;
; ; Show blob
;
; TVSCL, mask
;
; ; Get boundary
;
; Boundary, mask, x, y
;
; ; Plot Boundary in red (using 24-bit color)
;
; PlotS, x, y, Color=255, /Device
;
; LICENSE:
;
; Please give me and Dave Eberly of Magic Software credit for the code and
; note any changes you make to it. Don't remove this notice, etc...
;
; "AS IS", no warranty, express or implied. Authors are not liable.
;
; Enjoy.
;
;########################################################### ################
PRO Boundary, mask, X, Y
; Create distance map of mask, padded by two on each side.
dims = Size(mask, /Dimensions)
distMap = BytArr(dims[0] + 4 , dims[1] + 4)
distMap[2,2] = Morph_Distance(mask, Neighbor_Sampling=0)
; Set initial coordinate of starting boundary point.
startIndex = Max(Where(distMap EQ 1))
x0 = startIndex MOD (dims[0] + 4)
y0 = startIndex/(dims[0] + 4)
; Set up x, y directional tables.
dX = [-1, 0, +1, +1, +1, 0, -1, -1]
dY = [-1, -1, -1, 0, +1, +1, +1, 0]
; Determine direction from background to start point.
iCx = x0
iCy = y0
FOR idir = 0, 7 DO BEGIN
iNx = iCx + dX[idir]
iNy = iCy + dY[idir]
IF distMap[iNx, iNy] NE 0 THEN BEGIN
idir = (idir + 1) MOD 8
BREAK
ENDIF
ENDFOR
; Initialize arrays to hold found boundary vertex coordinates.
xList = [x0 - 2]
yList = [y0 - 2]
; Traverse boundary in CW order.
; Mark traversed points as -1 in distance map.
distMap[x0, y0] = -1 ; Already visited first point.
WHILE(1) DO BEGIN
iNbr = iDir
FOR i=0, 7 DO BEGIN
iNx = iCx + dx[iNbr]
iNy = iCy + dy[iNbr]
iNbr = (iNbr + 1) MOD 8
IF distMap[iNx, iNy] EQ 1 THEN BREAK
ENDFOR
IF i EQ 8 THEN BREAK
IF iNx EQ x0 AND iNy EQ y0 THEN BREAK
; Update vertex coords lists.
xList = [xList, iNx-2]
yList = [yList, iNy-2]
; Mark visited pixels.
distMap[iNx, iNy] = -1
; Set new point coords, new direction.
iCx = iNx
iCy = iNy
iDir = (i + 5 + iDir) MOD 8
ENDWHILE
; Output vectors.
x = xList
y = yList
END;--------------------------------------------------------
-
Attachment: boundary.pro
(Size: 3.17KB, Downloaded 100 times)
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
