diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index d013aba..0000000 Binary files a/.DS_Store and /dev/null differ diff --git a/.gitignore b/.gitignore index e16642c..482c16f 100644 --- a/.gitignore +++ b/.gitignore @@ -117,3 +117,7 @@ venv.bak/ .idea/ .vs/ *.~lock.* + + +# macOS +*.DS_Store diff --git a/docs/searchindex.js b/docs/searchindex.js index a52c11b..8a80b12 100644 --- a/docs/searchindex.js +++ b/docs/searchindex.js @@ -1 +1 @@ -Search.setIndex({docnames:["index","main/dev/CONTRIBUTING","main/user_guide/tools/cells_to_brainrender"],envversion:{"sphinx.domains.c":1,"sphinx.domains.changeset":1,"sphinx.domains.citation":1,"sphinx.domains.cpp":1,"sphinx.domains.javascript":1,"sphinx.domains.math":2,"sphinx.domains.python":1,"sphinx.domains.rst":1,"sphinx.domains.std":1,"sphinx.ext.viewcode":1,sphinx:56},filenames:["index.rst","main/dev/CONTRIBUTING.md","main/user_guide/tools/cells_to_brainrender.md"],objects:{},objnames:{},objtypes:{},terms:{"20gb":[],"2gb":[],"case":1,"default":2,"function":1,"import":2,"new":1,"true":[1,2],"try":1,The:[],These:1,abl:[],activ:[],add:1,add_cells_from_fil:2,adding:1,advanc:[],all:[1,2],allen:[],allthough:[],also:[],altern:[],amap:[0,1],analysi:0,ani:1,anoth:1,appropri:[],approv:1,atla:2,atlas:[],automat:1,avail:[],avoid:[],bash:[],befor:1,being:[],below:[],black:1,brain:0,brainrend:0,branch:1,build:1,call:2,can:2,carri:[],cell:0,cellfind:[0,1],cellfinder_download:[],cells_in_standard_spac:2,centr:0,chang:[1,2],check:1,choos:[],clean:1,cli:1,clone:1,cluster:[],code:1,collect:0,com:1,command:1,commit:1,companion:1,compat:[],comput:[],conf:[],connect:[],consist:1,contribut:0,convert:0,correct:[],creat:[],cudatoolkit:[],current:[],data:[0,2],date:1,defin:2,design:1,dev:1,develop:1,dimens:2,disk:[],document:[],done:[],driver:[],easier:1,edit:1,either:1,elsewher:[],end:2,ensur:1,error:1,etc:[],everyth:[],exist:[],exit:1,exported_cel:2,extent:2,extract:1,fetch:1,file:[],filenam:2,first:2,fork:1,format:0,found:[],from:0,futur:[],git:1,github:1,has:2,have:[],hdf:2,henc:1,here:[],highli:[],histori:1,hook:1,how:[],http:1,identifi:2,imag:0,imlib:1,infer:[],initi:1,instal:1,instruct:[],internet:[],invers:1,isol:[],issu:1,jupyt:2,keep:1,kei:2,latest:[],learn:[],like:[],line:1,link:[],linux:1,machin:[],mai:1,main:1,maintain:1,make:1,master:1,max:2,maximum:2,merg:1,miniconda3:[],miniconda:[],minim:[],modifi:[],mous:[],must:2,name:[],necessari:1,need:[],neuro:1,onc:[],onli:[],origin:[],other:[],out:[],output:2,packag:1,particularli:0,path:[],pip:[0,1],pixel:2,pleas:1,points_to_brainrend:2,pre:1,prevent:1,previous:[],process:0,push:1,py37:1,pytest:1,python:[],recommend:[],refer:[],registr:[],relev:1,rememb:[],remot:1,render:2,repositori:1,run:[1,2],sainsburi:0,sainsburywellcomecentr:1,save:[],scene:2,second:2,see:2,set:1,should:[1,2],size:2,softwar:1,sourc:[],space:2,specifi:2,squash:1,still:[],style:1,submit:1,suppli:[],system:[],tag:1,target:1,tell:[],test:[],thi:1,third:2,time:[],tmp:[],tool:1,twice:[],ubuntu:[],unit_test:1,updat:[],upon:1,upstream:1,use:1,used:[],uses:1,using:[],version:1,via:1,want:[],wellcom:0,when:[],where:[],which:[],whole:0,wide:[],work:1,workflow:1,would:[],x86_64:[],xml:2,you:[],your:1,your_usernam:1},titles:["neuro","Contributing","Converting cells from cellfinder to brainrender format"],titleterms:{For:[],The:2,about:[],addit:[],also:2,argument:2,atla:[],brainrend:2,cell:2,cellfind:2,classif:[],conda:[],contribut:1,convert:2,cuda:[],cudnn:[],depend:1,develop:0,download:[],environ:[],file:2,follow:2,format:[1,2],from:2,gpu:[],guid:[],instal:0,introduct:0,mai:2,model:[],need:2,neuro:0,option:2,posit:2,pull:1,releas:1,request:1,requir:[],set:[],setup:1,support:[],test:1,thi:2,tool:0,train:[],travi:1,used:2,user:[],visualis:2}}) \ No newline at end of file +Search.setIndex({docnames:["index","main/dev/CONTRIBUTING","main/user_guide/tools/cells_to_brainrender"],envversion:{"sphinx.domains.c":1,"sphinx.domains.changeset":1,"sphinx.domains.citation":1,"sphinx.domains.cpp":1,"sphinx.domains.javascript":1,"sphinx.domains.math":2,"sphinx.domains.python":1,"sphinx.domains.rst":1,"sphinx.domains.std":1,"sphinx.ext.viewcode":1,sphinx:56},filenames:["index.rst","main/dev/CONTRIBUTING.md","main/user_guide/tools/cells_to_brainrender.md"],objects:{},objnames:{},objtypes:{},terms:{"case":1,"default":2,"function":1,"import":2,"new":1,"true":[1,2],"try":1,These:1,add:1,add_cells_from_fil:2,adding:1,all:[1,2],amap:[0,1],analysi:0,ani:1,anoth:1,approv:1,atla:2,automat:1,befor:1,black:1,brain:0,brainrend:0,branch:1,build:1,call:2,can:2,cell:0,cellfind:[0,1],cells_in_standard_spac:2,centr:0,chang:[1,2],check:1,clean:1,cli:1,clone:1,code:1,collect:0,com:1,command:1,commit:1,companion:1,consist:1,contribut:0,convert:0,data:[0,2],date:1,defin:2,design:1,dev:1,develop:1,dimens:2,easier:1,edit:1,either:1,end:2,ensur:1,error:1,exit:1,exported_cel:2,extent:2,extract:1,fetch:1,filenam:2,first:2,fork:1,format:0,from:0,git:1,github:1,has:2,hdf:2,henc:1,histori:1,hook:1,http:1,identifi:2,imag:0,imlib:1,initi:1,instal:1,invers:1,issu:1,jupyt:2,keep:1,kei:2,line:1,linux:1,mai:1,main:1,maintain:1,make:1,master:1,max:2,maximum:2,merg:1,must:2,necessari:1,neuro:1,output:2,packag:1,particularli:0,pip:[0,1],pixel:2,pleas:1,points_to_brainrend:2,pre:1,prevent:1,process:0,push:1,py37:1,pytest:1,relev:1,remot:1,render:2,repositori:1,run:[1,2],sainsburi:0,sainsburywellcomecentr:1,scene:2,second:2,see:2,set:1,should:[1,2],size:2,softwar:1,space:2,specifi:2,squash:1,style:1,submit:1,tag:1,target:1,thi:1,third:2,tool:1,unit_test:1,upon:1,upstream:1,use:1,uses:1,version:1,via:1,wellcom:0,whole:0,work:1,workflow:1,xml:2,your:1,your_usernam:1},titles:["neuro","Contributing","Converting cells from cellfinder to brainrender format"],titleterms:{The:2,also:2,argument:2,brainrend:2,cell:2,cellfind:2,contribut:1,convert:2,depend:1,develop:0,file:2,follow:2,format:[1,2],from:2,instal:0,introduct:0,mai:2,need:2,neuro:0,option:2,posit:2,pull:1,releas:1,request:1,setup:1,test:1,thi:2,tool:0,travi:1,used:2,visualis:2}}) \ No newline at end of file diff --git a/neuro/injection_finder/.DS_Store b/neuro/injection_finder/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/neuro/injection_finder/.DS_Store and /dev/null differ diff --git a/neuro/injection_finder/extraction.py b/neuro/injection_finder/extraction.py index c24aae7..0c8095b 100644 --- a/neuro/injection_finder/extraction.py +++ b/neuro/injection_finder/extraction.py @@ -1,20 +1,21 @@ import os import numpy as np +from pathlib import Path from skimage.filters import gaussian as gaussian_filter from skimage.filters import threshold_otsu from skimage import measure from brainio import brainio +from imlib.IO.surfaces import marching_cubes_to_obj +from imlib.image.orient import reorient_image -from registration import get_registered_image -from utils import reorient_image, marching_cubes_to_obj, get_largest_component -from parsers import extraction_parser +from neuro.injection_finder.registration import get_registered_image +from neuro.injection_finder.parsers import extraction_parser -# For logging +import neuro as package_for_log import logging from fancylog import fancylog -import fancylog as package class Extractor: @@ -22,7 +23,6 @@ def __init__( self, img_filepath, registration_folder, - logging, overwrite=False, gaussian_kernel=2, percentile_threshold=99.95, @@ -30,21 +30,25 @@ def __init__( obj_path=None, overwrite_registration=False, ): - """ - Extractor processes a downsampled.nii image to extract the location of the injection site. - This is done by registering the image to the allen CCF, blurring, thresholding and finally a - marching cube algorithm to extract the surface of the injection site. - - :param img_filepath: str, path to .nii file - :param registration_folder: str, path to the registration folder [from cellfinder or amap] - :param logging: instance of fancylog logger - :param overwrite: bool, if False it will avoid overwriting files - :gaussian_kernel: float, size of kernel used for smoothing - :param percentile_threshold: float, in range [0, 1] percentile to use for thresholding - :param threshold_type: str, either ['otsu', 'percentile'], type of threshold used - :param obj_path: path to .obj file destination. - :param overwrite_registration: if false doesn't overwrite the registration step + Extractor processes a downsampled.nii image to extract the location of + the injection site. + This is done by registering the image to the allen CCF, blurring, + thresholding and finally a marching cube algorithm to extract the + surface of the injection site. + + :param img_filepath: str, path to .nii file + :param registration_folder: str, path to the registration folder + [from cellfinder or amap] + :param overwrite: bool, if False it will avoid overwriting files + :gaussian_kernel: float, size of kernel used for smoothing + :param percentile_threshold: float, in range [0, 1] percentile to use + for thresholding + :param threshold_type: str, either ['otsu', 'percentile'], + type of threshold used + :param obj_path: path to .obj file destination. + :param overwrite_registration: if false doesn't overwrite the + registration step """ # Get arguments @@ -72,14 +76,14 @@ def setup(self): self.img_filepath.split(".")[0] + "_thresholded.nii" ) - # Get path to obj file and check if it existsts + # Get path to obj file and check if it exists if self.obj_path is None: self.obj_path = self.img_filepath.split(".")[0] + ".obj" if os.path.isfile(self.obj_path) and not self.overwrite: self.logging.warning( - "A file exists already at {}. \ - Analysis will not run as overwrite is set disabled".format( + "A file exists already at {}." + "Analysis will not run as overwrite is set disabled".format( self.obj_path ) ) @@ -88,7 +92,6 @@ def setup(self): image = get_registered_image( self.img_filepath, self.registration_folder, - self.logging, overwrite=self.overwrite_registration, ) return image @@ -103,12 +106,12 @@ def extract(self, image, voxel_size=10): # Gaussian filter kernel_shape = [self.gaussian_kernel, self.gaussian_kernel, 6] - filtered = gaussian_filter(image, kernel_shape) + image = gaussian_filter(image, kernel_shape) self.logging.info("Filtering completed") # Thresholding if self.threshold_type.lower() == "otsu": - thresh = threshold_otsu(filtered) + thresh = threshold_otsu(image) self.logging.info( "Thresholding with {} threshold type".format( self.threshold_type @@ -119,21 +122,20 @@ def extract(self, image, voxel_size=10): self.threshold_type.lower() == "percentile" or self.threshold_type.lower() == "perc" ): - thresh = np.percentile(filtered.ravel(), self.percentile_threshold) + thresh = np.percentile(image.ravel(), self.percentile_threshold) self.logging.info( - "Thresholding with {} threshold type. {}th percentile [{}]".format( + "Thresholding with {} threshold type. " + "{}th percentile [{}]".format( self.threshold_type, self.percentile_threshold, thresh ) ) else: - raise valueError( + raise ValueError( "Unrecognised thresholding type: " + self.threshold_type ) - binary = filtered > thresh - oriented_binary = reorient_image( - binary, invert_axes=[2,], orientation="coronal" - ) + binary = image > thresh + binary = keep_n_largest_objects(binary) # Save thresholded image if not os.path.isfile(self.thresholded_savepath) or self.overwrite: @@ -144,10 +146,14 @@ def extract(self, image, voxel_size=10): ) brainio.to_nii(binary.astype(np.int16), self.thresholded_savepath) + binary = reorient_image( + binary, invert_axes=[2,], orientation="coronal" + ) + # apply marching cubes self.logging.info("Extracting surface from thresholded image") verts, faces, normals, values = measure.marching_cubes_lewiner( - oriented_binary, 0, step_size=1 + binary, 0, step_size=1 ) # Scale to atlas spacing @@ -159,8 +165,42 @@ def extract(self, image, voxel_size=10): faces = faces + 1 marching_cubes_to_obj((verts, faces, normals, values), self.obj_path) - # Keep only the largest connected component - get_largest_component(self.obj_path) + +def keep_n_largest_objects(numpy_array, n=1, connectivity=None): + """ + Given an input binary numpy array, return a "clean" array with only the + n largest connected components remaining + + Inspired by stackoverflow.com/questions/47540926 + + TODO: optimise + + :param numpy_array: Binary numpy array + :param n: How many objects to keep + :param connectivity: Labelling connectivity (see skimage.measure.label) + :return: "Clean" numpy array with n largest objects + """ + + labels = measure.label(numpy_array, connectivity=connectivity) + assert labels.max() != 0 # assume at least 1 CC + n_largest_objects = get_largest_non_zero_object(labels) + if n > 1: + i = 1 + while i < n: + labels[n_largest_objects] = 0 + n_largest_objects += get_largest_non_zero_object(labels) + i += 1 + return n_largest_objects + + +def get_largest_non_zero_object(label_image): + """ + In a labelled (each object assigned an int) numpy array. Return the + largest object with a value >= 1. + :param label_image: Output of skimage.measure.label + :return: Boolean numpy array or largest object + """ + return label_image == np.argmax(np.bincount(label_image.flat)[1:]) + 1 def main(): @@ -168,23 +208,30 @@ def main(): # Get output directory if args.output_directory is None: - outdir = os.get_cwd() + outdir = os.getcwd() elif not os.path.isdir(args.output_directory): raise ValueError("Output directory invalid") else: outdir = args.output_directory + if args.obj_path is None: + args.obj_path = Path(args.img_filepath).with_suffix(".obj") + else: + args.obj_path = Path(args.obj_path) + # Start log - log_name = "injection_finder_{}".format( - os.path.split(args.registration_folder)[-1] + fancylog.start_logging( + outdir, + package_for_log, + filename="injection_finder", + verbose=args.debug, + log_to_file=args.save_log, ) - fancylog.start_logging(outdir, package, filename=log_name, verbose=True) # Start extraction Extractor( args.img_filepath, args.registration_folder, - logging, overwrite=args.overwrite, gaussian_kernel=args.gaussian_kernel, percentile_threshold=args.percentile_threshold, diff --git a/neuro/injection_finder/parsers.py b/neuro/injection_finder/parsers.py index 3b9af16..5532e1e 100644 --- a/neuro/injection_finder/parsers.py +++ b/neuro/injection_finder/parsers.py @@ -38,8 +38,8 @@ def extraction_parser(): "--obj-path", dest="obj_path", type=str, - default=False, - help="Path to output .obj file. Optional.", + default=None, + help="Path to output .obj file. Will default to the image directory.", ) parser.add_argument( @@ -57,7 +57,8 @@ def extraction_parser(): dest="percentile_threshold", type=float, default=99.995, - help="Float in range [0, 100]. The percentile number of pixel intensity values for tresholding", + help="Float in range [0, 100]. The percentile number of pixel " + "intensity values for tresholding", ) parser.add_argument( @@ -66,7 +67,8 @@ def extraction_parser(): dest="threshold_type", type=str, default="otsu", - help="'otsu' or 'percentile'. Determines how the threshold value is computed", + help="'otsu' or 'percentile'. Determines how the threshold " + "value is computed", ) parser.add_argument( @@ -77,4 +79,17 @@ def extraction_parser(): default="False", help="If false skip running again the registration", ) + parser.add_argument( + "--debug", + dest="debug", + action="store_true", + help="Debug mode. Will increase verbosity of logging and save all " + "intermediate files for diagnosis of software issues.", + ) + parser.add_argument( + "--save-log", + dest="save_log", + action="store_true", + help="Save logging to file (in addition to logging to terminal).", + ) return parser diff --git a/neuro/injection_finder/registration.py b/neuro/injection_finder/registration.py index a979a55..3abe057 100644 --- a/neuro/injection_finder/registration.py +++ b/neuro/injection_finder/registration.py @@ -1,43 +1,15 @@ import os +import logging from brainio import brainio -from amap.tools import source_files -from amap.config.config import get_binary -from amap.tools.exceptions import RegistrationError +from imlib.register.niftyreg.transform import run_transform -from cellfinder.tools.system import ( - safe_execute_command, - SafeExecuteCommandError, -) - -PROGRAM_NAME = "reg_resample" DEFAULT_CONTROL_POINT_FILE = "inverse_control_point_file.nii" default_atlas_name = "brain_filtered.nii" -def prepare_segmentation_cmd( - program_path, - floating_image_path, - output_file_name, - destination_image_filename, - control_point_file, -): - cmd = "{} -cpp {} -flo {} -ref {} -res {}".format( - program_path, - control_point_file, - floating_image_path, - destination_image_filename, - output_file_name, - ) - return cmd - - -def get_registered_image(nii_path, registration_dir, logging, overwrite=False): - # get binaries - nifty_reg_binaries_folder = source_files.get_niftyreg_binaries() - program_path = get_binary(nifty_reg_binaries_folder, PROGRAM_NAME) - +def get_registered_image(nii_path, registration_dir, overwrite=False): # get file paths basedir = os.path.split(nii_path)[0] output_filename = os.path.join( @@ -45,11 +17,9 @@ def get_registered_image(nii_path, registration_dir, logging, overwrite=False): "{}_transformed.nii".format(os.path.split(nii_path)[1].split(".")[0]), ) if os.path.isfile(output_filename) and not overwrite: - run = False + logging.info("Skipping registration as output file already exists") else: - run = True - if run: destination_image = os.path.join(registration_dir, default_atlas_name) control_point_file = os.path.join( registration_dir, DEFAULT_CONTROL_POINT_FILE @@ -58,19 +28,14 @@ def get_registered_image(nii_path, registration_dir, logging, overwrite=False): log_file_path = os.path.join(basedir, "registration_log.txt") error_file_path = os.path.join(basedir, "registration_err.txt") - reg_cmd = prepare_segmentation_cmd( - program_path, + logging.info("Running registration") + run_transform( nii_path, output_filename, destination_image, control_point_file, + log_file_path, + error_file_path, ) - logging.info("Running registration") - try: - safe_execute_command(reg_cmd, log_file_path, error_file_path) - except SafeExecuteCommandError as err: - raise RegistrationError("Registration failed; {}".format(err)) - else: - logging.info("Skipping registration as output file already exists") return brainio.load_any(output_filename) diff --git a/neuro/injection_finder/utils.py b/neuro/injection_finder/utils.py deleted file mode 100644 index 9eee265..0000000 --- a/neuro/injection_finder/utils.py +++ /dev/null @@ -1,68 +0,0 @@ -import numpy as np -from vtkplotter.analysis import extractLargestRegion - -# ----------------------- IMAGE MANIPULATION FUNCTIONS ----------------------- # -def reorient_image(image, invert_axes=None, orientation="saggital"): - """ - Reorients the image to the coordinate space of the atlas - - :param image_path: str - :param threshold: float - :param invert_axes: tuple (Default value = None) - :param image: - :param orientation: (Default value = "saggital") - - """ - if invert_axes is not None: - for axis in list(invert_axes): - image = np.flip(image, axis=axis) - - if orientation is not "saggital": - if orientation is "coronal": - transposition = (2, 1, 0) - elif orientation is "horizontal": - transposition = (1, 2, 0) - - image = np.transpose(image, transposition) - return image - - -# ------------------------- MARCHING CUBES FUNCTIONS ------------------------- # -def marching_cubes_to_obj(marching_cubes_out, output_file): - """ - Saves the output of skimage.measure.marching_cubes as an .obj file - - :param marching_cubes_out: tuple - :param output_file: str - - """ - - verts, faces, normals, _ = marching_cubes_out - with open(output_file, "w") as f: - for item in verts: - f.write(f"v {item[0]} {item[1]} {item[2]}\n") - for item in normals: - f.write(f"vn {item[0]} {item[1]} {item[2]}\n") - for item in faces: - f.write( - f"f {item[0]}//{item[0]} {item[1]}//{item[1]} " - f"{item[2]}//{item[2]}\n" - ) - f.close() - - -# ----------------------------- VTKPLOTTER UTILS ----------------------------- # -def get_center_of_mass(actor): - """ - Get the center of mass of a vtk actor - """ - return actor.centerOfMass() - -def get_largest_component(obj_filepath): - """ - Given a .obj file with multiple disconnected meshes in it, it - selects the largest of these and discards the rest. - """ - actor = load(obj_filepath) - actor = extractLargestRegion(actor).flipNormals() - save(actor, obj_filepath)