added week 4 solution

Week04/QUIZ_week4_solution.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "608c7b2a",
+   "metadata": {},
+   "source": [
+    "# QUIZ WEEK 4 SOLUTIONS"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a4d9315d",
+   "metadata": {},
+   "source": [
+    "## Question 1, features size"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "56f027e9",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "features_size = 2394000\n"
+     ]
+    }
+   ],
+   "source": [
+    "p = 5\n",
+    "h = 256\n",
+    "w = 384\n",
+    "\n",
+    "features_size = (h - p + 1) * (w - p + 1) * (p**2)\n",
+    "print(f'{features_size = }')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0649178f",
+   "metadata": {},
+   "source": [
+    "## Question 2, cluster probability"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "4329d615",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "cluster_probability = 0.2603369065849923\n"
+     ]
+    }
+   ],
+   "source": [
+    "cluster_probability = 1 - 483/653\n",
+    "print(f'{cluster_probability = }')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3179e4d3",
+   "metadata": {},
+   "source": [
+    "## Question 3, dominant segmentation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "ca50755a",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "dominant_segmentation = 3\n"
+     ]
+    }
+   ],
+   "source": [
+    "# dominant semgentation\n",
+    "import numpy as np\n",
+    "\n",
+    "P1 = np.array([[0.01, 0.11, 0.21, 0.31, 0.41],\n",
+    "    [0.04, 0.14, 0.24, 0.34, 0.44],\n",
+    "    [0.07, 0.17, 0.27, 0.37, 0.47]])\n",
+    "P2 = np.array([[0.23, 0.26, 0.29, 0.32, 0.35],\n",
+    "    [0.30, 0.33, 0.36, 0.39, 0.42],\n",
+    "    [0.37, 0.40, 0.43, 0.46, 0.49]])\n",
+    "\n",
+    "P = np.stack([P1, P2, 1-(P1+P2)])\n",
+    "S = np.argmax(P, axis=0) + 1  # segmentation into 1, 2, 3\n",
+    "counts = np.unique(S, return_counts=True)[1]  # how many times each label occurs\n",
+    "dominant_segmentation = np.argmax(counts) + 1\n",
+    "\n",
+    "print(f'{dominant_segmentation = }')\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.12"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "b35caa3764096924ce573458b7678f48285fe4447fc77b92472b406a014fb35c"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

Week04/feature_segmentatation_prob_smooth.py

+#%%
+# PREPARE
+"""
+Vedrana Andersen Dahl (vand@dtu.dk) 
+Anders Bjorholm Dahl (abda@dtu.dk)
+"""
+
+import skimage.io
+import numpy as np
+import matplotlib.pyplot as plt
+import sklearn.cluster
+import local_features as lf
+import scipy.ndimage
+
+def ind2labels(ind):
+    """ Helper function for transforming uint8 image into labeled image."""
+    return np.unique(ind, return_inverse=True)[1].reshape(ind.shape)
+
+path = '../data/week4/3labels/' # Change path to your directory
+
+#%% 
+# READ IN IMAGES
+training_image = skimage.io.imread(path + 'training_image.png')
+training_image = training_image.astype(float)
+training_labels = skimage.io.imread(path + 'training_labels.png')
+
+training_labels = ind2labels(training_labels)
+nr_labels = np.max(training_labels)+1 # number of labels in the training image
+
+fig, ax = plt.subplots(1, 2)
+ax[0].imshow(training_image, cmap=plt.cm.gray)
+ax[0].set_title('training image')
+ax[1].imshow(training_labels)
+ax[1].set_title('labels for training image')
+fig.tight_layout()
+plt.show()
+
+#%% 
+# TRAIN THE MODEL
+
+sigma = [1, 2 , 3]
+features = lf.get_gauss_feat_multi(training_image, sigma)
+features = features.reshape((features.shape[0], features.shape[1]*features.shape[2]))
+labels = training_labels.ravel()
+
+nr_keep = 15000 # number of features randomly picked for clustering 
+keep_indices = np.random.permutation(np.arange(features.shape[0]))[:nr_keep]
+
+features_subset = features[keep_indices]
+labels_subset = labels[keep_indices]
+
+nr_clusters = 1000 # number of feature clusters
+# for speed, I use mini-batches
+kmeans = sklearn.cluster.MiniBatchKMeans(n_clusters=nr_clusters, 
+                                         batch_size=2*nr_clusters, 
+                                         n_init='auto')
+kmeans.fit(features_subset)
+assignment = kmeans.labels_
+
+edges = np.arange(nr_clusters + 1) - 0.5 # histogram edges halfway between integers
+hist = np.zeros((nr_clusters, nr_labels))
+for l in range(nr_labels):
+    hist[:, l] = np.histogram(assignment[labels_subset == l], bins=edges)[0]
+sum_hist = np.sum(hist, axis=1)
+cluster_probabilities = hist/(sum_hist.reshape(-1, 1))
+
+fig, ax = plt.subplots(2, 1)
+legend_label = [f'label {x}' for x in range(nr_labels)]
+
+ax[0].plot(hist, '.', alpha=0.5, markersize=3)
+ax[0].set_xlabel('cluster id')
+ax[0].set_ylabel('number of features in cluster')
+ax[0].legend(legend_label)
+ax[0].set_title('features in clusters per label')
+ax[1].plot(cluster_probabilities, '.', alpha=0.5, markersize=3)
+ax[1].set_xlabel('cluster id')
+ax[1].set_ylabel('label probability for cluster')
+ax[1].legend(legend_label)
+ax[1].set_title('cluster probabilities')
+fig.tight_layout()
+plt.show()
+
+# Finished training
+
+#%% 
+# USE THE MODEL
+testing_image = skimage.io.imread(path + 'testing_image.png')
+testing_image = testing_image.astype(float)
+
+features_testing = lf.get_gauss_feat_multi(testing_image, sigma)
+features_testing = features_testing.reshape((features_testing.shape[0], features_testing.shape[1]*features_testing.shape[2]))
+labels = training_labels.ravel()
+
+assignment_testing = kmeans.predict(features_testing)
+
+probability_image = np.zeros((assignment_testing.size, nr_labels))
+for l in range(nr_labels):
+    probability_image[:, l] = cluster_probabilities[assignment_testing, l]
+probability_image = probability_image.reshape(testing_image.shape + (nr_labels, ))
+
+P_rgb = np.zeros(probability_image.shape[0:2]+(3, ))
+k = min(nr_labels, 3)
+P_rgb[:, :, :k] = probability_image[:, :, :k]
+
+fig, ax = plt.subplots(1, 2)
+ax[0].imshow(testing_image, cmap=plt.cm.gray)
+ax[0].set_title('testing image')
+ax[1].imshow(P_rgb)
+ax[1].set_title('probabilities for testing image as RGB')
+fig.tight_layout()
+plt.show()
+
+#%%
+# SMOOTH PROBABILITY MAP
+
+sigma = 3 # Gaussian smoothing parameter
+
+seg_im_max = np.argmax(P_rgb, axis = 2)
+c = np.eye(P_rgb.shape[2])
+P_rgb_max = c[seg_im_max]
+
+probability_smooth = np.zeros(probability_image.shape)
+for i in range(0, probability_image.shape[2]):
+    probability_smooth[:, :, i] = scipy.ndimage.gaussian_filter(probability_image[:, :, i], sigma, order=0)
+seg_im_smooth = np.argmax(probability_smooth, axis=2)
+
+probability_smooth_max = c[seg_im_smooth]
+
+P_rgb_smooth = np.zeros(probability_smooth_max.shape[0:2]+(3, ))
+k = min(nr_labels, 3)
+P_rgb_smooth[:, :, :k] = probability_smooth[:, :, :k]
+P_rgb_smooth_max = np.zeros(probability_smooth_max.shape[0:2]+(3, ))
+P_rgb_smooth_max[:, :, :k] = probability_smooth_max[:, :, :k]
+
+# Display result
+fig, ax = plt.subplots(2, 4, sharex=True, sharey=True)
+ax[0][0].imshow(P_rgb[:, :, 0])
+ax[0][1].imshow(P_rgb[:, :, 1])
+ax[0][2].imshow(P_rgb[:, :, 2])
+ax[0][3].imshow(P_rgb_max)
+ax[1][0].imshow(P_rgb_smooth[:, :, 0])
+ax[1][1].imshow(P_rgb_smooth[:, :, 1])
+ax[1][2].imshow(P_rgb_smooth[:, :, 2])
+ax[1][3].imshow(P_rgb_smooth_max)
+fig.tight_layout()
+plt.show()
+
+# %%

Week04/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, 2020
+    """
+      
+    r, c = im.shape
+    imfeat = np.zeros((r, c, 15))
+    imfeat[:, :, 0] = scipy.ndimage.gaussian_filter(im, sigma, order=0)
+    imfeat[:, :, 1] = scipy.ndimage.gaussian_filter(im, sigma, order=[0, 1])
+    imfeat[:, :, 2] = scipy.ndimage.gaussian_filter(im, sigma, order=[1, 0])
+    imfeat[:, :, 3] = scipy.ndimage.gaussian_filter(im, sigma, order=[0, 2])
+    imfeat[:, :, 4] = scipy.ndimage.gaussian_filter(im, sigma, order=[1, 1])
+    imfeat[:, :, 5] = scipy.ndimage.gaussian_filter(im, sigma, order=[2, 0])
+    imfeat[:, :, 6] = scipy.ndimage.gaussian_filter(im, sigma, order=[0, 3])
+    imfeat[:, :, 7] = scipy.ndimage.gaussian_filter(im, sigma, order=[1, 2])
+    imfeat[:, :, 8] = scipy.ndimage.gaussian_filter(im, sigma, order=[2, 1])
+    imfeat[:, :, 9] = scipy.ndimage.gaussian_filter(im, sigma, order=[3, 0])
+    imfeat[:, :, 10] = scipy.ndimage.gaussian_filter(im, sigma, order=[0, 4])
+    imfeat[:, :, 11] = scipy.ndimage.gaussian_filter(im, sigma, order=[1, 3])
+    imfeat[:, :, 12] = scipy.ndimage.gaussian_filter(im, sigma, order=[2, 2])
+    imfeat[:, :, 13] = scipy.ndimage.gaussian_filter(im, sigma, order=[3, 1])
+    imfeat[:, :, 14] = scipy.ndimage.gaussian_filter(im, sigma, order=[4, 0])
+
+    if normalize:
+        imfeat -= np.mean(imfeat, axis=(0, 1))
+        im_std = np.std(imfeat, axis=(0, 1))
+        im_std[im_std<10e-10] = 1
+        imfeat /= im_std
+    
+    return imfeat
+
+def get_gauss_feat_multi(im, sigma = [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 i in range(0, len(sigma)):
+        feat = get_gauss_feat_im(im, sigma[i], 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):
+            patches[i * block_size[0] * block_size[1] : (i+1) * block_size[0] * block_size[1], 
+                    :] = im2col(im[:, :, i], block_size, stepsize)
+        return patches
+
+#%%
+import skimage.io
+import matplotlib.pyplot as plt
+
+if __name__ == '__main__':
+    
+    filename = '../../../../Data/week3/3labels/training_image.png'
+    I = skimage.io.imread(filename)
+    I = I.astype(float)
+    
+    
+    # features based on gaussian derivatives
+    sigma = 1;
+    gf = get_gauss_feat_im(I, sigma)
+    
+    fig, ax = plt.subplots(3, 5)
+    for j in range(5):
+        for i in range(3):
+            ax[i][j].imshow(gf[:, :, 5 * i + j], cmap='jet')
+            print(np.std(gf[:, :, 5 * i + j]))
+            ax[i][j].set_title(f'layer {5 * i + j}')
+            
+            
+    # features based on image patches
+    I = I[300:320, 400:420] # smaller image such that we can see the difference in layers
+    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}')
+   
+