[Scipy-svn] r3547 - trunk/scipy/ndimage
scipy-svn@scip...
scipy-svn@scip...
Fri Nov 16 19:15:10 CST 2007
Author: tom.waite
Date: 2007-11-16 19:15:06 -0600 (Fri, 16 Nov 2007)
New Revision: 3547
Modified:
trunk/scipy/ndimage/segmenter.py
Log:
Move C code to Python. Start with setup functions.
Modified: trunk/scipy/ndimage/segmenter.py
===================================================================
--- trunk/scipy/ndimage/segmenter.py 2007-11-16 12:50:30 UTC (rev 3546)
+++ trunk/scipy/ndimage/segmenter.py 2007-11-17 01:15:06 UTC (rev 3547)
@@ -1,4 +1,4 @@
-
+import math
import numpy as N
import scipy.ndimage.segment as S
@@ -282,6 +282,7 @@
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
@@ -290,4 +291,203 @@
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
+
More information about the Scipy-svn
mailing list