#NEWLINE# def vertebral_detection( fname , fname_seg , init_disc , verbose , laplacian = 0 ) : #NEWLINE# #INDENT# shift_AP = 32 #NEWLINE# size_AP = 11 #NEWLINE# size_RL = 1 #NEWLINE# size_IS = 19 #NEWLINE# smooth_factor = [ 9 , 3 , 1 ] #NEWLINE# thr_corr = 0.2 #NEWLINE# fig_anat_straight = 1 #NEWLINE# fig_pattern = 2 #NEWLINE# contrast_template = 't2' #NEWLINE# if ( contrast_template == 't1' ) : #NEWLINE# #INDENT# contrast_template = 'T1' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( contrast_template == 't2' ) : #NEWLINE# #INDENT# contrast_template = 'T2' #NEWLINE# #UNINDENT# #UNINDENT# from os import path #NEWLINE# path_script = path . dirname ( __file__ ) #NEWLINE# path_sct = slash_at_the_end ( path . dirname ( path_script ) , 1 ) #NEWLINE# folder_template = 'data/template/' #NEWLINE# fname_template_list = glob ( ( ( ( ( path_sct + folder_template ) + '*' ) + contrast_template ) + '.nii.gz' ) ) #NEWLINE# fname_template = list_distance [ 0 ] #NEWLINE# fname_disc_list = glob ( ( ( path_sct + folder_template ) + '*_disc.nii.gz' ) ) #NEWLINE# fname_disc = fname_disc_list [ 0 ] #NEWLINE# data_template = Image ( fname_template ) . data #NEWLINE# data_disc_template = Image ( fname_disc ) . data #NEWLINE# if laplacian : #NEWLINE# #INDENT# printv ( '\nApplying Laplacian filter to template data...' , verbose ) #NEWLINE# from sct_maths import laplacian #NEWLINE# data_template = laplacian ( data_template . astype ( float ) , [ 1 , 1 , 1 ] ) #NEWLINE# #UNINDENT# list_disc_value_template = sorted ( data_disc_template [ data_disc_template . nonzero ( ) ] ) #NEWLINE# list_disc_z_template = [ int ( np . where ( ( data_disc_template == list_disc_value_template [ i ] ) ) [ 2 ] ) for i in range ( len ( list_disc_value_template ) ) ] #NEWLINE# printv ( ( '\nDisc values from template: ' + str ( list_disc_value_template ) ) , verbose ) #NEWLINE# list_distance_template = ( np . diff ( list_disc_z_template ) * ( - 1 ) ) . tolist ( ) #NEWLINE# printv ( ( 'Distances between discs (in voxel): ' + str ( list_distance_template ) ) , verbose ) #NEWLINE# if ( verbose == 2 ) : #NEWLINE# #INDENT# import matplotlib.pyplot as plt #NEWLINE# plt . ion ( ) #NEWLINE# #UNINDENT# img = Image ( fname ) #NEWLINE# data = img . data #NEWLINE# from scipy.ndimage.filters import gaussian_filter #NEWLINE# data = gaussian_filter ( data , smooth_factor , output = None , mode = 'reflect' ) #NEWLINE# ( nx , ny , nz , nt , px , py , pz , pt ) = img . dim #NEWLINE# z = range ( nz ) #NEWLINE# xc = int ( round ( ( nx / 2 ) ) ) #NEWLINE# yc = int ( round ( ( ny / 2 ) ) ) #NEWLINE# if ( verbose == 2 ) : #NEWLINE# #INDENT# plt . matshow ( np . mean ( data [ ( xc - size_RL ) : ( xc + size_RL ) , : , : ] , axis = 0 ) . transpose ( ) , fignum = fig_anat_straight , cmap = plt . cm . gray , origin = 'lower' ) #NEWLINE# plt . title ( 'Anatomical image' ) #NEWLINE# plt . autoscale ( enable = False ) #NEWLINE# #UNINDENT# printv ( '\nDetect intervertebral discs...' , verbose ) #NEWLINE# current_z = init_disc [ 0 ] #NEWLINE# current_disc = init_disc [ 1 ] #NEWLINE# list_disc_z = [ ] #NEWLINE# list_disc_value = [ ] #NEWLINE# direction = 'superior' #NEWLINE# search_next_disc = True #NEWLINE# while search_next_disc : #NEWLINE# #INDENT# printv ( ( ( ( ( ( 'Current disc: ' + str ( current_disc ) ) + ' (z=' ) + str ( current_z ) ) + '). Direction: ' ) + direction ) , verbose ) #NEWLINE# try : #NEWLINE# #INDENT# current_z_template = int ( np . where ( ( data_disc_template == current_disc ) ) [ 2 ] ) #NEWLINE# #UNINDENT# except TypeError : #NEWLINE# #INDENT# printv ( 'WARNING: Reached the bottom of the template. Stop searching.' , verbose , 'warning' ) #NEWLINE# break #NEWLINE# #UNINDENT# pattern = data_template [ ( xc - size_RL ) : ( ( xc + size_RL ) + 1 ) , ( ( yc + shift_AP ) - size_AP ) : ( ( ( yc + shift_AP ) + size_AP ) + 1 ) , ( current_z_template - size_IS ) : ( ( current_z_template + size_IS ) + 1 ) ] #NEWLINE# pattern1d = pattern . ravel ( ) #NEWLINE# if ( verbose == 2 ) : #NEWLINE# #INDENT# plt . figure ( fig_anat_straight ) #NEWLINE# plt . scatter ( ( yc + shift_AP ) , current_z , c = 'red' , s = 50 ) #NEWLINE# #UNINDENT# range_z = range ( ( - 10 ) , 10 ) #NEWLINE# I_corr = np . zeros ( len ( range_z ) ) #NEWLINE# allzeros = 0 #NEWLINE# ind_I = 0 #NEWLINE# for iz in range_z : #NEWLINE# #INDENT# if ( ( ( current_z + iz ) + size_IS ) > nz ) : #NEWLINE# #INDENT# padding_size = ( ( current_z + iz ) + size_IS ) #NEWLINE# data_chunk3d = data [ ( xc - size_RL ) : ( ( xc + size_RL ) + 1 ) , ( ( yc + shift_AP ) - size_AP ) : ( ( ( yc + shift_AP ) + size_AP ) + 1 ) , ( ( current_z + iz ) - size_IS ) : ( ( ( ( current_z + iz ) + size_IS ) + 1 ) - padding_size ) ] #NEWLINE# data_chunk3d = np . pad ( data_chunk3d , ( ( 0 , 0 ) , ( 0 , 0 ) , ( 0 , padding_size ) ) , 'constant' , constant_values = 0 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( ( current_z - iz ) - size_IS ) < 0 ) : #NEWLINE# #INDENT# padding_size = abs ( ( ( current_z - iz ) - size_IS ) ) #NEWLINE# data_chunk3d = data [ ( xc - size_RL ) : ( ( xc + size_RL ) + 1 ) , ( ( yc + shift_AP ) - size_AP ) : ( ( ( yc + shift_AP ) + size_AP ) + 1 ) , ( ( ( current_z - iz ) - size_IS ) + padding_size ) : ( ( ( current_z - iz ) + size_IS ) + 1 ) ] #NEWLINE# data_chunk3d = np . pad ( data_chunk3d , ( ( 0 , 0 ) , ( 0 , 0 ) , ( padding_size , 0 ) ) , 'constant' , constant_values = 0 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# data_chunk3d = data [ ( xc - size_RL ) : ( ( xc + size_RL ) + 1 ) , ( ( yc + shift_AP ) - size_AP ) : ( ( ( yc + shift_AP ) + size_AP ) + 1 ) , ( ( current_z + iz ) - size_IS ) : ( ( ( current_z + iz ) + size_IS ) + 1 ) ] #NEWLINE# #UNINDENT# #UNINDENT# data_chunk1d = data_chunk3d . ravel ( ) #NEWLINE# if ( ( data_chunk1d . size == pattern1d . size ) and np . any ( data_chunk1d ) ) : #NEWLINE# #INDENT# I_corr [ ind_I ] = calc_MI ( data_chunk1d , pattern1d , 32 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# allzeros = 1 #NEWLINE# #UNINDENT# ind_I = ( ind_I + 1 ) #NEWLINE# #UNINDENT# if allzeros : #NEWLINE# #INDENT# printv ( '.. WARNING: Data contained zero. We probably hit the edge of the image.' , verbose ) #NEWLINE# #UNINDENT# I_corr_adj = I_corr #NEWLINE# if np . any ( I_corr_adj ) : #NEWLINE# #INDENT# ind_peak = [ i for i in range ( len ( I_corr_adj ) ) if ( I_corr_adj [ i ] == max ( I_corr_adj ) ) ] [ 0 ] #NEWLINE# printv ( ( ( ( ( '.. Peak found: z=' + str ( ind_peak ) ) + ' (correlation = ' ) + str ( I_corr_adj [ ind_peak ] ) ) + ')' ) , verbose ) #NEWLINE# if ( I_corr_adj [ ind_peak ] < thr_corr ) : #NEWLINE# #INDENT# printv ( '.. WARNING: Correlation is too low. Using adjusted template distance.' , verbose ) #NEWLINE# ind_peak = range_z . index ( 0 ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# printv ( '.. WARNING: Correlation vector only contains zeros. Using adjusted template distance.' , verbose ) #NEWLINE# ind_peak = range_z . index ( 0 ) #NEWLINE# #UNINDENT# if ( verbose == 2 ) : #NEWLINE# #INDENT# plt . figure ( fig_pattern , figsize = ( 20 , 7 ) ) #NEWLINE# plt . subplot ( 141 ) #NEWLINE# plt . imshow ( np . flipud ( np . mean ( pattern [ : , : , : ] , axis = 0 ) . transpose ( ) ) , origin = 'upper' , cmap = plt . cm . gray , interpolation = 'none' ) #NEWLINE# plt . title ( 'Template pattern' ) #NEWLINE# plt . subplot ( 142 ) #NEWLINE# iz = 0 #NEWLINE# data_chunk3d = data [ ( xc - size_RL ) : ( ( xc + size_RL ) + 1 ) , ( ( yc + shift_AP ) - size_AP ) : ( ( ( yc + shift_AP ) + size_AP ) + 1 ) , ( ( current_z + iz ) - size_IS ) : ( ( ( current_z + iz ) + size_IS ) + 1 ) ] #NEWLINE# plt . imshow ( np . flipud ( np . mean ( data_chunk3d [ : , : , : ] , axis = 0 ) . transpose ( ) ) , origin = 'upper' , cmap = plt . cm . gray , interpolation = 'none' ) #NEWLINE# plt . title ( 'Subject at iz=0' ) #NEWLINE# plt . subplot ( 143 ) #NEWLINE# iz = range_z [ ind_peak ] #NEWLINE# data_chunk3d = data [ ( xc - size_RL ) : ( ( xc + size_RL ) + 1 ) , ( ( yc + shift_AP ) - size_AP ) : ( ( ( yc + shift_AP ) + size_AP ) + 1 ) , ( ( current_z + iz ) - size_IS ) : ( ( ( current_z + iz ) + size_IS ) + 1 ) ] #NEWLINE# plt . imshow ( np . flipud ( np . mean ( data_chunk3d [ : , : , : ] , axis = 0 ) . transpose ( ) ) , origin = 'upper' , cmap = plt . cm . gray , interpolation = 'none' ) #NEWLINE# plt . title ( ( 'Subject at iz=' + str ( iz ) ) ) #NEWLINE# plt . subplot ( 144 ) #NEWLINE# plt . plot ( I_corr_adj ) #NEWLINE# plt . title ( 'MI between template and subject pattern' ) #NEWLINE# ( plt . plot ( ind_peak , I_corr_adj [ ind_peak ] , 'ro' ) , plt . draw ( ) ) #NEWLINE# plt . axvline ( x = range_z . index ( 0 ) , linewidth = 1 , color = 'black' , linestyle = 'dashed' ) #NEWLINE# plt . axhline ( y = thr_corr , linewidth = 1 , color = 'r' , linestyle = 'dashed' ) #NEWLINE# ( plt . figure ( fig_pattern ) , plt . savefig ( ( ( '../fig_pattern_MI' + str ( current_disc ) ) + '.png' ) ) , plt . close ( ) ) #NEWLINE# #UNINDENT# current_z = ( current_z + range_z [ ind_peak ] ) #NEWLINE# if ( verbose == 2 ) : #NEWLINE# #INDENT# ( plt . figure ( fig_anat_straight ) , plt . scatter ( ( yc + shift_AP ) , current_z , c = 'yellow' , s = 50 ) ) #NEWLINE# ( plt . text ( ( ( yc + shift_AP ) + 4 ) , current_z , ( ( str ( current_disc ) + '/' ) + str ( ( current_disc + 1 ) ) ) , verticalalignment = 'center' , horizontalalignment = 'left' , color = 'yellow' , fontsize = 15 ) , plt . draw ( ) ) #NEWLINE# #UNINDENT# if ( direction == 'superior' ) : #NEWLINE# #INDENT# list_disc_z . insert ( 0 , current_z ) #NEWLINE# list_disc_value . insert ( 0 , current_disc ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( direction == 'inferior' ) : #NEWLINE# #INDENT# list_disc_z . append ( current_z ) #NEWLINE# list_disc_value . append ( current_disc ) #NEWLINE# #UNINDENT# #UNINDENT# if ( len ( list_disc_z ) > 1 ) : #NEWLINE# #INDENT# list_distance_current = ( np . diff ( list_disc_z ) * ( - 1 ) ) . tolist ( ) #NEWLINE# index_disc_identified = [ i for ( i , j ) in enumerate ( list_disc_value_template ) if ( j in list_disc_value [ : ( - 1 ) ] ) ] #NEWLINE# list_distance_template_identified = [ list_distance_template [ i ] for i in index_disc_identified ] #NEWLINE# list_subject_to_template_distance = [ ( float ( list_distance_current [ i ] ) / list_distance_template_identified [ i ] ) for i in range ( len ( list_distance_current ) ) ] #NEWLINE# correcting_factor = np . mean ( list_subject_to_template_distance ) #NEWLINE# printv ( ( '.. correcting factor: ' + str ( correcting_factor ) ) , verbose ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# correcting_factor = 1 #NEWLINE# #UNINDENT# list_distance = [ int ( round ( ( list_distance_template [ i ] * correcting_factor ) ) ) for i in range ( len ( list_distance_template ) ) ] #NEWLINE# if ( direction == 'superior' ) : #NEWLINE# #INDENT# try : #NEWLINE# #INDENT# approx_distance_to_next_disc = list_distance [ list_disc_value_template . index ( ( current_disc - 1 ) ) ] #NEWLINE# #UNINDENT# except ValueError : #NEWLINE# #INDENT# printv ( ( 'WARNING: Disc value not included in template. Using previously-calculated distance: ' + str ( approx_distance_to_next_disc ) ) ) #NEWLINE# #UNINDENT# current_z = ( current_z + approx_distance_to_next_disc ) #NEWLINE# current_disc = ( current_disc - 1 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( direction == 'inferior' ) : #NEWLINE# #INDENT# try : #NEWLINE# #INDENT# approx_distance_to_next_disc = list_distance [ list_disc_value_template . index ( current_disc ) ] #NEWLINE# #UNINDENT# except : #NEWLINE# #INDENT# printv ( ( 'WARNING: Disc value not included in template. Using previously-calculated distance: ' + str ( approx_distance_to_next_disc ) ) ) #NEWLINE# #UNINDENT# current_z = ( current_z - approx_distance_to_next_disc ) #NEWLINE# current_disc = ( current_disc + 1 ) #NEWLINE# #UNINDENT# #UNINDENT# if ( ( current_z >= nz ) or ( current_disc == 1 ) ) : #NEWLINE# #INDENT# printv ( '.. Switching to inferior direction.' , verbose ) #NEWLINE# direction = 'inferior' #NEWLINE# current_disc = ( init_disc [ 1 ] + 1 ) #NEWLINE# current_z = ( init_disc [ 0 ] - list_distance [ list_disc_value_template . index ( current_disc ) ] ) #NEWLINE# #UNINDENT# if ( current_z <= 0 ) : #NEWLINE# #INDENT# search_next_disc = False #NEWLINE# #UNINDENT# #UNINDENT# upper_disc = min ( list_disc_value ) #NEWLINE# if ( not ( upper_disc == 1 ) ) : #NEWLINE# #INDENT# printv ( ( 'Adding top disc based on adjusted template distance: #' + str ( ( upper_disc - 1 ) ) ) , verbose ) #NEWLINE# approx_distance_to_next_disc = list_distance [ list_disc_value_template . index ( ( upper_disc - 1 ) ) ] #NEWLINE# next_z = ( max ( list_disc_z ) + approx_distance_to_next_disc ) #NEWLINE# printv ( ( '.. approximate distance: ' + str ( approx_distance_to_next_disc ) ) , verbose ) #NEWLINE# if ( next_z > nz ) : #NEWLINE# #INDENT# list_disc_z . insert ( 0 , nz ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# list_disc_z . insert ( 0 , next_z ) #NEWLINE# #UNINDENT# list_disc_value . insert ( 0 , ( upper_disc - 1 ) ) #NEWLINE# #UNINDENT# seg = Image ( fname_seg ) #NEWLINE# for iz in range ( nz ) : #NEWLINE# #INDENT# try : #NEWLINE# #INDENT# ind_above_iz = max ( [ i for i in range ( len ( list_disc_z ) ) if ( list_disc_z [ i ] > iz ) ] ) #NEWLINE# #UNINDENT# except ValueError : #NEWLINE# #INDENT# vertebral_level = 0 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# vertebral_level = ( list_disc_value [ ind_above_iz ] + 1 ) #NEWLINE# #UNINDENT# ind_nonzero = np . nonzero ( seg . data [ : , : , iz ] ) #NEWLINE# seg . data [ ( ind_nonzero [ 0 ] , ind_nonzero [ 1 ] , iz ) ] = vertebral_level #NEWLINE# if ( verbose == 2 ) : #NEWLINE# #INDENT# plt . figure ( fig_anat_straight ) #NEWLINE# plt . scatter ( int ( round ( ( ny / 2 ) ) ) , iz , c = vertebral_level , vmin = min ( list_disc_value ) , vmax = max ( list_disc_value ) , cmap = 'prism' , marker = '_' , s = 200 ) #NEWLINE# #UNINDENT# #UNINDENT# seg . file_name += '_labeled' #NEWLINE# seg . save ( ) #NEWLINE# if ( verbose == 2 ) : #NEWLINE# #INDENT# ( plt . figure ( fig_anat_straight ) , plt . savefig ( '../fig_anat_straight_with_labels.png' ) ) #NEWLINE# plt . close ( )