Automatic Change Detection Algorithm of Man-made Objects (ACDA)

For the automatic change detection of man-made objects, at first the Chan-Vese model will be used. Code for the model in Python found on GitHub and modified to detect man-made objects and compare them by using multitemporal data. This approach is widely applied in medical applications (example given bellow).

levelset.gif

The Chan-Vese Model

The Chan-Vese level sets model aims at minimizing the energy:

equation

where I is the image, C is the boundary of the segmented region, c-1 and c-2 are the averages of I respectively inside and outside C, and kappa is the curvature of C.

The next example highlights how the segmentation does not depend on image edges in an experiment similar to one of Chan and Vese’s examples. The method identifies the clusters
from their darker average intensity, even though the cluster boundaries are not well defined. The initial boundary is a circle and μ = 0.3.

1.png

Effect of Varying μ

While Chan–Vese has quite a few tuning parameters, the most important is μ. Parameter μ adjusts the length penalty, which balances between fitting the input image more accurately (smaller μ) vs. producing a smoother boundary (larger μ). In the example below, there are two groups of circles. Depending on μ, the circles are either segmented individually or as two clusters.

2.png

Effect of Varying ν

Parameter ν sets the penalty (or reward, if ν < 0) for area inside C. Note that this parameter is meaningful only when there is a prescribed inside vs. outside of the segmentation boundary. In the figure below, the evolution is shown over 850 iterations with five different values of ν where the initialization is a circle. When ν is too negative, the boundary expands to fill the full domain, and when ν is too positive, the boundary shrinks until it vanishes.

3.png

Results of the Automatic Change Detection algorithm

These two images are high resolution Quickbird data. The dimensions of the images is the same (about 200*200 pixels). It’s important the dimensions of the images to be exactly the same for the algorithm to work. In the example bellow, the initial boundary is a rectangle and after 500 iterations the algorithm detects the man-made objects with setting parameters μ=1 and v=0. The algorithm converts the segmentation boundaries into bin-images(black and white). By subtracting those two bin-images, the change appears as white (Third .GIF).

im6a-6b

The Change Detection Algorithm (ACDA) in Python

For this algorithm to work Python, Numpy, SciPy and Matplotlib must be installed to your computer.

# Buildings Detection Python
# Version 0.0.6
#********************************************************
#*Copyright 2017 Falagas A., Ouzoun M.
#*
#*Permission is hereby granted, free of charge,
#*to any person obtaining a copy of this software
#*and associated documentation files (the "Software"),
#*to deal in the Software without restriction,
#*including without limitationthe rights to use,
#*copy, modify, merge, publish, distribute, sublicense,
#*and/or sell copies of the Software, and to permit
#*persons to whom the Software is furnished to do so,
#*subject to the following conditions:
#*
#*The above copyright notice and this permission notice
#*shall be included in all copies or substantial
#*portions of the Software.
#*
#*THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY
#*OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT
#*LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#*FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
#*IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
#*BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
#*WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
#*ARISING FROM, OUT OF OR IN CONNECTION WITH THE
#*SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#********************************************************
#*Code for the Chan-Vese model created by:
#*Kevin Keraudren
#********************************************************

import numpy as np
import scipy.ndimage as nd
import matplotlib.pyplot as plt
import scipy.misc
import os
import sys
import warnings

eps = np.finfo(float).eps
warnings.simplefilter('ignore')
def chanvese(i, I, init_mask, max_its=200, alpha=0.2, thresh=0, color='r', display=False):
I = I.astype(np.float)

# Create a signed distance map (SDF) from mask
phi = mask2phi(init_mask)

if display:
plt.ion()
fig, axes = plt.subplots(ncols=1)
show_curve_and_phi(fig, I, phi, color)
plt.savefig('1_'+imnam+imform, bbox_inches='tight')
# Main loop
its = 0
stop = False
prev_mask = init_mask
c = 0
saveaf=1

while (its < max_its and not stop):
# Get the curve's narrow band
idx = np.flatnonzero(np.logical_and(phi <= 1.2, phi >= -1.2))
if len(idx) > 0:
# Intermediate output
if display:
if np.mod(its, 50) == 0:
print('iteration: {0}'.format(its))
show_curve_and_phi(fig, I, phi, color)
else:
if np.mod(its, 10) == 0:
print('iteration: {0}'.format(its))

