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 np
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)

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

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

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

Gaussian Convolution, Blurring using Python

import matplotlib.pyplot as plt
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()

Results:

Before Image



After Image