[Scipy-svn] r3506 - in trunk/scipy/ndimage: segment tests

scipy-svn@scip... scipy-svn@scip...
Wed Nov 7 14:36:41 CST 2007


Author: oliphant
Date: 2007-11-07 14:36:40 -0600 (Wed, 07 Nov 2007)
New Revision: 3506

Added:
   trunk/scipy/ndimage/segment/objectdata.py
   trunk/scipy/ndimage/tests/test_segment.py
Log:
Add missing files to ndimage.segment

Added: trunk/scipy/ndimage/segment/objectdata.py
===================================================================
--- trunk/scipy/ndimage/segment/objectdata.py	2007-11-07 18:45:41 UTC (rev 3505)
+++ trunk/scipy/ndimage/segment/objectdata.py	2007-11-07 20:36:40 UTC (rev 3506)
@@ -0,0 +1,25 @@
+
+import numpy as N
+
+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)]
+                   )
+

Added: trunk/scipy/ndimage/tests/test_segment.py
===================================================================
--- trunk/scipy/ndimage/tests/test_segment.py	2007-11-07 18:45:41 UTC (rev 3505)
+++ trunk/scipy/ndimage/tests/test_segment.py	2007-11-07 20:36:40 UTC (rev 3506)
@@ -0,0 +1,142 @@
+
+import numpy as N
+from numpy.testing import *
+import scipy.ndimage.segment as S
+
+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, 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):
+	# 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):
+	# 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):
+	# 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):
+	#
+	# 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):
+	#
+	# pass raw image, labeled mask and the partially filled ROIList
+	# VoxelMeasures will fill the texture (Law's, co-occurence, Gabor) features in the list
+	#
+	S.texture_measures(rawImage, labeledEdges, ROIList)
+	return 
+
+def segment_regions():
+	# 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():
+	# 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):
+	# 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='junk.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='junk.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)
+
+
+class TestSegment(NumpyTestCase):
+    def test1(self):
+	image = get_slice(filename)
+	sourceImage = image.copy()
+	edges, objects = sobel(image)
+	get_shape_mask(edges, objects)
+	get_voxel_measures(sourceImage, edges, objects)
+	get_texture_measures(sourceImage, edges, objects)
+
+    def test2(self):
+	sourceImage, labeledMask, ROIList = segment_regions()
+
+    def test3(self):
+	regionMask, numberRegions = grow_regions()
+	print regionMask.max()
+	#save_slice(regionMask, 'regionMask.raw')
+
+    
+if __name__ == "__main__":
+    NumpyTest().run()



More information about the Scipy-svn mailing list