# Find interior and exterior mean
upts = np.flatnonzero(phi <= 0)&nbsp; # interior points vpts = np.flatnonzero(phi > 0)&nbsp; # exterior points
u = np.sum(I.flat[upts]) / (len(upts) + eps)&nbsp; # interior mean
v = np.sum(I.flat[vpts]) / (len(vpts) + eps)&nbsp; # exterior mean

# Force from image information
F = (I.flat[idx] - u)**2 - (I.flat[idx] - v)**2
# Force from curvature penalty
curvature = get_curvature(phi, idx)

# Gradient descent to minimize energy
dphidt = F / np.max(np.abs(F)) + alpha * curvature

# Maintain the CFL condition
dt = 0.45 / (np.max(np.abs(dphidt)) + eps)

# Evolve the curve
phi.flat[idx] += dt * dphidt

# Keep SDF smooth
phi = sussman(phi, 0.5)

new_mask = phi <= 0
c = convergence(prev_mask, new_mask, thresh, c)
if its==(250*saveaf):
show_curve_and_phi(fig, I, phi, color)
plt.savefig(str(its)+'_'+imnam+imform, bbox_inches='tight')
saveaf=saveaf+1

if c <= 5:
its = its + 1
prev_mask = new_mask
else:
stop = True

else:
break

# Final output
if display:
show_curve_and_phi(fig, I, phi, color)
plt.savefig(str(its)+'_'+imnam+imform, bbox_inches='tight')

# Make mask from SDF
seg = phi <= 0&nbsp; # Get mask from levelset
scipy.misc.imsave(imnam+'_seg'+imform, seg)
fig.axes[0].cla()
fig.axes[0].imshow(seg, cmap='gray')
fig.axes[0].set_axis_off()
plt.draw()
plt.savefig('end_'+imnam+imform, bbox_inches='tight')

return seg, phi, its

def IBD(imi, imj):
imi=nd.imread(imi)
imj=nd.imread(imj)
imi=np.array(imi, dtype=int)
imj=np.array(imj, dtype=int)
s=np.abs(imi-imj)
x, y = s.shape
k=0
l=0
while k<x:
while l<y:
if s[k,l]<=255 and s[k,l]>10:
s[k,l]=1
else:
s[k,l]=0
l=l+1
k=k+1

scipy.misc.imsave(imnam+'_det'+imform,s)
return s

# -----------------------------------------------------------------------

def bwdist(a):
"""
Intermediary function. 'a' has only True/False vals,
so we convert them into 0/1 values - in reverse.
True is 0, False is 1, distance_transform_edt wants it that way.
"""
return nd.distance_transform_edt(a == 0)

# Displays the image with curve superimposed
def show_curve_and_phi(fig, I, phi, color):
fig.axes[0].cla()
fig.axes[0].imshow(I, cmap='gray')
fig.axes[0].contour(phi, 0, colors=color)
fig.axes[0].set_axis_off()
plt.draw()
plt.pause(0.0001)

def im2double(a):
a = a.astype(np.float)
a /= np.abs(a).max()
return a

# Converts a mask to a SDF
def mask2phi(init_a):
phi = bwdist(init_a) - bwdist(1 - init_a) + im2double(init_a) - 0.5
return phi

# Compute curvature along SDF
def get_curvature(phi, idx):
dimy, dimx = phi.shape
yx = np.array([np.unravel_index(i, phi.shape) for i in idx])&nbsp; # subscripts
y = yx[:, 0]
x = yx[:, 1]

# Get subscripts of neighbors
ym1 = y - 1
xm1 = x - 1
yp1 = y + 1
xp1 = x + 1

# Bounds checking
ym1[ym1 < 0] = 0
xm1[xm1 < 0] = 0 yp1[yp1 >= dimy] = dimy - 1
xp1[xp1 >= dimx] = dimx - 1

# Get indexes for 8 neighbors
idup = np.ravel_multi_index((yp1, x), phi.shape)
iddn = np.ravel_multi_index((ym1, x), phi.shape)
idlt = np.ravel_multi_index((y, xm1), phi.shape)
idrt = np.ravel_multi_index((y, xp1), phi.shape)
idul = np.ravel_multi_index((yp1, xm1), phi.shape)
idur = np.ravel_multi_index((yp1, xp1), phi.shape)
iddl = np.ravel_multi_index((ym1, xm1), phi.shape)
iddr = np.ravel_multi_index((ym1, xp1), phi.shape)

