modified microscopy analysis

code/local_features.py

-"""
-Helping functions for extracting features from images
-"""
-
-import skimage.io
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy.ndimage
-
-def get_gauss_feat_im(im, sigma=1, normalize=False):
-    """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: a 3D array of size (r,c,15) with a 15-dimentional feature
-            vector for every image pixel.
-    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))
-        imfeat *= (1/np.std(imfeat, axis=(0,1)))
-    
-    return imfeat
-
-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
-
-#%%
-
-if __name__ == '__main__':
-    
-    filename = '../example_data/training_image.png'
-    I = skimage.io.imread(filename)
-    I = I.astype(np.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])
-            ax[i][j].set_title(f'layer {5*i+j}')
-            
-            
-    # features based on image patches
-    I = I[I.shape[0]//2:I.shape[0]//2+20,I.shape[1]//2:I.shape[1]//2+20] # 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()
-    ax.imshow(I)
-            
-    fig,ax = plt.subplots(3,3)
-    for j in range(3):
-        for i in range(3):
-            ax[i][j].imshow(pf[3*i+j])
-            ax[i][j].set_title(f'layer {3*i+j}')
-   
-    

code/microscopy_analysis.py

@@ -10,7 +10,6 @@ import os
 import skimage.io as io
 import skimage.transform
 import numpy as np
-import local_features as lf
 
 import matplotlib.pyplot as plt
 from math import ceil
@@ -310,37 +309,6 @@ def rgb2cmy(img):
 
 #%%Functions for feature analysis and visualisation
 
-# computes the max projection image and the features into a list
-def compute_max_img_feat(in_dir, ext, base_name, sigma, sc_fac, save = False, abs_intensity = True):
-    dir_list = [dI for dI in os.listdir(in_dir) if (os.path.isdir(os.path.join(in_dir,dI)) and dI[0:len(base_name)]==base_name)]
-    dir_list.sort()
-    max_img_list = []
-    for d in dir_list:
-        image_dir = in_dir + d + '/'
-        max_img = get_max_img(image_dir, ext)
-        
-        if save:
-            dir_name_out_img = in_dir.replace('data','maxProjImages')
-            os.makedirs(dir_name_out_img, exist_ok = True)
-
-            io.imsave( dir_name_out_img + d + '.png', max_img.astype('uint8'))
-            
-        max_img_list += [skimage.transform.rescale(max_img, sc_fac, multichannel=True)]
-    
-    feat_list = []
-    for max_img in max_img_list:
-        r,c = max_img.shape[:2]
-        feat = np.zeros((r*c,45))
-        for i in range(0,3):
-            feat_tmp = lf.get_gauss_feat_im(max_img[:,:,i], sigma)
-            feat[:,i*15:(i+1)*15] = feat_tmp.reshape(-1,feat_tmp.shape[2])
-        if not(abs_intensity):
-            feat_list += [np.concatenate((feat[:,1:15],feat[:,16:30],feat[:,31:45]),axis = 1)]
-        else:    
-            feat_list += [feat]
-    return max_img_list, feat_list
-
-
 # computes a histogram of features from a kmeans object
 def compute_assignment_hist(feat_list, kmeans, background_feat = None):
     assignment_list = []