import numpy as npOnce the code is finished, the last u will have the final values. More detailed code along with the sparse matrix version is here.

import scipy.linalg

# Number of internal points

N = 200

# Calculate Spatial Step-Size

h = 1/(N+1.0)

k = h/2

x = np.linspace(0,1,N+2)

x = x[1:-1] # get rid of the '0' and '1' at each end

# Initial Conditions

u = np.transpose(np.mat(10*np.sin(np.pi*x)))

# second derivative matrix

I2 = -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 = 1

NumOfTimeSteps = 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)

## 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:

## 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

Code

## Saturday, January 15, 2011

### Discrete Laplace Transform - del2

I ported Matlab del2() call into Python. Here it is for reference:

import numpy as np

def 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

Code

## Friday, January 14, 2011

### Gaussian Convolution, Blurring using Python

import matplotlib.pyplot as pltResults:

import numpy as np

import scipy.signal as signal

def 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()

Before Image

After Image

Subscribe to:
Posts (Atom)