[Scipy-svn] r3541 - trunk/scipy/ndimage
scipy-svn@scip...
scipy-svn@scip...
Thu Nov 15 15:57:14 CST 2007
Author: tom.waite
Date: 2007-11-15 15:56:38 -0600 (Thu, 15 Nov 2007)
New Revision: 3541
Added:
trunk/scipy/ndimage/segmenter.py
Log:
Python segment execution
Added: trunk/scipy/ndimage/segmenter.py
===================================================================
--- trunk/scipy/ndimage/segmenter.py 2007-11-15 06:48:24 UTC (rev 3540)
+++ trunk/scipy/ndimage/segmenter.py 2007-11-15 21:56:38 UTC (rev 3541)
@@ -0,0 +1,293 @@
+
+import numpy as N
+import scipy.ndimage.segment as S
+
+# make sure this is local to use as default
+inputname = 'slice112.raw'
+
+import os
+filename = os.path.join(os.path.split(__file__)[0],inputname)
+
+
+def shen_castan(image, IIRFilter=0.8, scLow=0.3, window=7, lowThreshold=220+2048,
+ highThreshold=600+2048, dust=16):
+ """
+ labeledEdges, ROIList = shen_castan(image, [default])
+
+ implements Shen-Castan edge finding
+
+ Inputs - image, IIR filter, shen_castan_low, window, low_threshold, high_threshold, dust
+ - image is the numarray 2D image
+ - IIR filter is filter parameter for exponential filter
+ - shen_castan_low is edge threshold is range (0.0, 1.0]
+ - window is search window for edge detection
+ - low_ and high_ threshold are density values
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
+ labeledEdges, numberObjects = S.shen_castan_edges(scLow, IIRFilter, window,
+ lowThreshold, highThreshold, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=S.objstruct)
+ # return the bounding box for each connected edge
+ S.get_object_stats(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def sobel(image, sLow=0.3, tMode=1, lowThreshold=220+2048, highThreshold=600+2048, BPHigh=10.0,
+ apearture=21, dust=16):
+ """
+ labeledEdges, ROIList = sobel(image, [default])
+
+ implements sobel magnitude edge finding
+
+ Inputs - image, sobel_low, tMode, low_threshold, high_threshold,
+ high_filter_cutoff, filter_aperature, dust
+ - image is the numarray 2D image
+ - sobel_low is edge threshold is range (0.0, 1.0]
+ - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
+ - low_ and high_ threshold are density values
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - aperature is odd filter kernel length
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
+ # get sobel edge points. return edges that are labeled (1..numberObjects)
+ labeledEdges, numberObjects = S.sobel_edges(sLow, tMode, lowThreshold,
+ highThreshold, BPHigh, apearture, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=S.objstruct)
+ # return the bounding box for each connected edge
+ S.get_object_stats(labeledEdges, ROIList)
+ # thin (medial axis transform) of the sobel edges as the sobel produces a 'band edge'
+ S.morpho_thin_filt(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def canny(image, cSigma=1.0, cLow=0.5, cHigh=0.8, tMode=1, lowThreshold=220+2048,
+ highThreshold=600+2048, BPHigh=10.0, apearture=21, dust=16):
+ """
+ labeledEdges, ROIList = canny(image, [default])
+
+ implements canny edge finding
+
+ Inputs - image, DG_sigma, canny_low, canny_high, tMode, low_threshold,
+ high_threshold, high_filter_cutoff, filter_aperature, dust
+ - image is the numarray 2D image
+ - DG_sigma is Gaussain sigma for the derivative-of-gaussian filter
+ - clow is low edge threshold is range (0.0, 1.0]
+ - chigh is high edge threshold is range (0.0, 1.0]
+ - tMode is threshold mode: 1 for ave, 2 for mode (histogram peak)
+ - low_ and high_ threshold are density values
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - high_filter_cutoff is digital high frequency cutoff in range (0.0, 180.0]
+ - aperature is odd filter kernel length
+ - dust is blob filter. blob area (length x width of bounding box) under this
+ size threshold are filtered (referred to as dust and blown away)
+
+ Outputs - labeledEdges, ROIList[>dust]
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList[>dust] is a blob feature list. Only values
+ with bounding box area greater than dust threshold are returned
+
+ """
+ # get canny edge points. return edges that are labeled (1..numberObjects)
+ labeledEdges, numberObjects = S.canny_edges(cSigma, cLow, cHigh, tMode, lowThreshold, highThreshold,
+ BPHigh, apearture, image)
+ # allocated struct array for edge object measures. for now just the rect bounding box
+ ROIList = N.zeros(numberObjects, dtype=S.objstruct)
+ # return the bounding box for each connected edge
+ S.get_object_stats(labeledEdges, ROIList)
+ return labeledEdges, ROIList[ROIList['Area']>dust]
+
+def get_shape_mask(labeledEdges, ROIList):
+ """
+ get_shape_mask(labeledEdges, ROIList)
+
+ takes labeled edge image plus ROIList (blob descriptors) and generates
+ boundary shape features and builds labeled blob masks. 'labeledEdges'
+ is over-written by 'labeledMask'. Adds features to ROIList structure
+
+ Inputs - labeledEdges, ROIList
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. edge image input is over-written with mask image.
+ ROIList added to.
+
+ """
+
+ # pass in Sobel morph-thinned labeled edge image (LEI) and ROIList
+ # GetShapeMask will augment the ROI list
+ # labeledEdges is the original edge image and overwritten as mask image
+ # maskImage is the mask that is used for blob texture / pixel features
+ S.build_boundary(labeledEdges, ROIList)
+ return
+
+def get_voxel_measures(rawImage, labeledEdges, ROIList):
+ """
+ get_voxel_measures(rawImage, labeledEdges, ROIList)
+
+ takes raw 2D image, labeled blob mask and ROIList. computes voxel features
+ (moments, histogram) for each blob. Adds features to ROIList structure.
+
+ Inputs - rawImage, labeledEdges, ROIList
+ - rawImage is the original source 2D image
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. ROIList added to.
+
+ """
+ #
+ # pass raw image, labeled mask and the partially filled ROIList
+ # VoxelMeasures will fill the voxel features in the list
+ #
+ S.voxel_measures(rawImage, labeledEdges, ROIList)
+ return
+
+def get_texture_measures(rawImage, labeledEdges, ROIList):
+ """
+ get_texture_measures(rawImage, labeledEdges, ROIList)
+
+ takes raw 2D image, labeled blob mask and ROIList. computes 2D
+ texture features using 7x7 Law's texture filters applied
+ to segmented blobs. TEM (texture energy metric) is computed
+ for each Law's filter image and stored in TEM part of ROIList.
+
+ Inputs - rawImage, labeledEdges, ROIList
+ - rawImage is the original source 2D image
+ - labeledEdges is boundary (edges) of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is a blob feature list.
+
+ Output - no return. ROIList added to.
+ """
+ #
+ # pass raw image, labeled mask and the partially filled ROIList
+ # VoxelMeasures will fill the texture (Law's, sub-edges, co-occurence, Gabor) features in the list
+ #
+ S.texture_measures(rawImage, labeledEdges, ROIList)
+ return
+
+def segment_regions():
+ """
+ sourceImage, labeledMask, ROIList = segment_regions()
+
+ Inputs - No Input
+
+ Outputs - sourceImage, labeledMask, ROIList
+ - sourceImage is raw 2D image (default cardiac CT slice for demo
+ - labeledMask is mask of segmented 'blobs',
+ numerically labeled by blob number
+ - ROIList is numerical Python structure of intensity, shape and
+ texture features for each blob
+
+ High level script calls Python functions:
+ get_slice() - a cardiac CT slice demo file
+ sobel() - sobel magnitude edge finder,
+ returns connected edges
+ get_shape_mask() - gets segmented blob boundary and mask
+ and shape features
+ get_voxel_measures() - uses masks get object voxel moment
+ and histogram features
+ get_texture_measures() - uses masks get object 2D texture features
+ """
+ # get slice from the CT volume
+ image = get_slice(filename)
+ # need a copy of original image as filtering will occur on the extracted slice
+ sourceImage = image.copy()
+ # Sobel is the first level segmenter. Sobel magnitude and MAT (medial axis transform)
+ # followed by connected component analysis. What is returned is labeled edges and the object list
+ labeledMask, ROIList = sobel(image)
+ # From the labeled edges and the object list get the labeled mask for each blob object
+ get_shape_mask(labeledMask, ROIList)
+ # Use the labeled mask and source image (raw) to get voxel features
+ get_voxel_measures(sourceImage, labeledMask, ROIList)
+ # Use the labeled mask and source image (raw) to get texture features
+ get_texture_measures(sourceImage, labeledMask, ROIList)
+ return sourceImage, labeledMask, ROIList
+
+def grow_regions():
+ """
+ regionMask, numberRegions = region_grow()
+ Inputs - No Input
+ Outputs - regionMask, numberRegions
+ - regionMask is the labeled segment masks from 2D image
+ - numberRegions is the number of segmented blobs
+
+ High level script calls Python functions:
+ get_slice() - a cardiac CT slice demo file
+ region_grow() - "grows" connected blobs. default threshold
+ and morphological filter structuring element
+ """
+ # get slice from the CT volume
+ image = get_slice(filename)
+ regionMask, numberRegions = region_grow(image)
+ return regionMask, numberRegions
+
+
+def region_grow(image, lowThreshold=220+2048, highThreshold=600+2048, open=7, close=7):
+ """
+ regionMask, numberRegions = region_grow(image, [defaults])
+
+ Inputs - image, low_threshold, high_threshold, open, close
+ - image is the numarray 2D image
+ - low_ and high_ threshold are density values
+ - open is open morphology structuring element
+ odd size. 0 to turn off. max is 11
+ - close is close morphology structuring element
+ odd size. 0 to turn off. max is 11
+
+ Outputs - regionMask, numberRegions
+ - regionMask is the labeled segment masks from 2D image
+ - numberRegions is the number of segmented blobs
+ """
+ # morphology filters need to be clipped to 11 max and be odd
+ regionMask, numberRegions = S.region_grow(lowThreshold, highThreshold, close, open, image)
+ return regionMask, numberRegions
+
+
+def get_slice(imageName='slice112.raw', bytes=2, rows=512, columns=512):
+ # get a slice alrady extracted from the CT volume
+ #image = open(imageName, 'rb')
+ #slice = image.read(rows*columns*bytes)
+ #values = struct.unpack('h'*rows*columns, slice)
+ #ImageSlice = N.array(values, dtype=float).reshape(rows, columns)
+
+ ImageSlice = N.fromfile(imageName, dtype=N.uint16).reshape(rows, columns);
+
+ # clip the ends for this test CT image file as the spine runs off the end of the image
+ ImageSlice[505:512, :] = 0
+ return (ImageSlice).astype(float)
+
+def get_slice2(image_name='slice112.raw', bytes=2, shape=(512,512)):
+ import mmap
+ file = open(image_name, 'rb')
+ mm = mmap.mmap(file.fileno(), 0, access=mmap.ACCESS_READ)
+ slice = N.frombuffer(mm, dtype='u%d' % bytes).reshape(shape)
+ slice = slice.astype(float)
+ slice[505:512,:] = 0
+ return slice
+
+def save_slice(mySlice, filename='junk.raw', bytes=4):
+ # just save the slice to a fixed file
+ slice = mySlice.astype('u%d' % bytes)
+ slice.tofile(filename)
+
+
More information about the Scipy-svn
mailing list