Initial commit

ClickCells/README

+———————————————————————
+ClickCells Python tools
+———————————————————————
+Python scripts to help create training data for cell density napping in tracer-based histological slices.
+
+********************
+About
+********************
+
+
+>> select_zones.py 
+Inputs an image and allows the selection of a zone of interest.
+
+>> split_zones.py
+Slice an image into many sub-images. 
+
+>> click_cells.npy
+Loads .npy sub-images created from  slice_zone.py in order for the user to click on visible cells. These cells' coordinates are then be exported into a .txt file.
+
+
+********************
+Dependencies
+********************
+Python 3.*
+Pillow
+Numpy
+Matplotlib
+Pandas
+
+
+
+********************
+Authors
+********************
+
+- Nadir Basma
+- Saad Jbabdi
+FMRIB, University of Oxford, UK
\ No newline at end of file

ClickCells/click_cells.py

+#!/usr/bin/env python3
+from PIL import Image
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import os.path as op
+import re
+import glob
+import argparse
+import shutil
+import random
+
+
+class CellPicker(object):
+    def __init__(self, figure, axis):
+        self.points = []
+        self.figure = figure
+        self.axis   = axis
+        self.figure.canvas.mpl_connect('button_press_event', self.on_press)
+        self.figure.canvas.mpl_connect('key_press_event', self.on_key_undo)
+        self.figure.canvas.mpl_connect('key_press_event', self.on_key_quit)
+        self.figure.canvas.mpl_connect('key_press_event', self.on_key_dunno)
+
+        
+    def on_press(self, ev):
+        self.points.append((ev.xdata, ev.ydata))
+        points_array = np.array(self.points)
+        self.axis.plot(points_array[:, 0],
+                       points_array[:, 1],
+                       marker='o', linestyle='None', markersize=5, color="red")
+        self.figure.canvas.draw()    
+
+    def on_key_quit(self, event):
+        if event.key == 'q':
+            plt.close()
+            exit()
+    
+    def on_key_undo(self, event):
+        if event.key == 'u':
+            self.points.pop(-1)
+            points_array = np.array(self.points)
+            self.axis.clear()
+            self.axis.plot(points_array[:, 0],
+                           points_array[:, 1], marker='o',
+                           linestyle='None', markersize=5, color="red")
+            self.figure.canvas.draw()
+
+    def on_key_dunno(self, event):
+        if event.key == 'n':
+            plt.close()
+            
+
+        
+def main():
+    
+    # Parse command line arguments
+    parser = argparse.ArgumentParser(
+        "Click on cells and save to txt file"
+    )
+    parser.add_argument("-i","--input_folder",
+                        required=True,
+                        help="Input folder (_splitted).")
+    parser.add_argument("-o","--output_cell_coordinates",
+                        required=True,
+                        help="Output file name.")
+    parser.add_argument("--shuffle", action='store_true', default=False,
+                        dest='shuffle',help="Load sub-images in random order.")
+    parser.add_argument("--empty_zone", action='store_true', default=False,
+                        dest="empty_zone",help="Entire zone is empty")
+    args = parser.parse_args()
+    
+    # Find Numpy array files
+    infolder = args.input_folder 
+    outfile = args.output_cell_coordinates
+    
+    files = glob.glob(op.join(infolder,'*_w_*_h_*.npy'))
+    if args.shuffle == True:
+        random.shuffle(files)
+
+    create_header = True
+    if op.exists(outfile):
+        print("File {} exists. Overwrite/Append/Exit?[O,A,E]".format(outfile))
+        response = input()
+        if response.upper() == "O":
+            os.remove(outfile)
+        elif response.upper() == "E":
+            print("Exiting without doing anything")
+            exit()
+        elif response.upper() == "A":
+            create_header = False
+
+    if create_header == True:
+        with open(outfile,'w') as f:
+            f.write('Sub-Image-File\tXcoord\tYcoord\n')
+    
+    for file in files:
+        #First, get h and w values corresponding to sub-image from filename.
+        #This is in order to export x,y coordinates relative to original brain_slice image.
+        res = re.findall("w_(\d+).(\d+)_h_(\d+).(\d+)", file)[0]
+        w = round(float(res[0]+"."+res[1]))
+        h = round(float(res[2]+"."+res[3]))
+
+        if args.empty_zone == True:
+            with open(outfile,'a') as f:
+                f.write('%s\tNaN\tNaN\n' %file)
+        else:
+            #prepare figure and load .npy files in as an image, ready to show
+            fig = plt.figure()
+            ax = fig.add_subplot(111)
+            ax.set_title("Click on cells. \n"
+                         "'q'=quit, 'u'=undo, 'n'=dunno \nClose figure if you can't see any cells")
+            im = ax.imshow(np.load(file))        
+            
+            ax.set_xlim(ax.get_xlim())
+            ax.set_ylim(ax.get_ylim())
+            
+            
+            #Instantiate CellPicker class
+            p = CellPicker(fig,ax)
+            
+            
+            #Once past this we move to the next image...
+            plt.show()
+            
+            # Write x,y coordinates to the output file,
+            # adding h and w values in order to translate into coordinates
+            # of original image.
+            
+            with open(outfile,'a') as f:
+                if not p.points:
+                    f.write('%s\tNaN\tNaN\n' %file)
+                else:        
+                    for point in p.points:
+                        f.write('%20s \t %12f \t %12f \n' % (os.path.abspath(file), point[1] + w, point[0] + h))
+                
+                        
+            
+
+if __name__ == '__main__':
+    main()

