; $Id: region_grow.pro,v 1.7 2004/03/30 19:26:08 griz Exp $ ; ; Copyright (c) 2001-2004, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ; ;+ ; NAME: ; REGION_GROW ; ; PURPOSE: ; This function performs region growing for a given region within ; an N-dimensional array. REGION_GROW finds all pixels within the ; array that are connected neighbors to the region pixels and that ; fall within provided constraints. ; ; The constraints are specified either as a threshold range (a ; minimum and maximum pixel value) or as a multiple of the standard ; deviation of the region pixel values. ; ; If the threshold is used, the region is grown to include all ; connected neighboring pixels that fall within the given threshold. ; ; If the standard deviation multiple is used, the region is grown ; to include all connected neighboring pixels that fall within the ; range of the mean (of the region's pixel values) plus or minus the ; multiple times the standard deviation. ; ; CATEGORY: ; Image Processing. ; ; CALLING SEQUENCE: ; Result = REGION_GROW(Array, ROIPixels) ; ; INPUTS: ; Array: An N-dimensional array of data values. The region will ; be grown according to the data values within this array. ; ; ROIPixels: A vector of indices into Array that represent ; the initial region that is to be grown. ; ; KEYWORD PARAMETERS: ; ALL_NEIGHBORS: Set this keyword to indicate that all adjacent ; neighbors to a given pixel should be considered during ; region growing. (This is sometimes called 8-neighbor ; searching when the array is 2-dimensional.) The default ; is to search only the neighbors that are exactly one unit ; in distance from the current pixel (sometimes called ; 4-neighbor searching when the array is 2-dimensional). ; ; NAN: Set this keyword to a non-zero value to cause the routine ; to check for occurrences of the IEEE floating-point value, ; NaN, in the input array. Elements with the value NaN are ; treated as missing data. ; ; STDDEV_MULTIPLIER: Set this keyword to a scalar value that ; serves as the multiplier of the standard deviation of the ; original region pixel values. The pixel values of the ; grown region must fall within the range of: ; ; Mean +/- StdDevMultiplier*StdDev ; ; This keyword is mutually exclusive of THRESHOLD. ; ; THRESHOLD: Set this keyword to a two-element vector, [min,max], ; of the inclusive range within which the pixel values of the ; grown region must fall. The default is the range of pixel ; values within the initial region. This keyword is mutually ; exclusive of STDDEV_MULTIPLIER. ; ; OUTPUTS: ; This function returns the vector of indices into Array that represent ; pixels within the grown region. (Note: the grown region will not ; include pixels at the edges of the input array.) If no pixels fall ; within the grown region, this function returns the value -1. ; ; EXAMPLE: ; Grow a pre-defined region within an image of human red blood cells. ; ; ; Load an image. ; fname = FILEPATH('rbcells.jpg', SUBDIR=['examples','data']) ; READ_JPEG, fname, img ; imgDims = SIZE(img, /DIMENSIONS) ; ; ; Define original region pixels. ; x = FINDGEN(16*16) MOD 16 + 276. ; y = LINDGEN(16*16) / 16 + 254. ; roiPixels = x + y * imgDims[0] ; ; ; Grow the region. ; newROIPixels = REGION_GROW(img, roiPixels) ; ; ; Load a greyscale color table. ; LOADCT, 0 ; ; ; Set the topmost color table entry to red. ; topClr = !D.TABLE_SIZE-1 ; TVLCT, 255, 0, 0, topClr ; ; ; Show the results. ; tmpImg = BYTSCL(img, TOP=(topClr-1)) ; tmpImg[rOIPixels] = topClr ; WINDOW, 0, XSIZE=imgDims[0], YSIZE=imgDims[1], $ ; TITLE='Original Region' ; TV, tmpImg ; ; tmpImg = BYTSCL(img, TOP=(topClr-1)) ; tmpImg[newROIPixels] = topClr ; WINDOW, 2, XSIZE=imgDims[0], YSIZE=imgDims[1], $ ; TITLE='Grown Region' ; TV, tmpImg ; ; MODIFICATION HISTORY: ; Written by: DLD, February 2001. ; ;- function REGION_GROW, array, roiPixels, $ ALL_NEIGHBORS=allNeighbors, $ NAN=nan, $ STDDEV_MULTIPLIER=stdDevMult, $ THRESHOLD=threshold arrDims = SIZE(array, /DIMENSIONS) ; Extract the image pixel values that fall within the input region. roiPixelData = array[roiPixels] nROIPixels = N_ELEMENTS(roiPixels) if (TOTAL(FINITE(roiPixelData)) eq 0) then $ MESSAGE, $ 'Unable to grow region that consists entirely of missing data.' ; Set the initial threshold range for the grown region. if (N_ELEMENTS(threshold) eq 2) then begin thresh_lo = threshold[0] thresh_hi = threshold[1] endif else $ thresh_lo = MIN(roiPixelData, MAX=thresh_hi, NAN=nan) ; If requested, compute a threshold range based on the standard ; deviation and mean. if (N_ELEMENTS(stdDevMult) ne 0) then begin if (N_ELEMENTS(threshold) eq 2) then begin errStr = 'STDDEV_MULTIPLIER and THRESHOLD are mutually ' + $ 'exclusive. Using the THRESHOLD value.' MESSAGE, /CONTINUE, errStr endif else begin if (nROIPixels gt 1) then begin mean = MEAN(roiPixelData, /DOUBLE, NAN=nan) sdev = STDDEV(roiPixelData, /DOUBLE, NAN=nan) endif else begin mean = roiPixelData sdev = 0.0d endelse offset = stdDevMult * sdev if (SIZE(array, /TYPE) le 3) then begin thresh_lo = ROUND(mean - offset) thresh_hi = ROUND(mean + offset) endif else begin thresh_lo = mean - offset thresh_hi = mean + offset endelse endelse endif ; Threshold the original image. if (KEYWORD_SET(nan)) then begin ; Temporarily overwrite array NaN values with 0s. iNan = WHERE(~FINITE(array), nNan) if (nNan gt 0) then begin oldNanVals = array[iNan] array[iNan] = 0 endif endif else $ nNan = 0 threshArray = (array ge thresh_lo and array le thresh_hi) ; Replace NaN values in input array, and set result of threshArray ; accordingly. if (nNan gt 0) then begin array[iNan] = oldNanVals threshArray[iNan] = 0b endif ; Label the regions within the thresholded image. labelArray = LABEL_REGION(threshArray, ALL_NEIGHBORS=allNeighbors) ; Determine which labels fall within the ROI. if (nROIPixels gt 1) then begin labels = WHERE(HISTOGRAM(labelArray[roiPixels], MIN=0) ne 0, nLabels) endif else begin nLabels = 1 labels = labelArray[roiPixels] endelse ; Compute a total of all labeled pixels for the labels within the ROI ; (excluding pixels with a label of 0). labelHist = HISTOGRAM(labelArray, REVERSE=r, MIN=0) nPixels = TOTAL(labelHist[labels]) - (labelHist[labels[0]] * $ (labels[0] eq 0)) ; Create a new array of pixel indices for the labeled pixels. if (nPixels gt 0) then begin growROIPixels = LONARR(nPixels) j = 0 for i=0L, nLabels-1 do begin if (r[labels[i]+1] le N_ELEMENTS(r) and labels[i] ne 0) then begin growROIPixels[j] = r[r[labels[i]]:r[labels[i]+1]-1] j = j + labelHist[labels[i]] endif endfor endif else $ growROIPixels = -1L return, growROIPixels end