[Scipy-svn] r3918 - in trunk/scipy/ndimage: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Feb 11 20:51:23 EST 2008
Author: chris.burns
Date: 2008-02-11 19:51:18 -0600 (Mon, 11 Feb 2008)
New Revision: 3918
Added:
trunk/scipy/ndimage/_segmenter.py
Removed:
trunk/scipy/ndimage/segmenter.py
Modified:
trunk/scipy/ndimage/tests/test_segment.py
Log:
Make segmenter.py private. Issue warning on import.
Copied: trunk/scipy/ndimage/_segmenter.py (from rev 3917, trunk/scipy/ndimage/segmenter.py)
===================================================================
--- trunk/scipy/ndimage/segmenter.py 2008-02-12 01:25:33 UTC (rev 3917)
+++ trunk/scipy/ndimage/_segmenter.py 2008-02-12 01:51:18 UTC (rev 3918)
@@ -0,0 +1,520 @@
+import math
+import numpy as N
+import scipy.ndimage._segment as S
+
+# Issue warning regarding heavy development status of this module
+import warnings
+_msg = "The segmenter code is under heavy development and therefore public \
+API will change in the future. The NIPY group is actively working on this \
+code, and has every intention of generalizing this for the Scipy community.\
+Use this module minimally, if at all, until it this warning is removed."
+warnings.warn(_msg, UserWarning)
+
+# WARNING: _objstruct data structure mirrors a corresponding data structure
+# in ndImage_Segmenter_structs.h that is built into the _segment.so library.
+# These structs must match!
+_objstruct = N.dtype([('L', 'i'),
+ ('R', 'i'),
+ ('T', 'i'),
+ ('B', 'i'),
+ ('Label', 'i'),
+ ('Area', 'i'),
+ ('cX', 'f'),
+ ('cY', 'f'),
+ ('curveClose', 'i'),
+ ('cXB', 'f'),
+ ('cYB', 'f'),
+ ('bLength', 'f'),
+ ('minRadius', 'f'),
+ ('maxRadius', 'f'),
+ ('aveRadius', 'f'),
+ ('ratio', 'f'),
+ ('compactness', 'f'),
+ ('voxelMean', 'f'),
+ ('voxelVar', 'f'),
+ ('TEM', 'f', 20)]
+ )
+
+
+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=_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=_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=_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(filename):
+ """
+ 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(filename):
+ """
+ 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)
+ # this is for the test CT as spine runs off back of image
+ 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)
+
+def build_d_gauss_kernel(gWidth=21, sigma=1.0):
+
+ """
+ build the derivative of Gaussian kernel for Canny edge filter
+ DGFilter = build_d_gauss_kernel(gWidth, sigma)
+ Inputs:
+ gWdith is width of derivative of Gaussian kernel
+ sigma is sigma term of derivative of Gaussian kernel
+ Output:
+ DGFilter (a struct). Use in Canny filter call
+
+ """
+ kernel = N.zeros((1+2*(gWidth-1)), dtype=float)
+ indices = range(1, gWidth)
+
+ i = 0
+ kernel[gWidth-1] = math.exp(((-i*i)/(2.0 * sigma * sigma)))
+ kernel[gWidth-1] *= -(i / (sigma * sigma))
+ for i in indices:
+ kernel[gWidth-1+i] = math.exp(((-i*i)/(2.0 * sigma * sigma)))
+ kernel[gWidth-1+i] *= -(i / (sigma * sigma))
+ kernel[gWidth-1-i] = -kernel[gWidth-1+i]
+
+ DGFilter= {'kernelSize' : gWidth, 'coefficients': kernel}
+
+ return DGFilter
+
+def build_2d_kernel(aperature=21, hiFilterCutoff=10.0):
+
+ """
+ build flat FIR filter with sinc kernel
+ this is bandpass, but low cutoff is 0.0
+ Use in Sobel and Canny filter edge find as image pre-process
+
+ FIRFilter = build_2d_kernel(aperature, hiFilterCutoff)
+ Inputs:
+ aperature is number of FIR taps in sinc kernel
+ hiFilterCutoff is digital frequency cutoff in range (0.0, 180.0)
+ Output:
+ FIRFilter (a struct)
+
+ """
+
+ rad = math.pi / 180.0
+ HalfFilterTaps = (aperature-1) / 2
+ kernel = N.zeros((aperature), dtype=N.float32)
+ LC = 0.0
+ HC = hiFilterCutoff * rad
+ t2 = 2.0 * math.pi
+ t1 = 2.0 * HalfFilterTaps + 1.0
+ indices = range(-HalfFilterTaps, HalfFilterTaps+1, 1)
+ j = 0
+ for i in indices:
+ if i == 0:
+ tLOW = LC
+ tHIGH = HC
+ else:
+ tLOW = math.sin(i*LC)/i
+ tHIGH = math.sin(i*HC)/i
+ # Hamming window
+ t3 = 0.54 + 0.46*(math.cos(i*t2/t1))
+ t4 = t3*(tHIGH-tLOW)
+ kernel[j] = t4
+ j += 1
+
+ # normalize the kernel
+ sum = kernel.sum()
+ kernel /= sum
+
+ FIRFilter= {'kernelSize' : aperature, 'coefficients': kernel}
+
+ return FIRFilter
+
+
+def build_laws_kernel():
+
+ """
+ build 6 length-7 Law's texture filter masks
+ mask names are: 'L', 'S', 'E', 'W', 'R', 'O'
+
+ LAWSFilter = build_laws_kernel()
+
+ Inputs:
+ None
+
+ Output:
+ LAWSFilter (a struct)
+
+ """
+ aperature = (6, 7)
+ coefficients = N.zeros((aperature), dtype=N.float32)
+ names = ('L', 'E', 'S', 'W', 'R', 'O' )
+
+ coefficients[0, :] = ( 1.0, 6.0, 15.0, 20.0, 15.0, 6.0, 1.0 )
+ coefficients[1, :] = (-1.0, -4.0, -5.0, 0.0, 5.0, 4.0, 1.0 )
+ coefficients[2, :] = (-1.0, -2.0, 1.0, 4.0, 1.0, -2.0, -1.0 )
+ coefficients[3, :] = (-1.0, 0.0, 3.0, 0.0, -3.0, 0.0, 1.0 )
+ coefficients[4, :] = ( 1.0, -2.0, -1.0, 4.0, -1.0, -2.0, 1.0 )
+ coefficients[5, :] = (-1.0, 6.0, -15.0, 20.0, -15.0, 6.0, -1.0 )
+
+ LAWSFilter= {'numKernels' : 6, 'kernelSize' : 7, 'coefficients': coefficients, 'names': names}
+
+ return LAWSFilter
+
+def build_morpho_thin_masks():
+
+ """
+ build 2 sets (J and K) of 8 3x3 morphology masks (structuring elements)
+ to implement thinning (medial axis transformation - MAT)
+
+ MATFilter = build_morpho_thin_masks()
+
+ Inputs:
+ None
+
+ Output:
+ MATFilter (a struct)
+
+ """
+
+ # (layers, rows, cols)
+ shape = (8, 3, 3)
+ J_mask = N.zeros((shape), dtype=N.ushort)
+ K_mask = N.zeros((shape), dtype=N.ushort)
+
+ # load the 8 J masks for medial axis transformation
+ J_mask[0][0][0] = 1;
+ J_mask[0][0][1] = 1;
+ J_mask[0][0][2] = 1;
+ J_mask[0][1][1] = 1;
+
+ J_mask[1][0][1] = 1;
+ J_mask[1][1][1] = 1;
+ J_mask[1][1][2] = 1;
+
+ J_mask[2][0][0] = 1;
+ J_mask[2][1][0] = 1;
+ J_mask[2][2][0] = 1;
+ J_mask[2][1][1] = 1;
+
+ J_mask[3][0][1] = 1;
+ J_mask[3][1][0] = 1;
+ J_mask[3][1][1] = 1;
+
+ J_mask[4][0][2] = 1;
+ J_mask[4][1][1] = 1;
+ J_mask[4][1][2] = 1;
+ J_mask[4][2][2] = 1;
+
+ J_mask[5][1][0] = 1;
+ J_mask[5][1][1] = 1;
+ J_mask[5][2][1] = 1;
+
+ J_mask[6][1][1] = 1;
+ J_mask[6][2][0] = 1;
+ J_mask[6][2][1] = 1;
+ J_mask[6][2][2] = 1;
+
+ J_mask[7][1][1] = 1;
+ J_mask[7][1][2] = 1;
+ J_mask[7][2][1] = 1;
+
+
+ # load the 8 K masks for medial axis transformation
+ K_mask[0][2][0] = 1;
+ K_mask[0][2][1] = 1;
+ K_mask[0][2][2] = 1;
+
+ K_mask[1][1][0] = 1;
+ K_mask[1][2][0] = 1;
+ K_mask[1][2][1] = 1;
+
+ K_mask[2][0][2] = 1;
+ K_mask[2][1][2] = 1;
+ K_mask[2][2][2] = 1;
+
+ K_mask[3][1][2] = 1;
+ K_mask[3][2][1] = 1;
+ K_mask[3][2][2] = 1;
+
+ K_mask[4][0][0] = 1;
+ K_mask[4][1][0] = 1;
+ K_mask[4][2][0] = 1;
+
+ K_mask[5][0][1] = 1;
+ K_mask[5][0][2] = 1;
+ K_mask[5][1][2] = 1;
+
+ K_mask[6][0][0] = 1;
+ K_mask[6][0][1] = 1;
+ K_mask[6][0][2] = 1;
+
+ K_mask[7][0][0] = 1;
+ K_mask[7][0][1] = 1;
+ K_mask[7][1][0] = 1;
+
+ MATFilter = {'number3x3Masks' : 8, 'jmask' : J_mask, 'kmask' : K_mask}
+
+ return MATFilter
+
Deleted: trunk/scipy/ndimage/segmenter.py
===================================================================
--- trunk/scipy/ndimage/segmenter.py 2008-02-12 01:25:33 UTC (rev 3917)
+++ trunk/scipy/ndimage/segmenter.py 2008-02-12 01:51:18 UTC (rev 3918)
@@ -1,512 +0,0 @@
-import math
-import numpy as N
-import scipy.ndimage._segment as S
-
-# WARNING: _objstruct data structure mirrors a corresponding data structure
-# in ndImage_Segmenter_structs.h that is built into the _segment.so library.
-# These structs must match!
-_objstruct = N.dtype([('L', 'i'),
- ('R', 'i'),
- ('T', 'i'),
- ('B', 'i'),
- ('Label', 'i'),
- ('Area', 'i'),
- ('cX', 'f'),
- ('cY', 'f'),
- ('curveClose', 'i'),
- ('cXB', 'f'),
- ('cYB', 'f'),
- ('bLength', 'f'),
- ('minRadius', 'f'),
- ('maxRadius', 'f'),
- ('aveRadius', 'f'),
- ('ratio', 'f'),
- ('compactness', 'f'),
- ('voxelMean', 'f'),
- ('voxelVar', 'f'),
- ('TEM', 'f', 20)]
- )
-
-
-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=_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=_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=_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(filename):
- """
- 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(filename):
- """
- 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)
- # this is for the test CT as spine runs off back of image
- 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)
-
-def build_d_gauss_kernel(gWidth=21, sigma=1.0):
-
- """
- build the derivative of Gaussian kernel for Canny edge filter
- DGFilter = build_d_gauss_kernel(gWidth, sigma)
- Inputs:
- gWdith is width of derivative of Gaussian kernel
- sigma is sigma term of derivative of Gaussian kernel
- Output:
- DGFilter (a struct). Use in Canny filter call
-
- """
- kernel = N.zeros((1+2*(gWidth-1)), dtype=float)
- indices = range(1, gWidth)
-
- i = 0
- kernel[gWidth-1] = math.exp(((-i*i)/(2.0 * sigma * sigma)))
- kernel[gWidth-1] *= -(i / (sigma * sigma))
- for i in indices:
- kernel[gWidth-1+i] = math.exp(((-i*i)/(2.0 * sigma * sigma)))
- kernel[gWidth-1+i] *= -(i / (sigma * sigma))
- kernel[gWidth-1-i] = -kernel[gWidth-1+i]
-
- DGFilter= {'kernelSize' : gWidth, 'coefficients': kernel}
-
- return DGFilter
-
-def build_2d_kernel(aperature=21, hiFilterCutoff=10.0):
-
- """
- build flat FIR filter with sinc kernel
- this is bandpass, but low cutoff is 0.0
- Use in Sobel and Canny filter edge find as image pre-process
-
- FIRFilter = build_2d_kernel(aperature, hiFilterCutoff)
- Inputs:
- aperature is number of FIR taps in sinc kernel
- hiFilterCutoff is digital frequency cutoff in range (0.0, 180.0)
- Output:
- FIRFilter (a struct)
-
- """
-
- rad = math.pi / 180.0
- HalfFilterTaps = (aperature-1) / 2
- kernel = N.zeros((aperature), dtype=N.float32)
- LC = 0.0
- HC = hiFilterCutoff * rad
- t2 = 2.0 * math.pi
- t1 = 2.0 * HalfFilterTaps + 1.0
- indices = range(-HalfFilterTaps, HalfFilterTaps+1, 1)
- j = 0
- for i in indices:
- if i == 0:
- tLOW = LC
- tHIGH = HC
- else:
- tLOW = math.sin(i*LC)/i
- tHIGH = math.sin(i*HC)/i
- # Hamming window
- t3 = 0.54 + 0.46*(math.cos(i*t2/t1))
- t4 = t3*(tHIGH-tLOW)
- kernel[j] = t4
- j += 1
-
- # normalize the kernel
- sum = kernel.sum()
- kernel /= sum
-
- FIRFilter= {'kernelSize' : aperature, 'coefficients': kernel}
-
- return FIRFilter
-
-
-def build_laws_kernel():
-
- """
- build 6 length-7 Law's texture filter masks
- mask names are: 'L', 'S', 'E', 'W', 'R', 'O'
-
- LAWSFilter = build_laws_kernel()
-
- Inputs:
- None
-
- Output:
- LAWSFilter (a struct)
-
- """
- aperature = (6, 7)
- coefficients = N.zeros((aperature), dtype=N.float32)
- names = ('L', 'E', 'S', 'W', 'R', 'O' )
-
- coefficients[0, :] = ( 1.0, 6.0, 15.0, 20.0, 15.0, 6.0, 1.0 )
- coefficients[1, :] = (-1.0, -4.0, -5.0, 0.0, 5.0, 4.0, 1.0 )
- coefficients[2, :] = (-1.0, -2.0, 1.0, 4.0, 1.0, -2.0, -1.0 )
- coefficients[3, :] = (-1.0, 0.0, 3.0, 0.0, -3.0, 0.0, 1.0 )
- coefficients[4, :] = ( 1.0, -2.0, -1.0, 4.0, -1.0, -2.0, 1.0 )
- coefficients[5, :] = (-1.0, 6.0, -15.0, 20.0, -15.0, 6.0, -1.0 )
-
- LAWSFilter= {'numKernels' : 6, 'kernelSize' : 7, 'coefficients': coefficients, 'names': names}
-
- return LAWSFilter
-
-def build_morpho_thin_masks():
-
- """
- build 2 sets (J and K) of 8 3x3 morphology masks (structuring elements)
- to implement thinning (medial axis transformation - MAT)
-
- MATFilter = build_morpho_thin_masks()
-
- Inputs:
- None
-
- Output:
- MATFilter (a struct)
-
- """
-
- # (layers, rows, cols)
- shape = (8, 3, 3)
- J_mask = N.zeros((shape), dtype=N.ushort)
- K_mask = N.zeros((shape), dtype=N.ushort)
-
- # load the 8 J masks for medial axis transformation
- J_mask[0][0][0] = 1;
- J_mask[0][0][1] = 1;
- J_mask[0][0][2] = 1;
- J_mask[0][1][1] = 1;
-
- J_mask[1][0][1] = 1;
- J_mask[1][1][1] = 1;
- J_mask[1][1][2] = 1;
-
- J_mask[2][0][0] = 1;
- J_mask[2][1][0] = 1;
- J_mask[2][2][0] = 1;
- J_mask[2][1][1] = 1;
-
- J_mask[3][0][1] = 1;
- J_mask[3][1][0] = 1;
- J_mask[3][1][1] = 1;
-
- J_mask[4][0][2] = 1;
- J_mask[4][1][1] = 1;
- J_mask[4][1][2] = 1;
- J_mask[4][2][2] = 1;
-
- J_mask[5][1][0] = 1;
- J_mask[5][1][1] = 1;
- J_mask[5][2][1] = 1;
-
- J_mask[6][1][1] = 1;
- J_mask[6][2][0] = 1;
- J_mask[6][2][1] = 1;
- J_mask[6][2][2] = 1;
-
- J_mask[7][1][1] = 1;
- J_mask[7][1][2] = 1;
- J_mask[7][2][1] = 1;
-
-
- # load the 8 K masks for medial axis transformation
- K_mask[0][2][0] = 1;
- K_mask[0][2][1] = 1;
- K_mask[0][2][2] = 1;
-
- K_mask[1][1][0] = 1;
- K_mask[1][2][0] = 1;
- K_mask[1][2][1] = 1;
-
- K_mask[2][0][2] = 1;
- K_mask[2][1][2] = 1;
- K_mask[2][2][2] = 1;
-
- K_mask[3][1][2] = 1;
- K_mask[3][2][1] = 1;
- K_mask[3][2][2] = 1;
-
- K_mask[4][0][0] = 1;
- K_mask[4][1][0] = 1;
- K_mask[4][2][0] = 1;
-
- K_mask[5][0][1] = 1;
- K_mask[5][0][2] = 1;
- K_mask[5][1][2] = 1;
-
- K_mask[6][0][0] = 1;
- K_mask[6][0][1] = 1;
- K_mask[6][0][2] = 1;
-
- K_mask[7][0][0] = 1;
- K_mask[7][0][1] = 1;
- K_mask[7][1][0] = 1;
-
- MATFilter = {'number3x3Masks' : 8, 'jmask' : J_mask, 'kmask' : K_mask}
-
- return MATFilter
-
Modified: trunk/scipy/ndimage/tests/test_segment.py
===================================================================
--- trunk/scipy/ndimage/tests/test_segment.py 2008-02-12 01:25:33 UTC (rev 3917)
+++ trunk/scipy/ndimage/tests/test_segment.py 2008-02-12 01:51:18 UTC (rev 3918)
@@ -1,6 +1,6 @@
from scipy.testing import *
-from scipy.ndimage.segmenter import *
+from scipy.ndimage._segmenter import *
inputname = 'slice112.raw'
More information about the Scipy-svn
mailing list