UNIST 20115009 Kyunghoon Kim

Recovering Line-Networks in Images by Junction-Point Processes

http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6619091&tag=1

This is about simulation of CVPR paper.



In [1]:

    
from PIL import Image # Load PIL(Python Image Library)
import matplotlib.pyplot as plt # Load pyplot library
import numpy as np # Load Numerical Library
import copy
import timeit
import math
import random



In [2]:

    
im = Image.open('road.jpg') # Load the image file
imarray = np.array(im) # Construct the array
print im.size



In [3]:

    
imgray = im.convert('L') # Convert the image to gray scale
plt.imshow(imgray, cmap=plt.get_cmap('gray'), interpolation='nearest') # interpolation='nearest' simply display the image without try to interpolate betwen pixels if the display resolution is not the same as the image resolution (which is most often the case). It will results in an image in which is image pixel is displayed as a square of multiple display pixels.

imgray_array = np.array(imgray) # Let imgray_array be the array of imgray.
print imgray_array # As we know, each pixels have a color value.









    



[[ 13  19  20 ...,  38  37  27]
 [ 13  19  21 ...,  36  37  27]
 [ 17  29  18 ...,  35  39  27]
 ..., 
 [151 241 230 ...,  97 106  75]
 [151 239 230 ..., 106 106  67]
 [150 234 229 ..., 113 117  77]]



In [4]:

    
plt.hist(imgray_array.flatten(), 256, [0,256]);



In [5]:

    
img = imgray_array
print im.size

Poisson Distribution $f(n;\lambda)=\frac{\lambda^n e^{-\lambda}}{n!}$

The Poisson parameter $\lambda_F$ which represents the expected number of junction-points is calculated as the average number of junction-points in the annotated image samples.



In [8]:

    
lambda_F = 10
def poisson_dist(n, lamb):
    return (lamb**n*math.exp(-lamb))/math.factorial(n)

the Occurrence probability of a k-junction-point in $x$

The Maximum Likelihood Estimate of each occurrence probability $p_k \in \mathcal{P}$ is given by the rate of k-junction-points in the annotated image samples.



In [8]:

$\alpha$ : the Collection of the parameter sets $\{\alpha_1 , \cdots , \alpha_k \}$ of the Dirichlet distribution of the k-junction-points with $k\geq 1$.

Since there is no closed-form solution for the Maximum Likelihood Estimate of Dirichlet distribution, we apply a fixed-point iteration to estimate the parameter set $\alpha$[1].

[1]: Minka, Thomas P. "Estimating a Dirichlet distribution.", 2000



In [8]:

$H_w$

$H_w$ is estimated by counting the number of lines with different widths from the annotated image samples.



In [8]:

Optimization

RJMCMC

Birth and death kernel

