Added material for chapter 5

Chapter03/local_features.py

+"""
+Helping functions for extracting features from images
+
+Vedrana Andersen Dahl (vand@dtu.dk) 
+Anders Bjorholm Dahl (abda@dtu.dk)
+"""
+
+#%%
+import numpy as np
+import scipy.ndimage
+
+def get_gauss_feat_im(im, sigma=1, normalize=True):
+    """Gauss derivative feaures for every image pixel.
+    Arguments:
+        image: a 2D image, shape (r, c).
+        sigma: standard deviation for Gaussian derivatives.
+        normalize: flag indicating normalization of features.
+    Returns:
+        imfeat: 3D array of size (r, c, 15) with a 15-dimentional feature
+             vector for every pixel in the image.
+    Author: vand@dtu.dk, 2022
+    """
+    orders = [0, 
+        [0, 1], [1, 0], 
+        [0, 2], [1, 1], [2, 0], 
+        [0, 3], [1, 2], [2, 1], [3, 0],
+        [0, 4], [1, 3], [2, 2], [3, 1], [4, 0]]
+
+    imfeat = [scipy.ndimage.gaussian_filter(im, sigma, o) for o in orders]    
+    imfeat = np.stack(imfeat, axis=2)
+   
+    if normalize:
+        imfeat -= np.mean(imfeat, axis=(0, 1))
+        std =  np.std(imfeat, axis=(0, 1))
+        std[std==0] = 1
+        imfeat /= std
+    
+    return imfeat
+
+def get_gauss_feat_multi(im, sigmas=[1, 2, 4], normalize = True):
+    '''Multi-scale Gauss derivative feaures for every image pixel.
+    Arguments:
+        image: a 2D image, shape (r, c).
+        sigma: list of standard deviations for Gaussian derivatives.
+        normalize: flag indicating normalization of features.
+    Returns:
+        imfeat: a a 3D array of size (r*c, n_scale, 15) with n_scale features in 
+            each pixels, and n_scale is length of sigma. Each pixel contains a 
+            feature vector and feature image is size (r, c, 15*n_scale).
+    Author: abda@dtu.dk, 2021
+
+    '''
+    imfeats = []
+    for s in sigmas:
+        feat = get_gauss_feat_im(im, s, normalize)
+        imfeats.append(feat.reshape(-1, feat.shape[2]))
+    
+    imfeats = np.asarray(imfeats).transpose(1,0,2)
+    return imfeats
+
+
+def im2col(im, patch_size=[3, 3], stepsize=1):
+    """Rearrange image patches into columns
+    Arguments:
+        image: a 2D image, shape (r,c).
+        patch size: size of extracted paches.
+        stepsize: patch step size.
+    Returns:
+        patches: a 2D array which in every column has a patch associated 
+            with one image pixel. For stepsize 1, number of returned column 
+            is (r-patch_size[0]+1)*(c-patch_size[0]+1) due to bounary. The 
+            length of columns is pathc_size[0]*patch_size[1].
+    """
+    
+    r, c = im.shape
+    s0, s1 = im.strides    
+    nrows =r - patch_size[0] + 1
+    ncols = c - patch_size[1] + 1
+    shp = patch_size[0], patch_size[1], nrows, ncols
+    strd = s0, s1, s0, s1
+
+    out_view = np.lib.stride_tricks.as_strided(im, shape=shp, strides=strd)
+    return out_view.reshape(patch_size[0]*patch_size[1], -1)[:, ::stepsize]
+
+
+def ndim2col(im, block_size=[3, 3], stepsize=1):
+    """Rearrange image blocks into columns for N-D image (e.g. RGB image)"""""
+    if(im.ndim == 2):
+        return im2col(im, block_size, stepsize)
+    else:
+        r, c, l = im.shape
+        patches = np.zeros((l * block_size[0] * block_size[1],
+                            (r - block_size[0] + 1) * (c - block_size[1] + 1)))
+        for i in range(l):
+            p = block_size[0] * block_size[1]
+            patches[i * p : (i + 1) * p, :] = im2col(
+                im[:, :, i], block_size, stepsize)
+        return patches
+
+#%%
+import skimage.io
+import matplotlib.pyplot as plt
+
+if __name__ == '__main__':
+    
+    #%% features based on gaussian derivatives
+    filename = '../../data/week3/3labels/training_image.png'
+    I = skimage.io.imread(filename).astype(float)/255
+    I = I[200:400, 200:400] # smaller image such that we can see 
+    fig, ax = plt.subplots()
+    ax.imshow(I)
+    
+    sigma = 5
+    gf = get_gauss_feat_im(I, sigma)
+    
+    fig,ax = plt.subplots(5, 3, figsize=(15, 25))
+
+
+    ax = ax.ravel()
+    for i, a in enumerate(ax):
+        a.imshow(gf[..., i], cmap='jet')
+        a.set_title(f'layer {i}')
+    plt.show()
+            
+            
+    #%% features based on image patches
+    I = skimage.io.imread(filename).astype(float)/255
+    I = I[300:320, 400:420] # smaller image such that we can see 
+    fig, ax = plt.subplots()
+    ax.imshow(I)
+
+    pf = im2col(I, [3, 3])
+    pf = pf.reshape((9, I.shape[0]-2, I.shape[1]-2))
+            
+    fig,ax = plt.subplots(3,3)
+    for j in range(3):
+        for i in range(3):
+            ax[i][j].imshow(pf[3*i+j], cmap='jet')
+            ax[i][j].set_title(f'layer {3*i+j}')
+            
+    plt.show()
+    
+
+# %%
+I.shape
+
+# %%

