[Scipy-svn] r4047 - trunk/scipy/ndimage
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Mar 20 13:44:28 EDT 2008
Author: tom.waite
Date: 2008-03-20 12:44:00 -0500 (Thu, 20 Mar 2008)
New Revision: 4047
Modified:
trunk/scipy/ndimage/_segmenter.py
Log:
updates and improvements
Modified: trunk/scipy/ndimage/_segmenter.py
===================================================================
--- trunk/scipy/ndimage/_segmenter.py 2008-03-20 17:43:13 UTC (rev 4046)
+++ trunk/scipy/ndimage/_segmenter.py 2008-03-20 17:44:00 UTC (rev 4047)
@@ -1,29 +1,28 @@
import math
-import numpy as N
-import pylab as P
+import numpy as NP
import scipy.ndimage._segment as S
-_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)]
- )
+_objstruct = NP.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)]
+ )
# Issue warning regarding heavy development status of this module
@@ -57,7 +56,7 @@
"""
[rows, cols] = magnitude.shape
- edge_image = N.zeros(rows*cols, dtype=N.int16).reshape(rows, cols)
+ edge_image = NP.zeros(rows*cols, dtype=NP.int16).reshape(rows, cols)
S.canny_hysteresis(magnitude, edge_image, canny_stats['low'], canny_stats['high'])
return edge_image
@@ -107,7 +106,7 @@
"""
[rows, cols] = horz_DGFilter.shape
- magnitude = N.zeros(rows*cols, dtype=N.float64).reshape(rows, cols)
+ magnitude = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
aveMag, canny_low, canny_high = S.canny_nonmax_supress(horz_DGFilter, vert_DGFilter,
magnitude, mode, img_means['x-dg']*thres,
img_means['y-dg']*thres, canny_l, canny_h)
@@ -145,9 +144,10 @@
"""
+ slice = slice.astype(NP.float64)
[rows, cols] = slice.shape
- horz_DGFilter = N.zeros(rows*cols, dtype=N.float64).reshape(rows, cols)
- vert_DGFilter = N.zeros(rows*cols, dtype=N.float64).reshape(rows, cols)
+ horz_DGFilter = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
+ vert_DGFilter = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
aveX, aveY = S.canny_filter(slice, horz_DGFilter, vert_DGFilter,
dg_kernel['coefficients'], dg_kernel['kernelSize'])
@@ -184,7 +184,7 @@
"""
if ROI==None:
- ROIList = N.zeros(1, dtype=_objstruct)
+ ROIList = NP.zeros(1, dtype=_objstruct)
[rows, cols] = label_image.shape
ROIList['L'] = 2
ROIList['R'] = cols-3
@@ -193,15 +193,15 @@
[rows, cols] = label_image.shape
# destination image
- thin_edge_image = N.zeros(rows*cols, dtype=N.uint16).reshape(rows, cols)
- mat_image = N.zeros(rows*cols, dtype=N.uint16).reshape(rows, cols)
+ thin_edge_image = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
+ mat_image = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
# scratch memory for thin
- input = N.zeros(rows*cols, dtype=N.uint8).reshape(rows, cols)
- cinput = N.zeros(rows*cols, dtype=N.uint8).reshape(rows, cols)
- erosion = N.zeros(rows*cols, dtype=N.uint8).reshape(rows, cols)
- dialation = N.zeros(rows*cols, dtype=N.uint8).reshape(rows, cols)
- hmt = N.zeros(rows*cols, dtype=N.uint8).reshape(rows, cols)
- copy = N.zeros(rows*cols, dtype=N.uint8).reshape(rows, cols)
+ input = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
+ cinput = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
+ erosion = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
+ dialation = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
+ hmt = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
+ copy = NP.zeros(rows*cols, dtype=NP.uint8).reshape(rows, cols)
number_regions = ROI.size
indices = range(0, number_regions)
@@ -247,6 +247,147 @@
return mat_image
+def texture_filter(raw_image, label_image, laws_kernel, ROI=None, verbose=0):
+ """
+ texture_images = texture_filter(raw_image, label_image, laws_kernel, ROI=None, verbose=1)
+ .
+ OR
+ .
+ texture_filter(raw_image, label_image, laws_kernel, ROI=None, verbose=0)
+
+ Parameters
+ ..........
+
+ raw_image : {nd_array}
+ raw double image
+
+ label_image : {nd_array}
+ an image with labeled regions from get_blobs() method
+
+ laws_kernel : {dictionary}
+ set of 6 length-7 Law's texture feature kernels
+
+ ROI : {dictionary}
+ Region of Interest structure that has blob bounding boxes
+
+ verbose : {0, 1}, optional
+ determines if return is to include Law's filter images
+
+ Returns
+ ..........
+
+ laws_image : {dictionary}
+ contains 21 Laws filtered regions for each ROI
+ returned if verbose=1
+
+
+ """
+ if ROI==None:
+ ROI= NP.zeros(1, dtype=_objstruct)
+ [rows, cols] = label_image.shape
+ ROI['L'] = 2
+ ROI['R'] = cols-3
+ ROI['B'] = 2
+ ROI['T'] = rows-3
+
+ laws_image_list = {}
+ number_regions = ROI.size
+ layers = laws_kernel['filters']
+ indices = range(0, number_regions)
+ filters = range(1, layers)
+ for i in indices:
+ left = ROI[i]['L']
+ right = ROI[i]['R']
+ bottom = ROI[i]['B']
+ top = ROI[i]['T']
+ Label = ROI[i]['Label']
+ rows = top-bottom
+ cols = right-left
+ label_region = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
+ source_region = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
+ laws_block = NP.zeros(layers*rows*cols, dtype=NP.float32).reshape(layers, rows, cols)
+ # load the labeled region
+ label_region[0:rows, 0:cols][label_image[bottom:top, left:right]==Label] = 1
+ source_region[0:rows, 0:cols] = raw_image[bottom:top, left:right]
+ S.laws_texture_metric(label_region, source_region, laws_block, laws_kernel['numKernels'],
+ laws_kernel['kernelSize'], laws_kernel['filters'],
+ laws_kernel['coefficients'][0], laws_kernel['coefficients'][1],
+ laws_kernel['coefficients'][2], laws_kernel['coefficients'][3],
+ laws_kernel['coefficients'][4], laws_kernel['coefficients'][5])
+
+ for j in filters:
+ # compute the energy measure for each filter in the ROI
+ mask_image = laws_block[j, :, :][label_region[:, :]>0]
+ mean = mask_image.mean()
+ # TBD: scale by mean, or autoscale the TEM vector
+ ROI[i]['TEM'][j-1] = mask_image.std()
+
+ # accumulate the 21 Law's filtered ROI's and optional
+ # return as image (3D)
+ laws_image_list[i] = laws_block
+
+ if verbose == 1:
+ return laws_image_list
+ else:
+ return
+
+
+
+def get_voxel_measures(label_image, raw_image, ROI=None):
+ """
+ mat_image = mat_filter(label_image, raw_image, ROI=None)
+
+ takes the ROI dictionary with the blob bounding boxes and gets the voxel measures
+ from each ROI in the raw data.
+
+ Parameters
+ ..........
+
+ label_image : {nd_array}
+ an image with labeled regions from get_blobs() method
+
+ raw_image : {nd_array}
+ the original double image (raw voxels) from which voxel measures are made
+
+ ROI : {dictionary}
+ Region of Interest structure that has blob bounding boxes
+
+
+ Returns
+ ..........
+
+ none
+
+ """
+ if ROI==None:
+ ROIList = NP.zeros(1, dtype=_objstruct)
+ [rows, cols] = label_image.shape
+ ROIList['L'] = 2
+ ROIList['R'] = cols-3
+ ROIList['B'] = 2
+ ROIList['T'] = rows-3
+
+ number_regions = ROI.size
+ indices = range(0, number_regions)
+ inflate = 1
+ for i in indices:
+ left = ROI[i]['L']
+ right = ROI[i]['R']
+ bottom = ROI[i]['B']
+ top = ROI[i]['T']
+ Label = ROI[i]['Label']
+ rows = top-bottom-1
+ cols = right-left-1
+ section= NP.zeros(rows*cols, dtype=raw_image.dtype).reshape(rows, cols)
+ section = raw_image[bottom:top, left:right] \
+ [label_image[bottom:top, left:right]==Label]
+ mask = section[section>0]
+ ROI[i]['voxelMean'] = mask.mean()
+ ROI[i]['voxelVar'] = mask.std()
+
+ return
+
+
def get_blob_regions(labeled_image, groups, dust=16):
"""
ROIList = get_blob_regions(labeled_image, groups, dust=16)
@@ -273,7 +414,7 @@
"""
- ROIList = N.zeros(groups, dtype=_objstruct)
+ ROIList = NP.zeros(groups, dtype=_objstruct)
# return the bounding box for each connected edge
S.get_blob_regions(labeled_image, ROIList)
@@ -305,7 +446,7 @@
"""
[rows, cols] = binary_edge_image.shape
- labeled_edge_image = N.zeros(rows*cols, dtype=N.uint16).reshape(rows, cols)
+ labeled_edge_image = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
groups = S.get_blobs(binary_edge_image, labeled_edge_image)
return labeled_edge_image, groups
@@ -339,7 +480,7 @@
"""
[rows, cols] = sobel_edge_image.shape
- sobel_edge = N.zeros(rows*cols, dtype=N.uint16).reshape(rows, cols)
+ sobel_edge = NP.zeros(rows*cols, dtype=NP.uint16).reshape(rows, cols)
S.sobel_edges(sobel_edge_image, sobel_edge, sobel_stats['ave_gt0'], sobel_stats['min_gt0'],
sobel_stats['max_gt0'], mode, sobel_threshold)
@@ -368,80 +509,125 @@
mean and nonzero min, max of sobel filtering
"""
+ filtered_slice = filtered_slice.astype(NP.float64)
[rows, cols] = filtered_slice.shape
- sobel_edge_image = N.zeros(rows*cols, dtype=N.float64).reshape(rows, cols)
+ sobel_edge_image = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
pAve, min_value, max_value = S.sobel_image(filtered_slice, sobel_edge_image)
- # replace this with numpy calls for the image stats. but C ext is faster
+ # can replace this with numpy calls for the image stats. but C ext is faster
# S.sobel_image(filtered_slice, sobel_edge_image)
- # pAve = sobel_edge_image[sobel_edge_image>0].mean()
- # min_value = sobel_edge_image[sobel_edge_image>0].min()
- # max_value = sobel_edge_image[sobel_edge_image>0].max()
+ # sobel_mask = sobel_edge_image[sobel_edge_image>0]
+ # pAve = sobel_mask.mean()
+ # min_value = sobel_mask.min()
+ # max_value = sobel_mask.max()
sobel_stats= {'ave_gt0' : pAve, 'min_gt0': min_value, 'max_gt0': max_value}
return sobel_edge_image, sobel_stats
-def pre_filter(slice, filter, low_threshold=2048+220, high_threshold=600+2048, conv_binary=0):
+def pre_filter(slice, filter, low_threshold=0, high_threshold=0, conv_binary=0):
"""
+ edge_filter = pre_filter(slice, filter, low_threshold=0, high_threshold=slice.max(), conv_binary=0)
+
take 2D image slice and filter and pre-filter and threshold prior to segmentation
+ Parameters
+ ..........
+ slice : {nd_array}
+ input 2D image. gets cast to int16
+
+ filter : {dictionary}
+ 2D filter kernel set from build_2d_kernel()
+
+ low_threshold : {int}, optional
+ default 0
+
+ high_threshold : {int}, optional
+ default max value of source image
+
+ conv_binary : {0, 1}, optional
+ flag to convert edge_filter image to binary valued. default
+ is binary conversion off
+
+ Returns
+ ..........
+ edge_filter : {nd_array}
+ filtered and thresholded image that can be (optional) binary.
+
"""
# make sure the input is 16 bits. this is input to edge machine
# so can handle raw and 8 bit scaled inputs
- slice = slice.astype(N.int16)
+ if high_threshold==0:
+ # default to the maximum value of the image
+ high_threshold = slice.max()
+
+ slice = slice.astype(NP.int16)
[rows, cols] = slice.shape
- edge_image = N.zeros(rows*cols, dtype=N.float64).reshape(rows, cols)
+ edge_image = NP.zeros(rows*cols, dtype=NP.float64).reshape(rows, cols)
S.edge_prefilter(low_threshold, high_threshold, filter['kernelSize'], filter['kernel'],
slice, edge_image)
if conv_binary == 1:
edge_image[edge_image>0] = 1
- edge_image = edge_image.astype(N.uint16)
+ edge_image = edge_image.astype(NP.uint16)
return edge_image
def get_max_bounding_box(ROI):
- max_area = ROI[:]['Area'].max()
- indices = range(0, ROI.size)
- for i in indices:
- if ROI[i]['Area'] == max_area:
- left = ROI[i]['L']
- right = ROI[i]['R']
- top = ROI[i]['T']
- bottom = ROI[i]['B']
+ """
+ bounding_box = get_max_bounding_box(ROI)
- bounding_box = {'left' : left, 'right' : right, 'top' : top, 'bottom' : bottom}
+ take an ROI structure and find the maximum area bounding box
+ Parameters
+ ..........
+
+ ROI : {dictionary}
+ the ROI is the automatically extracted blob regions of interest
+ and contains the rectangular bounding box of each blob.
+
+ Returns
+ ..........
+
+ bounding_box : {dictionary}
+ the Left, Right, Top and Bottom of the LARGEST bounding box in the ROI
+
+ """
+ max_index = ROI[:]['Area'].argmax()
+ bounding_box = {'left' : ROI[max_index]['L'], 'right' : ROI[max_index]['R'],
+ 'top' : ROI[max_index]['T'], 'bottom' : ROI[max_index]['B']}
+
return bounding_box
-def set_draw_bounding_box(bounding_box):
- x = N.zeros(5, dtype=N.uint16)
- y = N.zeros(5, dtype=N.uint16)
+def get_all_bounding_boxes(ROI):
+ """
+ measures = get_all_bounding_boxes(ROI)
- x[0] = bounding_box['left']
- x[1] = bounding_box['right']
- x[2] = bounding_box['right']
- x[3] = bounding_box['left']
- x[4] = bounding_box['left']
+ get all bounding boxes in the ROI (feature) dictionary
- y[0] = bounding_box['bottom']
- y[1] = bounding_box['bottom']
- y[2] = bounding_box['top']
- y[3] = bounding_box['top']
- y[4] = bounding_box['bottom']
+ Parameters
+ ..........
- return x, y
+ ROI : {dictionary}
+ the ROI is the automatically extracted blob regions of interest
+ and contains the rectangular bounding box of each blob.
-def get_all_bounding_boxes(ROI):
+ Returns
+ ..........
+
+ measures : {dictionary}
+ the Left, Right, Top and Bottom of all bounding boxes in the ROI
+
+ """
+
number = ROI.size
indices = range(0, ROI.size)
- _shortstruct = N.dtype([('left', 'i'),
+ _shortstruct = NP.dtype([('left', 'i'),
('right', 'i'),
('top', 'i'),
('bottom', 'i')])
- measures = N.zeros(number, dtype=_shortstruct)
+ measures = NP.zeros(number, dtype=_shortstruct)
for i in indices:
measures[i]['left'] = ROI[i]['L']
measures[i]['right'] = ROI[i]['R']
@@ -450,56 +636,34 @@
return measures
-def draw_all_bounding_boxes(measures):
- number = measures.size
- indices = range(0, measures.size)
- for i in indices:
- x, y = set_draw_bounding_box(measures[i])
- P.plot(x, y)
-def build_test_discs():
+def build_2d_kernel(aperature=21, hiFilterCutoff=10.0):
+ """
- radius = 50
- rows = 512
- cols = 512
- test_image = N.zeros(rows*cols, dtype=N.int16).reshape(rows, cols)
- y_indices = N.array(range(-radius, radius+1))
- center_x = rows / 4
- center_y = cols / 4
+ FIRFilter = build_2d_kernel(aperature, hiFilterCutoff)
- for i in y_indices:
- x = math.sqrt(float(radius)**2 - float(i)**2)
- test_image[1*center_y+i, 1*center_x-x:1*center_x+x] = 100
- test_image[1*center_y+i, 3*center_x-x:3*center_x+x] = 100
- test_image[3*center_y+i, 1*center_x-x:1*center_x+x] = 100
- test_image[3*center_y+i, 3*center_x-x:3*center_x+x] = 100
+ build hamming-windowed FIR filter with sinc kernel
- return test_image
+ Parameters
+ ..........
-def get_slice(imageName='slice112.raw', bytes=2, rows=512, columns=512):
- # clip the ends for this test CT image file as the spine runs off the end of the image
- ImageSlice = N.fromfile(imageName, dtype=N.uint16).reshape(rows, columns)
- ImageSlice[505:512, :] = 0
- return ImageSlice
+ aperature : {int}, optional
+ the number of coefficients in the filter. default is 21. needs to be ODD
-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
+ hiFilterCutoff : {float}
+ the upper cutoff in digital frequency units
- 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)
+ Returns
+ ..........
+ FIRFilter : {dictionary}
+ filter kernel
+
"""
rad = math.pi / 180.0
HalfFilterTaps = (aperature-1) / 2
- kernel = N.zeros((aperature), dtype=N.float64)
+ kernel = NP.zeros((aperature), dtype=NP.float64)
LC = 0.0
HC = hiFilterCutoff * rad
t2 = 2.0 * math.pi
@@ -530,17 +694,29 @@
def build_d_gauss_kernel(gWidth=20, 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 dictionary). Use in Canny filter call
+ build the derivative of Gaussian kernel for Canny edge filter
+
+ Parameters
+ ..........
+ gWdith : {int}, optional
+ width of derivative of Gaussian kernel.
+ default value is 20
+
+ sigma : {float}, optional
+ sigma term of derivative of Gaussian kernel
+ default value is 1.0
+
+ Returns
+ ..........
+
+ DGFilter : {dictionary}
+ filter kernel
+
"""
- kernel = N.zeros((1+gWidth), dtype=N.float64)
+ kernel = NP.zeros((1+gWidth), dtype=NP.float64)
indices = range(0, gWidth)
for i in indices:
@@ -554,23 +730,29 @@
def build_morpho_thin_masks():
"""
+ MATFilter = 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
+ Parameters
+ ..........
- Output:
- MATFilter (a struct)
+ None
+ Returns
+ ..........
+
+ MATFilter : {dictionary}
+ morphology filter kernels. there are 2 sets of 8 3x3 masks
+
"""
# (layers, rows, cols)
size = (8*3*3)
- J_mask = N.zeros(size, dtype=N.int16)
- K_mask = N.zeros(size, dtype=N.int16)
+ J_mask = NP.zeros(size, dtype=NP.int16)
+ K_mask = NP.zeros(size, dtype=NP.int16)
maskCols = 3
# load the 8 J masks for medial axis transformation
@@ -663,3 +845,42 @@
return MATFilter
+
+def build_laws_kernel():
+
+ """
+ LAWSFilter = build_laws_kernel()
+
+ build 6 length-7 Law's texture filter masks
+ mask names are: 'L', 'S', 'E', 'W', 'R', 'O'
+
+
+ Parameters
+ ..........
+
+ None
+
+ Returns
+ ..........
+
+ LAWSFilter : {dictionary}
+ a set of 6 length-7 Laws texture kernels
+
+ """
+ aperature = (6, 7)
+ coefficients = NP.zeros((aperature), dtype=NP.float64)
+ 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, 'filters' : 21,
+ 'coefficients': coefficients, 'names': names}
+
+ return LAWSFilter
+
+
More information about the Scipy-svn
mailing list