ClickCells/create_db.py

+#!/usr/bin/env python3
+
+import argparse
+import numpy as np
+import pandas as pd
+import re
+
+DB_IMAGE_RES = 64
+
+def check_imshape(shape):
+    sx, sy, _ = shape
+    if sx != sy:
+        return False
+    if (sx % DB_IMAGE_RES != 0) or (sy % DB_IMAGE_RES !=0):
+        return False
+    return True
+
+def append_file_content(fname,image_list,count_list):
+    df = pd.read_table(fname)
+    udf = df.groupby('Sub-Image-File').count()
+
+    for f in udf.index:
+        # Load Numpy array
+        im = np.load(f.strip())
+        if not check_imshape(im.shape):
+            print("Error: Bad Image dimensions. Must be square and multiple of {}".format(DB_IMAGE_RES))
+            
+        sizx, sizy, _ = im.shape
+
+        # Split into sub-zones
+        size_ratio = sizx//DB_IMAGE_RES
+
+        im2 = im.reshape(size_ratio,sizx//size_ratio,size_ratio,sizy//size_ratio,3)
+        im3 = im2.transpose(0,2,1,3,4).reshape(size_ratio**2,sizx//size_ratio,sizy//size_ratio,3)
+        image_list.append(im3)
+
+        res = re.findall("w_(\d+).(\d+)_h_(\d+).(\d+)", f)[0]
+        W = round(float(res[0]+"."+res[1]))
+        H = round(float(res[2]+"."+res[3]))
+
+        count_cells = np.zeros((size_ratio, size_ratio),dtype=int)
+        for indiv_cells in df.values[df['Sub-Image-File']==f]:            
+            if np.isnan(indiv_cells[1]):
+                pass
+            else:
+                w = float(indiv_cells[1])-W
+                h = float(indiv_cells[2])-H
+
+                count_cells[int(w//(sizx/size_ratio)),int(h//(sizy/size_ratio))] += 1
+
+        count_list.append(count_cells.flatten())
+
+    return len(udf.index)
+    
+def create_db(file_list,outfile):
+    image_list = []
+    count_list = []
+    total = 0
+    for f in file_list:
+        total += append_file_content(f,image_list,count_list)
+
+    shape = image_list[0].shape[1:]
+    count_list = np.array(count_list).flatten()
+    np.savez(outfile,counts=count_list,
+             images=np.array(image_list).reshape(-1,*shape))
+
+    
+
+    print("Created DB with {} Images from a list of {}.".format(len(count_list),total))
+    
+def main():
+        
+    # Parse command line arguments
+    parser = argparse.ArgumentParser(
+        "Create DB from clicked textfiles"
+    )
+    parser.add_argument("outfile",
+                        help="Output file name")    
+    parser.add_argument("file_list",
+                        help="List of clicky text files",
+                        nargs='+')
+    args = parser.parse_args()
+
+
+    create_db(args.file_list,args.outfile)
+
+if __name__ == '__main__':
+    main()

ClickCells/select_zones.py

+#!/usr/bin/env python3
+
+from PIL import Image
+
+import matplotlib as mpl
+
+mpl.use('macosx')
+
+from matplotlib.widgets import RectangleSelector
+import numpy as np
+import os
+import os.path as op
+import shutil
+import matplotlib.pyplot as plt
+import argparse
+import sys
+
+
+Image.MAX_IMAGE_PIXELS = 5282426800
+
+class ZoneSelector(object):
+    def __init__(self,ax,outbase,img,res):
+        self.box         = []
+        self.axis        = ax
+        self.zone_number = 0
+        self.image       = img
+        self.outbase     = outbase
+        self.res         = res
+        self.rs          = RectangleSelector(self.axis,
+                                             self.select_callback,
+                                             drawtype='box',
+                                             rectprops = dict(facecolor='orange',
+                                                              edgecolor = 'red', alpha=0.2, fill=True))
+        
+
+    def select_callback(self, eclick, erelease):
+        x1, y1 = eclick.xdata, eclick.ydata
+        x2, y2 = erelease.xdata, erelease.ydata
+        zone = np.array((np.min([x1,x2]),
+                         np.min([y1,y2]),
+                         np.max([x1,x2]),
+                         np.max([y1,y2])))*self.res
+        crop = self.image.crop(zone)        
+        crop.save(self.outbase + "_zone{:03d}.tiff".format(self.zone_number), "TIFF")
+        with open(self.outbase + "_zone{:03d}.txt".format(self.zone_number),"w") as fcoord:
+            fcoord.write('{} {} {} {}'.format(*zone))
+        self.zone_number +=1
+
+        
+
+def main():
+    parser = argparse.ArgumentParser(
+        "Select zone of interedt in an image"
+    )
+    required = parser.add_argument_group('Required arguments')
+    required.add_argument("-i","--input", required=True, help="Input file name of main image including extension.")
+    required.add_argument("-o","--output_folder", required=True, help="Name of output folder.")
+    args = parser.parse_args()
+
+    image_name = args.input
+    outfolder = args.output_folder
+
+    # Deal with existing output
+    if op.exists(outfolder):
+        print("Folder '{}' exists. Are you sure you want to delete it? [Y,N]".format(outfolder))
+        response = input()
+        if response.upper() == "Y":
+            shutil.rmtree(outfolder)
+        else:
+            print("Exiting without doing anything....")
+            sys.exit(1)
+            
+    os.mkdir(outfolder)
+
+    # Read image
+    im  = Image.open(image_name)        
+    image_base = op.basename(image_name)
+
+    # Resample image before showing it
+    res = 8
+    imr = np.asarray(im)[::res, ::res, :]
+
+    ax = plt.figure().add_subplot(111)
+    ax.imshow(imr)
+    
+    plt.title('CLick and drag to draw zones')
+
+    zs = ZoneSelector(ax,op.join(outfolder,image_base),im,res)
+
+    plt.show()
+
+
+if __name__ == '__main__':
+    main()
+
+
+

ClickCells/split_zones.py

+#!/usr/bin/env python3
+
+from PIL import Image
+import numpy as np
+import matplotlib.pyplot as plt
+import argparse
+import glob
+import os
+import os.path as op
+import shutil
+
+def main():
+    parser = argparse.ArgumentParser(
+        "Slice a big image into small sub-images. "
+    )
+    required = parser.add_argument_group('Required arguments')
+    required.add_argument("-i","--input_folder", required=True, help="Input folder name")
+    parser.add_argument("-w","--width", type=int, default=128, help="Width of sub-images. default=128.")
+    parser.add_argument('-t', "--height", type=int, default=128, help="Height of sub-images. default=128")
+    parser.add_argument("--overwrite", action='store_true', default=False,
+                        dest='overwrite',help="Overwrite existing slicings?")
+    args = parser.parse_args() 
+    
+    # Look inside folder for files called 'zone???.txt'
+    listing = os.listdir(args.input_folder)
+    zones = glob.glob(op.join(args.input_folder, '*_zone???.tiff'))
+    print(zones)
+    
+    for zone in zones:
+        zone_base = op.basename(zone)[:-5]
+        zone_image_name = zone
+        zone_coord_name = zone[:-5]+".txt"
+        zone_slices = zone[:-5]+"_splitted_{}x{}".format(args.width,args.height)
+        if not op.exists(zone_coord_name):
+            print("File {} does not have an associated coord file {}. "
+                  "Skipping this zone".format(zone_image_name,zone_coord_name))
+            continue
+        
+        if op.exists(zone_slices):
+            if args.overwrite == True:
+                shutil.rmtree(zone_slices)
+        
+        os.mkdir(zone_slices)            
+            
+        
+        with open(zone_coord_name,"r") as fcoord:
+            coordinates = fcoord.readline()
+            x0, y0, x1, y1 = [float(number) for number in coordinates.split()]
+
+            im_z = Image.open(zone_image_name)
+            im_z_np = np.array(im_z)
+
+            (w,h,d) = im_z_np.shape
+    
+            print('Slicing zone {}, please wait...'.format(zone_base))
+            print('\n')
+
+            dh = args.height
+            dw = args.width
+
+            n=0
+            for ih in range(0, dh*(h//dh)-dh, dh):
+                for iw in range(0, dw*(w//dw)-dw, dw):
+                    box = im_z_np[iw:iw+dw, ih:ih+dh, :]
+
+                    outfile = op.join(zone_slices,"{:05d}_w_{}_h_{}".format(n,iw+x0,ih+y0))
+                    np.save(outfile, box)
+                    n += 1
+
+                    
+    print('Finished!')
+
+
+if __name__ == '__main__':
+    main()

ClickCells/utils.py

+import glymur
+import PIL
+from PIL import Image
+import matplotlib.pyplot as plt
+
+fname = "/Users/saad/Desktop/FromJulia/mn96FS_c10_s4.jp2"
+fname = "/Users/saad/grot/mn38c23LY_Neu13_10d_L.tif"
+jp2 = glymur.Jp2k(fname)
+
+step = 2**3
+thumbnail = jp2[::step,::step]
+
+
+plt.show(thumbnail)
+
+
+# Convert previous Matlab database to new format
+import scipy.io as sio
+db = sio.loadmat("/Users/saad/data/Haber_CellCounting/celldb.mat",struct_as_record=False, squeeze_me=True)
+images = db['celldb'].images.data
+labels = -db['celldb'].images.label + 2
+
+import numpy as np
+images = data.transpose(3,0,1,2)
+
+outfile="/Users/saad/data/Haber_CellCounting/celldb.npz"
+shape = images.shape[1:]
+np.savez(outfile,counts=np.array(labels).flatten(),
+             images=np.array(images).reshape(-1,*shape))
+  

Models/sequential_cnn.py

+#!/usr/bin/env python3
+
+import numpy as np
+
+from __future__ import print_function
+import keras
+from keras.layers import Dense, Dropout, Flatten
+from keras.layers import Conv2D, MaxPooling2D
+from keras.models import Sequential
+
+model = Sequential()
+
+batch_size = 128
+num_classes = 2
+epochs = 12
+
+
+# input image dimensions
+img_rows, img_cols = 64, 64
+
+# Get data (use sklearn to separate train/test?)
+# dbname = "/Users/saad/grot/db.npz"
+dbname = "/Users/saad/data/Haber_CellCounting/celldb.npz"
+#dbname = "/Users/saad/Desktop/FromJulia/grotdb.npz"
+with np.load(dbname) as db:
+    images = db['images']
+    counts=db['counts']
+
+# Here make sure we have equal numbers of cell/nocell examplars?
+nwithcells = (counts>0).sum()
+nimages    = min(nwithcells, (counts==0).sum())
+
+indcells = np.where((counts>0))[0]
+np.random.shuffle(indcells)
+indnocells = np.where((counts==0))[0]
+np.random.shuffle(indnocells)
+
+rimages = np.concatenate((images[indcells[:nimages]],images[indnocells[:nimages]]), axis=0)
+rcounts = np.append(counts[indcells[:nimages]],counts[indnocells[:nimages]]) > 0
+        
+    
+from sklearn.model_selection import train_test_split
+x_train, x_test, y_train, y_test = train_test_split(rimages, rcounts, test_size=0.1)
+
+# convert class vectors to binary class matrices
+y_train = keras.utils.to_categorical(y_train, num_classes)
+y_test = keras.utils.to_categorical(y_test, num_classes)
+
+# CNN
+model = Sequential()
+
+input_shape = (img_rows, img_cols, 3)
+
+model.add(Conv2D(32, kernel_size=(3, 3),
+                 activation='relu',
+                 input_shape=input_shape))
+model.add(Conv2D(64, (3, 3),
+                 activation='relu'))
+model.add(MaxPooling2D(pool_size=(2, 2)))
+model.add(Dropout(0.25))
+model.add(Flatten())
+model.add(Dense(128, activation='relu'))
+model.add(Dropout(0.5))
+model.add(Dense(num_classes, activation='softmax'))
+
+model.compile(loss=keras.losses.categorical_crossentropy,
+              optimizer=keras.optimizers.Adadelta(),
+              metrics=['accuracy'])
+
+model.fit(x_train, y_train,
+          batch_size=batch_size,
+          epochs=epochs,
+          verbose=1,
+          validation_data=(x_test, y_test))
+
+score = model.evaluate(x_test, y_test, verbose=0)
+print('Test loss:', score[0])
+print('Test accuracy:', score[1])
+
+
+pred = model.predict(x_test)
+
+
+# Look at a few examples
+from random import randint
+i = randint(0,len(y_test)-1)
+plt.imshow(x_test[i])
+title = "Ground truth : {}\n".format("Cell!" if np.argmax(y_test[i]) == 1 else "No Cell :-(")
+title = title + "Prediction   : {}".format("Cell!" if np.argmax(pred[i]) == 1 else "No Cell :-(")
+plt.title(title)
+plt.show()
+
+    
+# Apply to entire image
+from PIL import Image
+Image.MAX_IMAGE_PIXELS = 5282426800
+
+image_name = "/Users/saad/grot/mn38c23LY_Neu13_10d_L.tif"
+im = np.asarray(Image.open(image_name))
+
+step = 40
+im_lr = im[::step,::Step,:]
+
+
+strided_im = np.lib.stride_tricks.as_strided(im_lr, (im_lr.shape[0] - 64, im_lr.shape[1] -64, 64, 64, 3),
+                                             (im_lr.strides[0], im_lr.strides[1], im_lr.strides[0], im_lr.strides[1], im_lr.strides[2]))
+
+strided_im = strided_im[:,64,64,3]