Chapter05/circles_modeling_EMPTY.py

+#%%
+
+import skimage.io
+import numpy as np
+import scipy.ndimage
+import matplotlib.pyplot as plt
+
+#%%
+
+def segmentation_energy(S, D, mu, beta):
+    # TODO -- add your code here
+    
+    # likelihood energy
+    U1 = 0
+    
+    # prior energy
+    U2 = 0
+    
+    return U1, U2
+
+def segmentation_histogram(ax, D, S, edges=None):
+    '''
+    Plot histogram for grayscale data and each segmentation label.
+    '''
+    if edges is None:
+        edges = np.linspace(D.min(), D.max(), 100)
+    ax.hist(D.ravel(), bins=edges, color = 'k')
+    centers = 0.5 * (edges[:-1] + edges[1:])
+    for k in range(S.max() + 1):
+        ax.plot(centers, np.histogram(D[S==k].ravel(), edges)[0])
+        
+
+path = '../../data/week5/'
+D = skimage.io.imread(path + 'noisy_circles.png').astype(float)
+
+# Ground-truth segmentation.
+GT = skimage.io.imread(path + 'noise_free_circles.png')
+(mu, S_gt) = np.unique(GT, return_inverse=True)
+S_gt = S_gt.reshape(D.shape)
+
+segmentations = [S_gt]  # list where I'll place different segmentations
+
+#%% Find some configurations (segmentations) using conventional methods.
+
+# Simple thresholding
+S_t = np.zeros(D.shape, dtype=int) + (D > 100) + (D > 160) # thresholded
+segmentations += [S_t]
+
+# Gaussian filtering followed by thresholding
+D_s = scipy.ndimage.gaussian_filter(D, sigma=1, truncate=3, mode='nearest')
+S_g = np.zeros(D.shape, dtype=int) + (D_s > 100) + (D_s > 160) 
+segmentations += [S_g]
+
+# Median filtering followed by thresholding
+D_m = scipy.ndimage.median_filter(D, size=(5, 5), mode='reflect')
+S_t = np.zeros(D.shape, dtype=int) + (D_m > 100) + (D_m > 160) # thresholded
+segmentations += [S_t]
+
+
+#%% visualization
+fig, ax = plt.subplots()
+ax.imshow(D, vmin=0, vmax=255, cmap=plt.cm.gray)
+plt.show()
+
+
+fig, ax = plt.subplots(3, len(segmentations), figsize=(10, 10))
+beta = 100
+
+for i, s in enumerate(segmentations):
+    ax[0][i].imshow(s)
+
+    V1, V2 = segmentation_energy(s, D, mu, beta)
+    ax[0][i].set_title(f'likelihood: {V1:.2g}\nprior: {V2}\nposterior: {V1+V2:.2g}')
+    
+    segmentation_histogram(ax[1][i], D, s)
+    ax[1][i].set_xlabel('Intensity')
+    ax[1][i].set_ylabel('Count')
+    
+    err = S_gt - s
+    ax[2][i].imshow(err, vmin=-2, vmax=2, cmap=plt.cm.bwr)
+    ax[2][i].set_title(f'Pixel error: {(err != 0).sum()}')
+
+fig.tight_layout()
+plt.show()   
+
+
+# %%

Chapter05/gender_labeling.py

+#%% Example of max flow min cut for gender labeling problem
+
+import maxflow
+import numpy as np
+
+#%%
+d = np.array([179, 174, 182, 162, 175, 165]) # heights (data)
+mu = [181, 165] # means of two classes
+beta = 100 # weight of the prior term
+w_s = (d-mu[0])**2 # source weight
+w_t = (d-mu[1])**2 # sink weights
+N = len(d) # number of graph nodes
+
+# Create a graph with integer capacities.
+g = maxflow.Graph[int]()
+# Add (non-terminal) nodes and retrieve an index for each node
+nodes = g.add_nodes(N)
+# Create edges between nodes
+for i in range(N-1):
+    g.add_edge(nodes[i], nodes[i+1], beta, beta)
+# Set the capacities of the terminal edges.
+for i in range(N):
+    g.add_tedge(nodes[i], (d[i]-mu[1])**2, (d[i]-mu[0])**2)
+# Run the max flow algorithm
+flow = g.maxflow()
+print(f'Maximum flow: {flow}')
+
+# displaying the results
+labeling = [g.get_segment(nodes[i]) for i in range(N)]
+gend = 'MF'
+
+for i in range(0,N):
+    print(f'Person {i} is estimated as {gend[labeling[i]]}') 
+
+
+
+
+
+
+
+
+# %%