# Get central derivatives of SDF at x,y
phi_x = -phi.flat[idlt] + phi.flat[idrt]
phi_y = -phi.flat[iddn] + phi.flat[idup]
phi_xx = phi.flat[idlt] - 2 * phi.flat[idx] + phi.flat[idrt]
phi_yy = phi.flat[iddn] - 2 * phi.flat[idx] + phi.flat[idup]
phi_xy = 0.25 * (- phi.flat[iddl] - phi.flat[idur] +
phi.flat[iddr] + phi.flat[idul])
phi_x2 = phi_x**2
phi_y2 = phi_y**2

# Compute curvature (Kappa)
curvature = ((phi_x2 * phi_yy + phi_y2 * phi_xx - 2 * phi_x * phi_y * phi_xy) /
(phi_x2 + phi_y2 + eps) ** 1.5) * (phi_x2 + phi_y2) ** 0.5

return curvature

# Level set re-initialization by the sussman method
def sussman(D, dt):
# forward/backward differences
a = D - np.roll(D, 1, axis=1)
b = np.roll(D, -1, axis=1) - D
c = D - np.roll(D, -1, axis=0)
d = np.roll(D, 1, axis=0) - D

a_p = np.clip(a, 0, np.inf)
a_n = np.clip(a, -np.inf, 0)
b_p = np.clip(b, 0, np.inf)
b_n = np.clip(b, -np.inf, 0)
c_p = np.clip(c, 0, np.inf)
c_n = np.clip(c, -np.inf, 0)
d_p = np.clip(d, 0, np.inf)
d_n = np.clip(d, -np.inf, 0)

a_p[a < 0] = 0 a_n[a > 0] = 0
b_p[b < 0] = 0 b_n[b > 0] = 0
c_p[c < 0] = 0 c_n[c > 0] = 0
d_p[d < 0] = 0 d_n[d > 0] = 0

dD = np.zeros_like(D)
D_neg_ind = np.flatnonzero(D < 0) D_pos_ind = np.flatnonzero(D > 0)

dD.flat[D_pos_ind] = np.sqrt(
np.max(np.concatenate(
([a_p.flat[D_pos_ind]**2], [b_n.flat[D_pos_ind]**2])), axis=0) +
np.max(np.concatenate(
([c_p.flat[D_pos_ind]**2], [d_n.flat[D_pos_ind]**2])), axis=0)) - 1
dD.flat[D_neg_ind] = np.sqrt(
np.max(np.concatenate(
([a_n.flat[D_neg_ind]**2], [b_p.flat[D_neg_ind]**2])), axis=0) +
np.max(np.concatenate(
([c_n.flat[D_neg_ind]**2], [d_p.flat[D_neg_ind]**2])), axis=0)) - 1

D = D - dt * sussman_sign(D) * dD
return D

def sussman_sign(D):
return D / np.sqrt(D**2 + 1)

# Convergence Test
def convergence(p_mask, n_mask, thresh, c):
diff = p_mask - n_mask
n_diff = np.sum(np.abs(diff))
if n_diff < thresh:
c = c + 1
else:
c = 0
return c

if __name__ == "__main__":
i=1
imn=2
cwd = os.getcwd()
print 'Current Working Directory', cwd
files = 0
for file in os.listdir(cwd):
print file
files = files+1
imform=str(raw_input("Press the image's format (ex .TIFF, .png, .JPEG):"))
f = 0
for file in os.listdir(cwd):
if file.endswith(imform):
print file
f = f+1
if f == 0:
print 'Could not find {} files'.format(imform)
sys.exit(1)
while i <= imn:
img = str(raw_input('Press image name:'))
imnam, imformat=img.split('.')
img = nd.imread(img, flatten=True)
mask = np.zeros(img.shape)
mask[10:100, 10:100] = 1
chanvese(i, img, mask, max_its=1500, display=True, alpha=1)
if i==1:
imi=imnam+'_seg'+imform
else:
imj=imnam+'_seg'+imform
i=i+1
det=IBD(imi, imj)

Special thanks to my friend and classmate M. Ouzoun.

References

  1. T. Chan and L. Vese, Active contours without edges, IEEE transactions on image processing 10(2) (2001), pp. 266-
    277
  2. Robert Crandall, Image Segmentation Using the Chan-Vese Algorithm, ECE 532 Project Fall (2009)
  3. Konstantinos Karantzalos and Demetre Argialas, A Region-based Level Set Segmentation for Automatic Detection of Man-made Objects from Aerial and Satellite Images, Photogrammetric Engineering & Remote Sensing Vol. 75, No. 6, June 2009, pp. 667–677
  4. Pascal Getreuer, Chan–Vese Segmentation, Image Processing On Line on 2012–08–08
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s