-

Author: BB

compscieng/compscieng_2_11/tst1.py

@@ -109,13 +109,16 @@ def exo1(x0,W):
         d[d==inf] = 0
         imageplot(d)
         plt.set_cmap('jet')        
-    plt.savefig('out-1.png')
+    plt.savefig('out-360.png')
     return D
 
 n = 40
 W = ones( (n,n) )
-print (W.shape)
 x0 = [int(n/2), int(n/2)]
-print (x0)
 D = exo1(x0,W)
-print (D)
+
+plt.figure()
+displ = lambda D: cos(2*pi*5*D/max(D.flatten()) )
+imageplot(displ(D))
+plt.savefig('out-390.png')
+

compscieng/compscieng_2_11/tst2.py

@@ -0,0 +1,139 @@
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from numpy import *
+from scipy import ndimage
+
+def imageplot(f, str='', sbpt=[]):
+    """
+        Use nearest neighbor interpolation for the display.
+    """
+    if sbpt != []:
+        plt.subplot(sbpt[0], sbpt[1], sbpt[2])
+    imgplot = plt.imshow(f, interpolation='nearest')
+    imgplot.set_cmap('gray')
+    plt.axis('off')
+
+def perform_dijstra_fm(W, pstart, niter=np.inf, method='dijstr', bound='sym', svg_rate=10):
+    n = W.shape[0]
+    neigh = np.array([[1, -1, 0, 0], [0, 0,  1, -1]])
+
+    def symmetrize(x,n):
+        if (x<0):
+            x = -x;
+        elif (x>=n):
+            x = 2*(n-1)-x
+        return x
+
+    if bound=='per':
+        boundary = lambda x: np.mod(x,n)
+    else:
+        boundary = lambda x: [symmetrize(x[0],n), symmetrize(x[1],n)] # todo
+
+    ind2sub1 = lambda k: [int( (k-np.fmod(k,n))/n ), np.fmod(k,n)]
+    sub2ind1 = lambda u: int( u[0]*n + u[1] )
+    Neigh = lambda k,i: sub2ind1(boundary(ind2sub1(k) + neigh[:,i]))
+    extract   = lambda x,I: x[I]
+    extract1d = lambda x,I: extract(x.flatten(),I)
+
+    nstart = pstart.shape[1]
+    I = list( np.zeros( (nstart, 1) ) )
+    for i in np.arange(0, nstart):
+        I[i] = int( sub2ind1(pstart[:,i]) )
+
+    D = np.zeros( (n,n) ) + np.inf # current distance
+    for i in np.arange(0, nstart):
+        D[pstart[0,i],pstart[1,i]] = 0
+
+    S = np.zeros( (n,n) )
+    for i in np.arange(0, nstart):
+        S[pstart[0,i],pstart[1,i]] = 1 # open
+
+    iter = 0
+    q = 100  # maximum number of saves
+    Dsvg = np.zeros( (n,n,q) )
+    Ssvg = np.zeros( (n,n,q) )
+    while ( not(I==[]) & (iter<=niter) ):
+        iter = iter+1;
+        if iter==niter:
+            break
+        j = np.argsort( extract1d(D,I)  )
+        if np.ndim(j)==0:
+            j = [j]
+        j = j[0]
+        i = I[j]
+        a = I.pop(j)
+        u = ind2sub1(i);
+        S[u[0],u[1]] = -1
+        J = []
+        for k in np.arange(0,4):
+            j = Neigh(i,k)
+            if extract1d(S,j)!=-1:
+                J.append(j)
+                if extract1d(S,j)==0:
+                    u = ind2sub1(j)
+                    S[u[0],u[1]] = 1
+                    I.append(j)
+        DNeigh = lambda D,k: extract1d(D,Neigh(j,k))
+        for j in J:
+            dx = min(DNeigh(D,0), DNeigh(D,1))
+            dy = min(DNeigh(D,2), DNeigh(D,3))
+            u = ind2sub1(j)
+            w = extract1d(W,j);
+            if method=='dijstr':
+                D[u[0],u[1]] = min(dx + w, dy + w)
+            else:
+                Delta = 2*w - (dx-dy)**2
+                if (Delta>=0):
+                    D[u[0],u[1]] = (dx + dy + np.sqrt(Delta))/ 2
+                else:
+                    D[u[0],u[1]] = min(dx + w, dy + w)
+        t = iter/svg_rate
+        if (np.mod(iter,svg_rate)==0) & (t<q):
+            print (t)
+            Dsvg[:,:,int(t-1)] = D
+            Ssvg[:,:,int(t-1)] = S
+
+    Dsvg = Dsvg[:,:,:int(t-1)]
+    Ssvg = Ssvg[:,:,:int(t-1)]
+    return (D,Dsvg,Ssvg);
+
+def exo1(x0,W):
+    n = W.shape[0]
+    pstart = transpose(array([x0]))
+    [D,Dsvg,Ssvg] = perform_dijstra_fm(W, pstart, inf,'dijstr', 'sym',n*6)
+    plt.figure();
+    for i in arange(0,4):
+        plt.subplot(2, 2, i+1)
+        d = Dsvg[:,:,i]
+        d[d==inf] = 0
+        imageplot(d)
+        plt.set_cmap('jet')        
+    plt.savefig('out-360.png')
+    return D
+
+def exo2(x0,W):
+    n = W.shape[0]
+    pstart = transpose(array([x0]))
+    [D,Dsvg,Ssvg] = perform_dijstra_fm(W, pstart, inf,'fm', 'sym',n*6)
+    plt.figure();
+    for i in arange(0,4):
+        plt.subplot(2, 2, i+1)
+        d = Dsvg[:,:,i]
+        d[d==inf] = 0
+        imageplot(d)
+        plt.set_cmap('jet')        
+    plt.savefig('out-450.png')
+    return D
+
+n = 40
+W = ones( (n,n) )
+x0 = [int(n/2), int(n/2)]
+
+D = exo2(x0,W)
+
+plt.figure()
+displ = lambda D: cos(2*pi*5*D/max(D.flatten()) )
+imageplot(displ(D))
+plt.savefig('out-480.png')
+