## Sunday, January 30, 2011

### Heat Equation solution using Finite Difference and Crank-Nicholson

I rewrote parts of this code so that it used dense (non-sparse) matrices instead of sparse matrices, for demonstration purposes. The code is below:
`import numpy as npimport scipy.linalg# Number of internal pointsN = 200# Calculate Spatial Step-Sizeh = 1/(N+1.0)k = h/2x = np.linspace(0,1,N+2)x = x[1:-1] # get rid of the '0' and '1' at each end# Initial Conditionsu = np.transpose(np.mat(10*np.sin(np.pi*x)))# second derivative matrixI2 = -2*np.eye(N)E = np.diag(np.ones((N-1)), k=1)D2 = (I2 + E + E.T)/(h**2)I = np.eye(N)TFinal = 1NumOfTimeSteps = int(TFinal/k)for i in range(NumOfTimeSteps):   # Solve the System: (I - k/2*D2) u_new = (I + k/2*D2)*u_old   A = (I - k/2*D2)   b = np.dot((I + k/2*D2), u)   u = scipy.linalg.solve(A, b)`
Once the code is finished, the last u will have the final values. More detailed code along with the sparse matrix version is here.

## Thursday, January 20, 2011

### Chan Vese Segmentation Algorithm

Attached is the Chan-Vese segmentation algorithm written in C++ by Robert Crandall. Just run make then run the binary "segment" which will process the test.pgm file and create CVoutEdges.pgm and CVoutRegions.pgm files. In the attached zip you can also find his project paper file ECE532_ProjectPaper.pdf.

Code

## Saturday, January 15, 2011

### Discrete Laplace Transform - del2

I ported Matlab del2() call into Python. Here it is for reference:
`import numpy as npdef del2(M):  dx = 1  dy = 1  rows, cols = M.shape  dx = dx * np.ones ((1, cols - 1))  dy = dy * np.ones ((rows-1, 1))  mr, mc = M.shape  D = np.zeros ((mr, mc))  if (mr >= 3):      ## x direction      ## left and right boundary      D[:, 0] = (M[:, 0] - 2 * M[:, 1] + M[:, 2]) / (dx[:,0] * dx[:,1])      D[:, mc-1] = (M[:, mc - 3] - 2 * M[:, mc - 2] + M[:, mc-1]) \          / (dx[:,mc - 3] * dx[:,mc - 2])         ## interior points      tmp1 = D[:, 1:mc - 1]      tmp2 = (M[:, 2:mc] - 2 * M[:, 1:mc - 1] + M[:, 0:mc - 2])      tmp3 = np.kron (dx[:,0:mc -2] * dx[:,1:mc - 1], np.ones ((mr, 1)))      D[:, 1:mc - 1] = tmp1 + tmp2 / tmp3  if (mr >= 3):      ## y direction      ## top and bottom boundary      D[0, :] = D[0,:]  + \          (M[0, :] - 2 * M[1, :] + M[2, :] ) / (dy[0,:] * dy[1,:])         D[mr-1, :] = D[mr-1, :] \          + (M[mr-3,:] - 2 * M[mr-2, :] + M[mr-1, :]) \          / (dy[mr-3,:] * dx[:,mr-2])         ## interior points      tmp1 = D[1:mr-1, :]      tmp2 = (M[2:mr, :] - 2 * M[1:mr - 1, :] + M[0:mr-2, :])      tmp3 = np.kron (dy[0:mr-2,:] * dy[1:mr-1,:], np.ones ((1, mc)))      D[1:mr-1, :] = tmp1 + tmp2 / tmp3     return D / 4`

### Level Sets and Image Segmentation with Python

I ported the Matlab code that was written by Li, Xu, Gui and Fox for their paper Level Set Evolution Without Re-initialization: A New Variational Formulation, into Python. Li et al. use level sets with a new approach and are able to segment an image succesfully. Some screenshots as well as the code itself is below.

Code

## Friday, January 14, 2011

### Gaussian Convolution, Blurring using Python

`import matplotlib.pyplot as pltimport numpy as npimport scipy.signal as signaldef gauss_kern():   """ Returns a normalized 2D gauss kernel array for convolutions """   h1 = 15   h2 = 15   x, y = np.mgrid[0:h2, 0:h1]   x = x-h2/2   y = y-h1/2   sigma = 1.5   g = np.exp( -( x**2 + y**2 ) / (2*sigma**2) );   return g / g.sum()Img = plt.imread("twoObj.bmp")Img = Img[::-1]g = gauss_kern()Img_smooth = signal.convolve(Img,g,mode='same')plt.imshow(Img_smooth)plt.show()`
Results:

Before Image

After Image