$$\frac{Q_{BD} (x' \rightarrow x )}{Q_{BD} (x \rightarrow x' )} = \frac{p_d}{p_b}\frac{\lambda_F}{n(x')}$$

$\lambda_F$ is the Poisson parameter representing the expected number of junction-points in the domain $F$.

$n(x')$ is the number of jucntion-points in the proposed configuration $x'$.

$p_d$ is the probability of choosing a death.

$p_b$ is the probability of choosing a birth. In practice, $p_d=p_b=0.5$.

Transition kernel

The translation kernel allows a junction-point to be moved without modifying the graph complexity. As the translation of a junction-point is proposed randomly, the proposition kernel ratio is simply given by $$\frac{Q_T (x'\rightarrow x)}{Q_T(x\rightarrow x')}=1$$

TEST



In [9]:

    
plt.imshow(img, cmap=plt.get_cmap('gray'), interpolation='nearest')









    Out[9]:





<matplotlib.image.AxesImage at 0xb34fdd8>



In [10]:

    
imgy, imgx = img.shape
print imgx, imgy
total = imgx*imgy

Image gradient Magnitude



In [11]:

    
from scipy.ndimage import filters 
imx = np.zeros(img.shape)
imy = np.zeros(img.shape)
#filters.sobel(img, 1, imx)

$I_x$



In [12]:

    
for j in xrange(imgy-1):
    for i in xrange(imgx-1):
        imx[j,i]=img[j,i+1]-img[j,i] # delta x is 1









    



-c:3: RuntimeWarning: overflow encountered in ubyte_scalars



In [13]:

    
plt.imshow(imx, cmap=plt.get_cmap('gray'), interpolation='nearest')









    Out[13]:





<matplotlib.image.AxesImage at 0xa048cc0>

$I_y$



In [14]:

    
for j in xrange(imgx-1):
    for i in xrange(imgy-1):
        imy[i, j]=img[i+1,j]-img[i, j] # delta y is 1



In [15]:

    
plt.imshow(imy, cmap=plt.get_cmap('gray'), interpolation='nearest')









    Out[15]:





<matplotlib.image.AxesImage at 0xa61f940>

Gradient Magnitude $|\nabla I|=\sqrt{I_x^2 + I_y^2}$



In [16]:

    
img_grad = np.zeros(img.shape)
for j in xrange(imgy-1):
    for i in xrange(imgx-1):
        img_grad[j, i] = math.sqrt(imx[j, i]**2 + imy[j, i]**2)



In [17]:

    
plt.imshow(img_grad, cmap=plt.get_cmap('gray'), interpolation='nearest')









    Out[17]:





<matplotlib.image.AxesImage at 0xa64a5c0>



In [18]:

    
fig, ax = plt.subplots(1, 4, figsize=(15,10))
ax[0].imshow(img, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[0].set_title('original image')
ax[1].imshow(imx, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[1].set_title('x-derivative')
ax[2].imshow(imy, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[2].set_title('y-derivative')
ax[3].imshow(img_grad, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[3].set_title('Gradient magnitude')









    Out[18]:





<matplotlib.text.Text at 0xb98d4e0>

But we use the sobel filter for the performance. http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.filters.sobel.html



In [19]:

    
from scipy.ndimage import filters

imx = np.zeros(img.shape)
filters.sobel(img, 1, imx)
imy = np.zeros(img.shape)
filters.sobel(img, 0, imy)
img_grad = np.sqrt(imx**2 + imy**2)

fig, ax = plt.subplots(1, 4, figsize=(15,10))
ax[0].imshow(img, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[0].set_title('original image')
ax[1].imshow(imx, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[1].set_title('x-derivative')
ax[2].imshow(imy, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[2].set_title('y-derivative')
ax[3].imshow(img_grad, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[3].set_title('Gradient magnitude')









    Out[19]:





<matplotlib.text.Text at 0xc7a7a20>



In [20]:

    
plt.imshow(img_grad, cmap=plt.get_cmap('gray'), interpolation='nearest')









    Out[20]:





<matplotlib.image.AxesImage at 0xc93b908>



In [21]:

    
plt.plot(img_grad[0]) # for first row of image









    Out[21]:





[<matplotlib.lines.Line2D at 0xd078c50>]

So I assume that the boundary $\partial G_x$ of the line region if img_grad value is bigger than 200

Random Choice



In [22]:

    
def rand(imgx, imgy):
    x = random.choice(range(imgx))
    y = random.choice(range(imgy))
    return x, y

========================================================================================

Initialize $X_0 = w^{(0)}$, $t=0$



In [340]:

    
sample = 10 # the number of randomly selected points

fig, ax = plt.subplots(1,1, figsize=(5,5))
ax.set_xlim([0,300])
# ax.set_ylim([0,200])

ax.imshow(img, cmap=plt.get_cmap('gray'), interpolation='nearest')
w = []
for i in xrange(sample):
    x, y = rand(imgx, imgy)
    w.append([x,y])
    ax.plot(x, y, '.r')

A. Chose a sub-kernel $Q_m$ according to probability $p_m$

$$Q(w\rightarrow .) = \sum_m p_m Q_m(w\rightarrow . )$$

$p_1=0.5$, $Q_{BD}$, Birth and death kernel
$p_2=0.5$, $Q_{T}$, Transition kernel



In [177]:

    
r = random.random()
if r>0.5:
    Q = Q_BD
else:
    Q = Q_T

A.1 Birth and death kernel

$$\frac{Q_{BD} (x' \rightarrow x )}{Q_{BD} (x \rightarrow x' )} = \frac{p_d}{p_b}\frac{\lambda_F}{n(x')}$$

$\lambda_F$ is the Poisson parameter representing the expected number of junction-points in the domain $F$.

$n(x')$ is the number of jucntion-points in the proposed configuration $x'$.

$p_d$ is the probability of choosing a death.

$p_b$ is the probability of choosing a birth. In practice, $p_d=p_b=0.5$.



In [417]:

    
def Q_BD(w):
    if random.random()>0.5:
        # Birth
        x, y = rand(imgx, imgy)
        w.append([x, y])
    else:
        # Death
        w.pop(w.index(random.choice(w)))
    return w

A.2 Transition kernel



In [416]:

    
def Q_T(w):
    return 1

In this kernel, Green ratio is just calculated by the density $h(.)$ and the temparature $T_t$.

$\lambda_F$, lambda_F



In [418]:

    
n = 11 # the number of junction-points
p = float(n)/total
lambda_F = n*p
print lambda_F









    



0.0018501529052

B. the probability density $h$

$h(x)\propto h_d(x)h_p(x)$

B.1 Data consistency $h_d$

We use the two conventional assumptions(similar colors, large gradients) to design the data consistency density $h_d$. $$h_d(x)\propto \prod_{p\in I}P(I_p | G_x )$$



In [311]:

    
plt.figure(figsize=(15,15))
plt.imshow(img_grad, interpolation='nearest')









    Out[311]:





<matplotlib.image.AxesImage at 0x33d93c88>

We must make 'Region-based intensity' and 'Boundary-based intensity'. But I don't know the algorithm. So just consider the simple approaches.

Boundary means highvalue-lowvalue-highvalue. So I will extract this pattern from each row.



In [46]:

    
plt.imshow(img_grad[80:100,60:80], interpolation='nearest')









    Out[46]:





<matplotlib.image.AxesImage at 0xf82b6a0>



In [49]:

    
plt.plot(img_grad[80:100,60:80]);

E.g., this pattern is boundary. I consider that the region is area between boundarys. So the region [13-15] is one region.



In [52]:

    
plt.plot(img_grad[80:100,60:80][0]);



In [64]:

    
plt.plot(img_grad[80])









    Out[64]:





[<matplotlib.lines.Line2D at 0x11c12208>]



In [342]:

    
row = img_grad[80]
left = row[0]
temp = []
for right in row[1:]:
    if right>left and right>200: # differential > 0 and threshold of intensity
        left = right
        temp.append(right)
    else:
        left = right
        temp.append(0)
plt.plot(temp)









    Out[342]:





[<matplotlib.lines.Line2D at 0x331217b8>]



In [344]:

    
interval = 8 # Take the interval [pixel, pixel+interval] at each pixel

index = 0
region = np.zeros(img.shape)
for value in temp:
    next_index = min(index+interval, len(temp))
    test_box = []
    if value != 0:
        d1 = temp[index]
        for d2 in temp[index+1:next_index]:
            if d2 > d1:
                test_box.append(1)
            else:
                test_box.append(0)
            d1 = d2

    """ Search pattern """
    pattern_box = array(test_box)
    pos = np.where(pattern_box == 1)[0]
    if pos.size > 0:
        d1 = pos[0]
        ox = 0
        for d2 in pos[1:]:
            if d2 != d1+1:
                ox = 1
            d1 = d2
        if ox == 1:
            print index, test_box
            for k in range(min(pos),max(pos)+1):
                region[0, index+k] = 1
    index += 1









    



71 [1, 0, 1, 0, 0, 0, 0]
280 [1, 0, 1, 0, 0, 0, 0]

B.1.1 Region-based intensity Function

I consider that the region is area between boundarys.



In [349]:

    
def region_based(interval):
    threshold = 300
    region = np.zeros(img.shape)
    column = 0
    for row in img_grad:
        left = row[0]
        temp = []
        for right in row[1:]:
            if right>left and right>threshold: # differential > 0
                left = right
                temp.append(right)
            else:
                left = right
                temp.append(0)

        index = 0
        for value in temp:
            next_index = min(index+interval, len(temp))
            test_box = []
            if value != 0:
                d1 = temp[index]
                for d2 in temp[index+1:next_index]:
                    if d2 > d1:
                        test_box.append(1)
                    else:
                        test_box.append(0)
                    d1 = d2

            """ Search pattern """
            pattern_box = array(test_box)
            pos = np.where(pattern_box == 1)[0]
            if pos.size > 0:
                d1 = pos[0]
                ox = 0
                for d2 in pos[1:]:
                    if d2 != d1+1:
                        ox = 1
                    d1 = d2
                if ox == 1:
                    for k in range(min(pos),max(pos)+1):
                        region[column, index+k] = 1

            index += 1
        column += 1
    return region



In [350]:

    
fig, ax = plt.subplots(1, 4, figsize=(15,10))
ax[0].imshow(region_based(4), cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[0].set_title('Interval 4')
ax[1].imshow(region_based(5), cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[1].set_title('Interval 5')
ax[2].imshow(region_based(8), cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[2].set_title('Interval 8')
ax[3].imshow(region_based(10), cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[3].set_title('Interval 10')









    Out[350]:





<matplotlib.text.Text at 0x3669d080>

If we combine with "interval 4" and "interval 8" for noise remove or consider vertical direction, then more good. But there is no time. We use "interval 8" as the region-based intensity.

Compare it with original image.



In [347]:

    
plt.imshow(img, cmap=plt.get_cmap('gray'), interpolation='nearest')









    Out[347]:





<matplotlib.image.AxesImage at 0x33a13a90>

B.1.2 Boundary-based intensity Function



In [348]:

    
fig, ax = plt.subplots(1, 4, figsize=(15,10))
ax[0].imshow(img_grad>200, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[0].set_title('|gradient| > 200')
ax[1].imshow(img_grad>300, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[1].set_title('|gradient| > 300')
ax[2].imshow(img_grad>400, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[2].set_title('|gradient| > 400')
ax[3].imshow(img_grad>500, cmap=plt.get_cmap('gray'), interpolation='nearest')
ax[3].set_title('|gradient| > 500')









    Out[348]:





<matplotlib.text.Text at 0x31de2470>

B.1.3 the normalized histogram $H_r$

The histograms $H_r$ is estimated by counting the number of pixels with different colors in line regions.



In [286]:

    
line_region = region_based(8)



In [295]:

    
histogram1 = np.zeros(img.shape)
index_row = 0
for row in line_region:
    index_column = 0
    for column in row:
        if line_region[index_row, index_column] > 0:
            histogram1[index_row, index_column] = img[index_row, index_column]
        index_column += 1
    index_row += 1



In [296]:

    
plt.imshow(histogram1, interpolation='nearest')









    Out[296]:





<matplotlib.image.AxesImage at 0x33d9e438>

line_region_dist is $H_r$



In [307]:

    
flat = histogram1.flatten()
flat = np.delete(flat, np.where(flat == 0)) # delete 0 for distribution
line_region_dist = plt.hist(flat, 256, normed=True);

B.1.4 the normalized histogram $H_n$

The histograms $H_n$ is estimated by counting the number of pixels with different colors in non-line regions.



In [298]:

    
histogram2 = np.zeros(img.shape)
index_row = 0
for row in line_region:
    index_column = 0
    for column in row:
        if line_region[index_row, index_column] == 0:
            histogram2[index_row, index_column] = img[index_row, index_column]
        index_column += 1
    index_row += 1



In [352]:

    
plt.imshow(histogram2, interpolation='nearest')









    Out[352]:





<matplotlib.image.AxesImage at 0x36c85940>

non-line_region_dist is $H_n$



In [308]:

    
flat = histogram2.flatten()
flat = np.delete(flat, np.where(flat == 0)) # delete 0 for distribution
nonline_region_dist = plt.hist(flat, 256, normed=True);

B.1.5 local likelihood $P(I_p | G_x )$ on each pixel p of the image

$I_p$ is the color intensity of the image at $p$.

The local likelihood at pixel $p$ is expressed by taking into account both color similarity on the lines and discontinuity on the line borders.



In [383]:

    
p=[0,0] # pixel
Ip = img[p[0], p[1]] # i_p , intensity of pixel p



In [384]:

    
G_x = w # Modify



In [385]:

    
Ip









    Out[385]:





13



In [387]:

    
if line_region[p[0], p[1]] > 0:
    nabla = img_grad[p]
else:
    nabla = 1
    
if p in G_x:
    """d_inside"""
    print 'inside'
    dvalue = line_region_dist[0][Ip]*nabla
else:
    print 'outside'
    """d_outside"""
    dvalue = nonline_region_dist[0][Ip]*nabla
print "P(I_p|G_x)=", dvalue









    



outside
P(I_p|G_x)= 0.00429698975466



In [404]:

    
def data_consistency(w, G):
    result = []
    for wi in w:
        i, j = wi # index
        Ip = img[j, i] # intensity
        
        if line_region[j, i] > 0:
            nabla = img_grad[j, i]
        else:
            nabla = 1

        if p in G_x:
            """d_inside"""
            print 'inside'
            dvalue = line_region_dist[0][Ip]*nabla
        else:
            print 'outside'
            """d_outside"""
            dvalue = nonline_region_dist[0][Ip]*nabla
        result.append(dvalue)
        print "P(I_p|G_x)=", dvalue
    return result



In [407]:

    
likelihood = data_consistency(w, G_x)









    



outside
P(I_p|G_x)= 0.0139859550702
outside
P(I_p|G_x)= 0.00290337145585
outside
P(I_p|G_x)= 0.00218997161242
outside
P(I_p|G_x)= 0.0180340844144
outside
P(I_p|G_x)= 0.00504357098617
outside
P(I_p|G_x)= 0.0
outside
P(I_p|G_x)= 0.0
outside
P(I_p|G_x)= 0.0238740087141
outside
P(I_p|G_x)= 0.00567401735944
outside
P(I_p|G_x)= 0.00285359937375

So Data consistency is



In [414]:

    
h_d = 1
for i in likelihood:
    h_d = h_d*i
    print h_d
print '\nh_d = ', h_d









    



0.0139859550702
4.06064227337e-05
8.89269130684e-08
1.60371545698e-09
8.08845274892e-12
0.0
0.0
0.0
0.0
0.0

h_d =  0.0

B.2 Shape priors $h_p$

Three different criteria are taken into account to characterize the shape of a graph $G_x$ from a junction-point configuration $x$:

$h_p \propto h_{connectivity}(x)h_{orientation}(x)h_{width}(x)$

B.2.1 Graph connectivity

$h_{connectivity} (x) = \frac{(\sum_{k\geq 1}n_k)!}{\prod_{k\geq1}n_k !}\prod_{k\geq1}p_k^{n_k}$

B.2.2 Edge orientation

This use Dirichlet law.

C. RJMCMC(Reversible Jump Marcov Chain Monte Carlo)

Just iterate the Algorithm.

In the top page, we already let the $X_0$.

So we iterate,

Choose a sub-kernel $Q_m$ according to probability $p_m$.
Perturb $w$ to $w'$ according to $Q_m(w\rightarrow .)
Compute the Green ratio

$$R=\frac{Q_{m} (w' \rightarrow w )}{Q_{m} (w \rightarrow w' )} (\frac{h(w')}{h(w)})^{\frac{1}{T_t}} $$

Choose $X_{t+1}=w'$ with probability min(1, R) and $X_{t+1}=w$ otherwise



In [ ]: