<s> import os <EOL> import getopt <EOL> import commands <EOL> import math <EOL> import sys <EOL> import scipy <EOL> import scipy . signal <EOL> import scipy . fftpack <EOL> import pylab as pl <EOL> import sct_utils as sct <EOL> from sct_utils import fsloutput <EOL> try : <EOL> import nibabel <EOL> except ImportError : <EOL> print '<STR_LIT>' <EOL> sys . exit ( <NUM_LIT:2> ) <EOL> try : <EOL> import numpy as np <EOL> except ImportError : <EOL> print '<STR_LIT>' <EOL> sys . exit ( <NUM_LIT:2> ) <EOL> class label_class : <EOL> def __init__ ( self , contrast ) : <EOL> self . input_path = '<STR_LIT>' <EOL> self . input_anat = '<STR_LIT>' <EOL> self . output_path = '<STR_LIT>' <EOL> self . output_labled_centerline = '<STR_LIT>' <EOL> self . output_labled_surface = '<STR_LIT>' <EOL> self . surface_do = <NUM_LIT:1> <EOL> self . segmentation_do = <NUM_LIT:0> <EOL> self . segmentation_interval = <NUM_LIT:30> <EOL> self . segmentation_nom_radius = <NUM_LIT:5> <EOL> self . segmentation_tolerance = <NUM_LIT> <EOL> self . segmentation_ratio_criteria = <NUM_LIT> <EOL> self . segmentation_num_angles = <NUM_LIT:64> <EOL> self . segmentation_update_multiplier = <NUM_LIT> <EOL> self . segmentation_shear_force_multiplier = <NUM_LIT:0.5> <EOL> self . segmentation_max_coeff_horizontal = <NUM_LIT:10> <EOL> self . segmentation_max_coeff_vertical = <NUM_LIT:10> <EOL> self . segmentation_centerline = '<STR_LIT>' <EOL> self . segmentation_surface = '<STR_LIT>' <EOL> self . segmentation_straightened = '<STR_LIT>' <EOL> self . log = '<STR_LIT>' <EOL> self . input_centerline = '<STR_LIT>' <EOL> self . input_surface = '<STR_LIT>' <EOL> self . shift_AP = <NUM_LIT> <EOL> self . size_AP = <NUM_LIT:6> <EOL> self . size_RL = <NUM_LIT:5> <EOL> self . verbose = <NUM_LIT:1> <EOL> if contrast == '<STR_LIT>' : <EOL> self . segmentation_image_type = <NUM_LIT:1> <EOL> else : <EOL> self . segmentation_image_type = <NUM_LIT:2> <EOL> def main ( ) : <EOL> contrast = '<STR_LIT>' <EOL> label = label_class ( contrast ) <EOL> try : <EOL> opts , args = getopt . getopt ( sys . argv [ <NUM_LIT:1> : ] , '<STR_LIT>' ) <EOL> except getopt . GetoptError as err : <EOL> print str ( err ) <EOL> usage ( ) <EOL> for opt , arg in opts : <EOL> if opt == '<STR_LIT>' : <EOL> usage ( ) <EOL> elif opt in ( '<STR_LIT>' ) : <EOL> label . input_anat = arg <EOL> elif opt in ( '<STR_LIT:-c>' ) : <EOL> label . output_labled_centerline = arg <EOL> elif opt in ( '<STR_LIT>' ) : <EOL> label . output_labled_surface = arg <EOL> if label . input_anat == '<STR_LIT>' or label . output_labled_centerline == '<STR_LIT>' : <EOL> print '<STR_LIT>' <EOL> usage ( ) <EOL> input_anat = label . input_path + label . input_anat + '<STR_LIT>' <EOL> if label . segmentation_do == <NUM_LIT:1> : <EOL> input_centerline = label . output_path + label . segmentation_centerline + '<STR_LIT>' <EOL> input_surface = label . output_path + label . segmentation_surface + '<STR_LIT>' <EOL> else : <EOL> input_centerline = label . input_path + label . input_centerline + '<STR_LIT>' <EOL> input_surface = label . input_path + label . input_surface + '<STR_LIT>' <EOL> output_centerline_vertebra = label . output_path + label . output_labled_centerline <EOL> output_surface_vertebra = label . output_path + label . output_labled_surface <EOL> surface_do = label . surface_do <EOL> sct . check_file_exist ( input_anat ) <EOL> if contrast == '<STR_LIT>' : <EOL> labeling_vertebrae_T1 ( label , input_anat , input_centerline , input_surface , output_centerline_vertebra , output_surface_vertebra , surface_do ) <EOL> else : <EOL> labeling_vertebrae_T2 ( label , input_anat , input_centerline , input_surface , output_centerline_vertebra , output_surface_vertebra , surface_do ) <EOL> def labeling_vertebrae_T1 ( label , input_anat , input_centerline , input_surface , output_centerline_vertebra , output_surface_vertebra , surface_do ) : <EOL> command = '<STR_LIT>' + input_anat <EOL> result = commands . getoutput ( command ) <EOL> orientation = result [ result . find ( '<STR_LIT>' ) + <NUM_LIT:15> ] + result [ result . find ( '<STR_LIT>' ) + <NUM_LIT:15> ] + result [ result . find ( '<STR_LIT>' ) + <NUM_LIT:15> ] <EOL> if orientation != '<STR_LIT>' : <EOL> print '<STR_LIT>' <EOL> sct . run ( sct . fsloutput + '<STR_LIT>' ) <EOL> sct . run ( sct . fsloutput + '<STR_LIT>' ) <EOL> anat_file = nibabel . load ( '<STR_LIT>' ) <EOL> anat = anat_file . get_data ( ) <EOL> hdr = anat_file . get_header ( ) <EOL> dims = hdr [ '<STR_LIT>' ] <EOL> scales = hdr [ '<STR_LIT>' ] <EOL> centerline_file = nibabel . load ( '<STR_LIT>' ) <EOL> centerline = centerline_file . get_data ( ) <EOL> else : <EOL> anat_file = nibabel . load ( input_anat ) <EOL> anat = anat_file . get_data ( ) <EOL> hdr = anat_file . get_header ( ) <EOL> dims = hdr [ '<STR_LIT>' ] <EOL> scales = hdr [ '<STR_LIT>' ] <EOL> centerline_file = nibabel . load ( input_centerline ) <EOL> centerline = centerline_file . get_data ( ) <EOL> shift_AP = label . shift_AP * scales [ <NUM_LIT:1> ] <EOL> size_AP = label . size_AP * scales [ <NUM_LIT:1> ] <EOL> size_RL = label . size_RL * scales [ <NUM_LIT:3> ] <EOL> np . uint16 ( anat ) <EOL> X , Y , Z = np . where ( centerline > <NUM_LIT:0> ) <EOL> j = np . argsort ( Y ) <EOL> y = Y [ j ] <EOL> x = X [ j ] <EOL> z = Z [ j ] <EOL> index = <NUM_LIT:0> <EOL> for i in range ( len ( y ) - <NUM_LIT:1> ) : <EOL> if y [ i ] == y [ i + <NUM_LIT:1> ] : <EOL> if index == <NUM_LIT:0> : <EOL> index_double = i <EOL> else : <EOL> index_double = np . resize ( index_double , index + <NUM_LIT:1> ) <EOL> index_double [ index ] = i <EOL> index = index + <NUM_LIT:1> <EOL> mask = np . ones ( len ( y ) , dtype = bool ) <EOL> mask [ index_double ] = False <EOL> y = y [ mask ] <EOL> x = x [ mask ] <EOL> z = z [ mask ] <EOL> x1 = np . round ( x - shift_AP / scales [ <NUM_LIT:1> ] ) <EOL> I = np . zeros ( ( len ( y ) , <NUM_LIT:1> ) ) <EOL> for index in range ( len ( y ) ) : <EOL> lim_plus = index + <NUM_LIT:5> <EOL> lim_minus = index - <NUM_LIT:5> <EOL> if lim_minus < <NUM_LIT:0> : <EOL> lim_minus = <NUM_LIT:0> <EOL> if lim_plus >= len ( x1 ) : <EOL> lim_plus = len ( x1 ) - <NUM_LIT:1> <EOL> Vx = x1 [ lim_plus ] - x1 [ lim_minus ] <EOL> Vz = z [ lim_plus ] - z [ lim_minus ] <EOL> Vy = y [ lim_plus ] - y [ lim_minus ] <EOL> d = Vx * x1 [ index ] + Vy * y [ index ] + Vz * z [ index ] <EOL> for i_slice_RL in range ( <NUM_LIT:2> * np . int ( round ( size_RL / scales [ <NUM_LIT:3> ] ) ) ) : <EOL> for i_slice_AP in range ( <NUM_LIT:2> * np . int ( round ( size_AP / scales [ <NUM_LIT:1> ] ) ) ) : <EOL> result = ( d - Vx * ( x1 [ index ] + i_slice_AP - size_AP - <NUM_LIT:1> ) - Vz * z [ index ] ) / Vy <EOL> if result > anat . shape [ <NUM_LIT:1> ] : <EOL> result = anat . shape [ <NUM_LIT:1> ] <EOL> I [ index ] = I [ index ] + anat [ np . int ( round ( x1 [ index ] + i_slice_AP - size_AP - <NUM_LIT:1> ) ) , np . int ( round ( result ) ) , np . int ( round ( z [ index ] + i_slice_RL - size_RL - <NUM_LIT:1> ) ) ] <EOL> start_centerline_y = y [ <NUM_LIT:0> ] <EOL> X = np . where ( I == <NUM_LIT:0> ) <EOL> mask2 = np . ones ( ( len ( y ) , <NUM_LIT:1> ) , dtype = bool ) <EOL> mask2 [ X , <NUM_LIT:0> ] = False <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( I ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . show ( ) <EOL> frequency = scipy . fftpack . fftfreq ( len ( I [ : , <NUM_LIT:0> ] ) , d = <NUM_LIT:1> ) <EOL> z = np . abs ( scipy . fftpack . fft ( I [ : , <NUM_LIT:0> ] , n = None , axis = - <NUM_LIT:1> , overwrite_x = False ) ) <EOL> Wn = np . amax ( frequency ) / <NUM_LIT:10> <EOL> N = <NUM_LIT:5> <EOL> b , a = scipy . signal . iirfilter ( N , Wn , rp = None , rs = None , btype = '<STR_LIT>' , analog = False , ftype = '<STR_LIT>' , output = '<STR_LIT>' ) <EOL> I_fit = scipy . signal . filtfilt ( b , a , I [ : , <NUM_LIT:0> ] , axis = - <NUM_LIT:1> , padtype = '<STR_LIT>' , padlen = None ) <EOL> pl . plot ( I [ : , <NUM_LIT:0> ] ) <EOL> pl . plot ( I_fit ) <EOL> pl . show ( ) <EOL> I_detrend = np . zeros ( ( len ( I [ : , <NUM_LIT:0> ] ) , <NUM_LIT:1> ) ) <EOL> I_detrend [ : , <NUM_LIT:0> ] = I [ : , <NUM_LIT:0> ] - I_fit <EOL> I_detrend = I_detrend / ( np . amax ( I_detrend ) ) <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( I_detrend [ : , <NUM_LIT:0> ] ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . show ( ) <EOL> info_1 = input ( '<STR_LIT>' ) <EOL> if info_1 == <NUM_LIT:0> : <EOL> level_start = input ( '<STR_LIT>' ) <EOL> else : <EOL> level_start = <NUM_LIT:2> <EOL> mean_distance_dict = scipy . io . loadmat ( '<STR_LIT>' ) <EOL> mean_distance = ( mean_distance_dict . values ( ) [ <NUM_LIT:2> ] ) . T <EOL> C1C2_distance = mean_distance [ <NUM_LIT:0> : <NUM_LIT:2> ] <EOL> mean_distance = mean_distance [ level_start - <NUM_LIT:1> : len ( mean_distance ) - <NUM_LIT:1> ] <EOL> space = np . linspace ( - <NUM_LIT:5> / scales [ <NUM_LIT:2> ] , <NUM_LIT:5> / scales [ <NUM_LIT:2> ] , round ( <NUM_LIT:11> / scales [ <NUM_LIT:2> ] ) , endpoint = True ) <EOL> pattern = ( np . sinc ( ( space * scales [ <NUM_LIT:2> ] ) / <NUM_LIT:15> ) ) ** ( <NUM_LIT:20> ) <EOL> xmax_pattern = np . argmax ( pattern ) <EOL> pattern1 = np . concatenate ( ( pattern , np . zeros ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - len ( pattern ) ) ) ) <EOL> corr_all = scipy . signal . correlate ( I_detrend [ : , <NUM_LIT:0> ] , pattern1 ) <EOL> loc_corr = np . arange ( - np . round ( ( len ( corr_all ) / <NUM_LIT:2> ) ) , np . round ( len ( corr_all ) / <NUM_LIT:2> ) + <NUM_LIT:2> ) <EOL> index_fp = <NUM_LIT:0> <EOL> count = <NUM_LIT:0> <EOL> for i in range ( len ( corr_all ) ) : <EOL> if corr_all [ i ] > <NUM_LIT:0.1> : <EOL> if i == <NUM_LIT:0> : <EOL> if corr_all [ i ] < corr_all [ i + <NUM_LIT:1> ] : <EOL> index_fp = i <EOL> count = count + <NUM_LIT:1> <EOL> elif i == ( len ( corr_all ) - <NUM_LIT:1> ) : <EOL> if corr_all [ i ] < corr_all [ i - <NUM_LIT:1> ] : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> else : <EOL> if corr_all [ i ] < corr_all [ i + <NUM_LIT:1> ] : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif corr_all [ i ] < corr_all [ i - <NUM_LIT:1> ] : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if i == <NUM_LIT:0> : <EOL> index_fp = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_fp = np . ones ( len ( corr_all ) , dtype = bool ) <EOL> mask_fp [ index_fp ] = False <EOL> value = corr_all [ mask_fp ] <EOL> loc_corr = loc_corr [ mask_fp ] <EOL> loc_corr = loc_corr - I_detrend . shape [ <NUM_LIT:0> ] <EOL> loc_first_peak = xmax_pattern - loc_corr [ np . amax ( np . where ( value > <NUM_LIT:1> ) ) ] <EOL> Mcorr1 = value [ np . amax ( np . where ( value > <NUM_LIT:1> ) ) ] <EOL> if xmax_pattern < loc_first_peak : <EOL> template_truncated = np . concatenate ( ( np . zeros ( ( loc_first_peak - xmax_pattern ) ) , pattern ) ) <EOL> else : <EOL> template_truncated = pattern [ ( xmax_pattern - loc_first_peak - <NUM_LIT:1> ) : ] <EOL> xend = np . amax ( np . where ( template_truncated > <NUM_LIT> ) ) <EOL> pixend = xend - loc_first_peak <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( template_truncated ) <EOL> pl . plot ( I_detrend ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . show ( ) <EOL> loc_peak_I = np . arange ( len ( I_detrend [ : , <NUM_LIT:0> ] ) ) <EOL> count = <NUM_LIT:0> <EOL> index_p = <NUM_LIT:0> <EOL> for i in range ( len ( I_detrend [ : , <NUM_LIT:0> ] ) ) : <EOL> if I_detrend [ i ] > <NUM_LIT> : <EOL> if i == <NUM_LIT:0> : <EOL> if I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i + <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = i <EOL> count = count + <NUM_LIT:1> <EOL> elif i == ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - <NUM_LIT:1> ) : <EOL> if I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i - <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> else : <EOL> if I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i + <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i - <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if i == <NUM_LIT:0> : <EOL> index_p = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_p = np . ones ( len ( I_detrend [ : , <NUM_LIT:0> ] ) , dtype = bool ) <EOL> mask_p [ index_p ] = False <EOL> value_I = I_detrend [ mask_p ] <EOL> loc_peak_I = loc_peak_I [ mask_p ] <EOL> count = <NUM_LIT:0> <EOL> for i in range ( len ( loc_peak_I ) - <NUM_LIT:1> ) : <EOL> if i == <NUM_LIT:0> : <EOL> if loc_peak_I [ i + <NUM_LIT:1> ] - loc_peak_I [ i ] < round ( <NUM_LIT:10> / scales [ <NUM_LIT:1> ] ) : <EOL> index = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if ( loc_peak_I [ i + <NUM_LIT:1> ] - loc_peak_I [ i ] ) < round ( <NUM_LIT:10> / scales [ <NUM_LIT:1> ] ) : <EOL> index = np . resize ( index , count + <NUM_LIT:1> ) <EOL> index [ len ( index ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif ( loc_peak_I [ i ] - loc_peak_I [ i - <NUM_LIT:1> ] ) < round ( <NUM_LIT:10> / scales [ <NUM_LIT:1> ] ) : <EOL> index = np . resize ( index , count + <NUM_LIT:1> ) <EOL> index [ len ( index ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_I = np . ones ( len ( value_I ) , dtype = bool ) <EOL> mask_I [ index ] = False <EOL> value_I = value_I [ mask_I ] <EOL> loc_peak_I = loc_peak_I [ mask_I ] <EOL> from scipy . interpolate import UnivariateSpline <EOL> fit = UnivariateSpline ( loc_peak_I , value_I ) <EOL> P = fit ( np . arange ( len ( I_detrend ) ) ) <EOL> for i in range ( len ( I_detrend ) ) : <EOL> if P [ i ] > <NUM_LIT:0.1> : <EOL> I_detrend [ i , <NUM_LIT:0> ] = I_detrend [ i , <NUM_LIT:0> ] / P [ i ] <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . xlim ( <NUM_LIT:0> , len ( I_detrend ) - <NUM_LIT:1> ) <EOL> pl . plot ( loc_peak_I , value_I ) <EOL> pl . plot ( I_detrend ) <EOL> pl . plot ( P , color = '<STR_LIT:y>' ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . show ( block = False ) <EOL> mean_distance_new = mean_distance <EOL> mean_ratio = np . zeros ( len ( mean_distance ) ) <EOL> L = np . round ( <NUM_LIT> * max ( mean_distance ) ) - np . round ( <NUM_LIT> * min ( mean_distance ) ) <EOL> corr_peak = np . zeros ( ( L , len ( mean_distance ) ) ) <EOL> for i_peak in range ( len ( mean_distance ) ) : <EOL> scale_min = np . round ( <NUM_LIT> * mean_distance_new [ i_peak ] ) - xmax_pattern - pixend <EOL> if scale_min < <NUM_LIT:0> : <EOL> scale_min = <NUM_LIT:0> <EOL> scale_max = np . round ( <NUM_LIT> * mean_distance_new [ i_peak ] ) - xmax_pattern - pixend <EOL> scale_peak = np . arange ( scale_min , scale_max + <NUM_LIT:1> ) <EOL> for i_scale in range ( len ( scale_peak ) ) : <EOL> template_resize_peak = np . concatenate ( [ template_truncated , np . zeros ( scale_peak [ i_scale ] ) , pattern ] ) <EOL> if len ( I_detrend [ : , <NUM_LIT:0> ] ) > len ( template_resize_peak ) : <EOL> template_resize_peak1 = np . concatenate ( ( template_resize_peak , np . zeros ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - len ( template_resize_peak ) ) ) ) <EOL> corr_template = scipy . signal . correlate ( I_detrend [ : , <NUM_LIT:0> ] , template_resize_peak ) <EOL> if len ( I_detrend [ : , <NUM_LIT:0> ] ) > len ( template_resize_peak ) : <EOL> val = np . dot ( I_detrend [ : , <NUM_LIT:0> ] , template_resize_peak1 . T ) <EOL> else : <EOL> I_detrend_2 = np . concatenate ( ( I_detrend [ : , <NUM_LIT:0> ] , np . zeros ( len ( template_resize_peak ) - len ( I_detrend [ : , <NUM_LIT:0> ] ) ) ) ) <EOL> val = np . dot ( I_detrend_2 , template_resize_peak . T ) <EOL> corr_peak [ i_scale , i_peak ] = val <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . xlim ( <NUM_LIT:0> , len ( I_detrend [ : , <NUM_LIT:0> ] ) ) <EOL> pl . plot ( I_detrend [ : , <NUM_LIT:0> ] ) <EOL> pl . plot ( template_resize_peak ) <EOL> pl . show ( block = False ) <EOL> pl . plot ( corr_peak [ : , i_peak ] , marker = '<STR_LIT:+>' , linestyle = '<STR_LIT:None>' , color = '<STR_LIT:r>' ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . show ( block = False ) <EOL> max_peak = np . amax ( corr_peak [ : , i_peak ] ) <EOL> index_scale_peak = np . where ( corr_peak [ : , i_peak ] == max_peak ) <EOL> good_scale_peak = scale_peak [ index_scale_peak ] [ <NUM_LIT:0> ] <EOL> Mcorr = Mcorr1 <EOL> Mcorr = np . resize ( Mcorr , i_peak + <NUM_LIT:2> ) <EOL> Mcorr [ i_peak + <NUM_LIT:1> ] = np . amax ( corr_peak [ : , <NUM_LIT:0> : ( i_peak + <NUM_LIT:1> ) ] ) <EOL> flag = <NUM_LIT:0> <EOL> if i_peak > <NUM_LIT:0> : <EOL> if ( Mcorr [ i_peak + <NUM_LIT:1> ] - Mcorr [ i_peak ] ) < <NUM_LIT> * np . mean ( Mcorr [ <NUM_LIT:1> : i_peak + <NUM_LIT:2> ] - Mcorr [ <NUM_LIT:0> : i_peak + <NUM_LIT:1> ] ) : <EOL> test = i_peak <EOL> template_resize_peak = np . concatenate ( ( template_truncated , np . zeros ( round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend ) , pattern ) ) <EOL> good_scale_peak = np . round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend <EOL> flag = <NUM_LIT:1> <EOL> if i_peak == <NUM_LIT:0> : <EOL> if ( Mcorr [ i_peak + <NUM_LIT:1> ] - Mcorr [ i_peak ] ) < <NUM_LIT> * Mcorr [ <NUM_LIT:0> ] : <EOL> template_resize_peak = np . concatenate ( ( template_truncated , np . zeros ( round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend ) , pattern ) ) <EOL> good_scale_peak = round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend <EOL> flag = <NUM_LIT:1> <EOL> if flag == <NUM_LIT:0> : <EOL> template_resize_peak = np . concatenate ( ( template_truncated , np . zeros ( good_scale_peak ) , pattern ) ) <EOL> mean_distance_new [ i_peak ] = good_scale_peak + xmax_pattern + pixend <EOL> mean_ratio [ i_peak ] = np . mean ( mean_distance_new [ : , <NUM_LIT:0> : i_peak ] / mean_distance [ : , <NUM_LIT:0> : i_peak ] ) <EOL> template_truncated = template_resize_peak <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( I_detrend [ : , <NUM_LIT:0> ] ) <EOL> pl . plot ( template_truncated ) <EOL> pl . xlim ( <NUM_LIT:0> , ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - <NUM_LIT:1> ) ) <EOL> pl . show ( ) <EOL> minpeakvalue = <NUM_LIT:0.5> <EOL> loc_disk = np . arange ( len ( template_truncated ) ) <EOL> count = <NUM_LIT:0> <EOL> index_disk = <NUM_LIT:0> <EOL> for i in range ( len ( template_truncated ) ) : <EOL> if template_truncated [ i ] >= minpeakvalue : <EOL> if i == <NUM_LIT:0> : <EOL> if template_truncated [ i ] < template_truncated [ i + <NUM_LIT:1> ] : <EOL> index_disk = i <EOL> count = count + <NUM_LIT:1> <EOL> elif i == ( len ( template_truncated ) - <NUM_LIT:1> ) : <EOL> if template_truncated [ i ] < template_truncated [ i - <NUM_LIT:1> ] : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> else : <EOL> if template_truncated [ i ] < template_truncated [ i + <NUM_LIT:1> ] : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif template_truncated [ i ] < template_truncated [ i - <NUM_LIT:1> ] : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if i == <NUM_LIT:0> : <EOL> index_disk = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_disk = np . ones ( len ( template_truncated ) , dtype = bool ) <EOL> mask_disk [ index_disk ] = False <EOL> loc_disk = loc_disk [ mask_disk ] <EOL> X1 = np . where ( loc_disk > I_detrend . shape [ <NUM_LIT:0> ] ) <EOL> mask_disk1 = np . ones ( len ( loc_disk ) , dtype = bool ) <EOL> mask_disk1 [ X1 ] = False <EOL> loc_disk = loc_disk [ mask_disk1 ] <EOL> loc_disk = loc_disk + start_centerline_y - <NUM_LIT:1> <EOL> for i in range ( len ( loc_disk ) ) : <EOL> Index = np . array ( np . where ( y == loc_disk [ i ] ) ) . T <EOL> lim_plus = Index + <NUM_LIT:5> <EOL> lim_minus = Index - <NUM_LIT:5> <EOL> if lim_minus < <NUM_LIT:1> : <EOL> lim_minus = <NUM_LIT:1> <EOL> if lim_plus > len ( x ) : <EOL> lim_plus = len ( x ) <EOL> Vx = x [ lim_plus ] - x [ lim_minus ] <EOL> Vz = z [ lim_plus ] - z [ lim_minus ] <EOL> Vy = y [ lim_plus ] - y [ lim_minus ] <EOL> d = Vx * x1 [ Index ] + Vy * y [ Index ] + Vz * z [ Index ] <EOL> intersection = np . ones ( len ( x ) ) <EOL> for j in range ( len ( x ) ) : <EOL> intersection [ j ] = np . abs ( ( Vx * x [ j ] + Vy * y [ j ] + Vz * z [ j ] - d ) ) <EOL> min_intersection = np . amin ( intersection ) <EOL> index_intersection = np . where ( min_intersection == np . amin ( intersection ) ) <EOL> loc_disk [ i ] = y [ index_intersection ] <EOL> center_disk = centerline <EOL> for i in range ( len ( loc_disk ) - <NUM_LIT:1> ) : <EOL> tmp = center_disk [ : , loc_disk [ i ] : loc_disk [ i + <NUM_LIT:1> ] , : ] <EOL> tmp [ np . where ( tmp == <NUM_LIT:1> ) ] = i + level_start <EOL> center_disk [ : , loc_disk [ i ] : loc_disk [ i + <NUM_LIT:1> ] , : ] = tmp <EOL> center_disk [ np . where ( center_disk == <NUM_LIT:1> ) ] = <NUM_LIT:0> <EOL> if level_start == <NUM_LIT:2> : <EOL> center_disk [ x [ <NUM_LIT:0> ] , round ( loc_disk [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:1> ] ) : loc_disk [ <NUM_LIT:0> ] , z [ <NUM_LIT:0> ] ] = <NUM_LIT:2> <EOL> center_disk [ x [ <NUM_LIT:0> ] , round ( loc_disk [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:1> ] ) : round ( loc_disk [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:1> ] - <NUM_LIT:1> ) , z [ <NUM_LIT:0> ] ] = <NUM_LIT:1> <EOL> if orientation != '<STR_LIT>' : <EOL> a = orientation [ <NUM_LIT:0> ] <EOL> b = orientation [ <NUM_LIT:1> ] <EOL> c = orinetation [ <NUM_LIT:2> ] <EOL> if a == '<STR_LIT:A>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:P>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:S>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:I>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:R>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:L>' : a = '<STR_LIT>' <EOL> if b == '<STR_LIT:A>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:P>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:S>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:I>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:R>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:L>' : b = '<STR_LIT>' <EOL> if c == '<STR_LIT:A>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:P>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:S>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:I>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:R>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:L>' : c = '<STR_LIT>' <EOL> command = fsloutput + '<STR_LIT>' + label . input_path + label . input_anat + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + label . input_path + label . input_anat + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> flag = <NUM_LIT:0> <EOL> if flag == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT>' , output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> else : <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_centerline_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_surface_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> else : <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_centerline_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_surface_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> def labeling_vertebrae_T2 ( label ) : <EOL> if orientation != '<STR_LIT>' : <EOL> input_anat_reorient = label . input_path + label . input_anat + '<STR_LIT>' <EOL> command = '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + input_anat_reorient <EOL> result = commands . getoutput ( command ) <EOL> if label . segmentation_do == <NUM_LIT:1> : <EOL> input_centerline_reorient = label . output_path + label . segmentation_centerline + '<STR_LIT>' <EOL> else : <EOL> input_centerline_reorient = label . input_path + label . input_centerline + '<STR_LIT>' <EOL> command = '<STR_LIT>' + input_centerline + '<STR_LIT:U+0020>' + input_centerline_reorient <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> if label . segmentation_do == <NUM_LIT:1> : <EOL> input_surface_reorinet = label . output_path + label . segmentation_surface + '<STR_LIT>' <EOL> else : <EOL> input_surface_reorinet = label . input_path + label . input_surface + '<STR_LIT>' <EOL> command = '<STR_LIT>' + input_surface + '<STR_LIT:U+0020>' + input_surface_reorient <EOL> command = fsloutput + '<STR_LIT>' + label . input_path + label . input_anat + '<STR_LIT>' + '<STR_LIT>' + label . input_path + label . input_anat + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if label . segmentation_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + label . output_path + label . segmentation_centerline + '<STR_LIT>' + '<STR_LIT>' + label . output_path + label . segmentation_centerline + '<STR_LIT>' <EOL> else : <EOL> command = fsloutput + '<STR_LIT>' + label . input_path + label . input_centerline + '<STR_LIT>' + '<STR_LIT>' + label . input_path + label . input_centerline + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> if label . segmentation_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + label . output_path + label . segmentation_surface + '<STR_LIT>' + '<STR_LIT>' + label . output_path + label . segmentation_surface + '<STR_LIT>' <EOL> else : <EOL> command = fsloutput + '<STR_LIT>' + label . input_path + label . input_surface + '<STR_LIT>' + '<STR_LIT>' + label . input_path + label . input_surface + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> anat_file = nibabel . load ( input_anat_reorient ) <EOL> anat = anat_file . get_data ( ) <EOL> hdr = anat_file . get_header ( ) <EOL> dims = hdr [ '<STR_LIT>' ] <EOL> scales = hdr [ '<STR_LIT>' ] <EOL> centerline_file = nibabel . load ( input_centerline_reorient ) <EOL> centerline = centerline_file . get_data ( ) <EOL> else : <EOL> anat_file = nibabel . load ( input_anat ) <EOL> anat = anat_file . get_data ( ) <EOL> hdr = anat_file . get_header ( ) <EOL> dims = hdr [ '<STR_LIT>' ] <EOL> scales = hdr [ '<STR_LIT>' ] <EOL> centerline_file = nibabel . load ( input_centerline ) <EOL> centerline = centerline_file . get_data ( ) <EOL> shift_AP = label . shift_AP * scales [ <NUM_LIT:1> ] <EOL> size_AP = label . size_AP * scales [ <NUM_LIT:1> ] <EOL> size_RL = label . size_RL * scales [ <NUM_LIT:3> ] <EOL> np . uint16 ( anat ) <EOL> X , Y , Z = np . where ( centerline > <NUM_LIT:0> ) <EOL> j = np . argsort ( Y ) <EOL> y = Y [ j ] <EOL> x = X [ j ] <EOL> z = Z [ j ] <EOL> index = <NUM_LIT:0> <EOL> for i in range ( len ( y ) - <NUM_LIT:1> ) : <EOL> if y [ i ] == y [ i + <NUM_LIT:1> ] : <EOL> if index == <NUM_LIT:0> : <EOL> index_double = i <EOL> else : <EOL> index_double = np . resize ( index_double , index + <NUM_LIT:1> ) <EOL> index_double [ index ] = i <EOL> index = index + <NUM_LIT:1> <EOL> mask = np . ones ( len ( y ) , dtype = bool ) <EOL> mask [ index_double ] = False <EOL> y = y [ mask ] <EOL> x = x [ mask ] <EOL> z = z [ mask ] <EOL> x1 = np . round ( x - shift_AP / scales [ <NUM_LIT:1> ] ) <EOL> I = np . zeros ( ( len ( y ) , <NUM_LIT:1> ) ) <EOL> for index in range ( len ( y ) ) : <EOL> lim_plus = index + <NUM_LIT:5> <EOL> lim_minus = index - <NUM_LIT:5> <EOL> if lim_minus < <NUM_LIT:0> : <EOL> lim_minus = <NUM_LIT:0> <EOL> if lim_plus >= len ( x1 ) : <EOL> lim_plus = len ( x1 ) - <NUM_LIT:1> <EOL> Vx = x1 [ lim_plus ] - x1 [ lim_minus ] <EOL> Vz = z [ lim_plus ] - z [ lim_minus ] <EOL> Vy = y [ lim_plus ] - y [ lim_minus ] <EOL> d = Vx * x1 [ index ] + Vy * y [ index ] + Vz * z [ index ] <EOL> for i_slice_RL in range ( <NUM_LIT:2> * np . int ( round ( size_RL / scales [ <NUM_LIT:3> ] ) ) ) : <EOL> for i_slice_AP in range ( <NUM_LIT:2> * np . int ( round ( size_AP / scales [ <NUM_LIT:1> ] ) ) ) : <EOL> result = ( d - Vx * ( x1 [ index ] + i_slice_AP - size_AP - <NUM_LIT:1> ) - Vz * z [ index ] ) / Vy <EOL> if result > anat . shape [ <NUM_LIT:1> ] : <EOL> result = anat . shape [ <NUM_LIT:1> ] <EOL> I [ index ] = I [ index ] + anat [ np . int ( round ( x1 [ index ] + i_slice_AP - size_AP - <NUM_LIT:1> ) ) , np . int ( round ( result ) ) , np . int ( round ( z [ index ] + i_slice_RL - size_RL - <NUM_LIT:1> ) ) ] <EOL> start_centerline_y = y [ <NUM_LIT:0> ] <EOL> X = np . where ( I == <NUM_LIT:0> ) <EOL> mask2 = np . ones ( ( len ( y ) , <NUM_LIT:1> ) , dtype = bool ) <EOL> mask2 [ X , <NUM_LIT:0> ] = False <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( I ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . show ( block = False ) <EOL> I_detrend = scipy . signal . detrend ( I , axis = <NUM_LIT:0> ) <EOL> I_detrend = I_detrend / abs ( ( np . amin ( I_detrend ) ) ) <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( I_detrend ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . show ( block = False ) <EOL> info_1 = input ( '<STR_LIT>' ) <EOL> if info_1 == <NUM_LIT:0> : <EOL> level_start = input ( '<STR_LIT>' ) <EOL> else : <EOL> level_start = <NUM_LIT:2> <EOL> mean_distance_dict = scipy . io . loadmat ( '<STR_LIT>' ) <EOL> mean_distance = ( mean_distance_dict . values ( ) [ <NUM_LIT:2> ] ) . T <EOL> C1C2_distance = mean_distance [ <NUM_LIT:0> : <NUM_LIT:2> ] <EOL> mean_distance = mean_distance [ level_start + <NUM_LIT:1> : len ( mean_distance ) - <NUM_LIT:1> ] <EOL> space = np . linspace ( - <NUM_LIT:5> / scales [ <NUM_LIT:2> ] , <NUM_LIT:5> / scales [ <NUM_LIT:2> ] , round ( <NUM_LIT:11> / scales [ <NUM_LIT:2> ] ) , endpoint = True ) <EOL> pattern = ( np . sinc ( ( space * scales [ <NUM_LIT:2> ] ) / <NUM_LIT:15> ) ) ** ( <NUM_LIT:20> ) <EOL> xmax_pattern = np . argmin ( pattern ) <EOL> pixend = len ( pattern ) - xmax_pattern <EOL> pattern1 = np . concatenate ( ( pattern , np . zeros ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - len ( pattern ) ) ) ) <EOL> corr_all = scipy . signal . correlate ( I_detrend [ : , <NUM_LIT:0> ] , pattern1 ) <EOL> loc_corr = np . arange ( - np . round ( ( len ( corr_all ) / <NUM_LIT:2> ) ) , np . round ( len ( corr_all ) / <NUM_LIT:2> ) + <NUM_LIT:2> ) <EOL> index_fp = <NUM_LIT:0> <EOL> count = <NUM_LIT:0> <EOL> for i in range ( len ( corr_all ) ) : <EOL> if corr_all [ i ] > <NUM_LIT:0.1> : <EOL> if i == <NUM_LIT:0> : <EOL> if corr_all [ i ] < corr_all [ i + <NUM_LIT:1> ] : <EOL> index_fp = i <EOL> count = count + <NUM_LIT:1> <EOL> elif i == ( len ( corr_all ) - <NUM_LIT:1> ) : <EOL> if corr_all [ i ] < corr_all [ i - <NUM_LIT:1> ] : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> else : <EOL> if corr_all [ i ] < corr_all [ i + <NUM_LIT:1> ] : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif corr_all [ i ] < corr_all [ i - <NUM_LIT:1> ] : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if i == <NUM_LIT:0> : <EOL> index_fp = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> index_fp = np . resize ( index_fp , count + <NUM_LIT:1> ) <EOL> index_fp [ len ( index_fp ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_fp = np . ones ( len ( corr_all ) , dtype = bool ) <EOL> mask_fp [ index_fp ] = False <EOL> value = corr_all [ mask_fp ] <EOL> loc_corr = loc_corr [ mask_fp ] <EOL> loc_corr = loc_corr - I_detrend . shape [ <NUM_LIT:0> ] <EOL> loc_first_peak = xmax_pattern - loc_corr [ np . amax ( np . where ( value > <NUM_LIT> ) ) ] <EOL> Mcorr1 = value [ np . amax ( np . where ( value > <NUM_LIT> ) ) ] <EOL> if loc_first_peak >= <NUM_LIT:0> : <EOL> template_truncated = pattern [ ( loc_first_peak + <NUM_LIT:1> ) : ] <EOL> else : <EOL> template_truncated = np . concatenate ( ( np . zeros ( abs ( loc_first_peak ) ) , pattern ) ) <EOL> xend = len ( template_truncated ) <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( template_truncated ) <EOL> pl . plot ( I_detrend ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . show ( block = False ) <EOL> I_detrend [ : , <NUM_LIT:0> ] = scipy . ndimage . filters . gaussian_filter1d ( I_detrend [ : , <NUM_LIT:0> ] , <NUM_LIT:10> ) <EOL> loc_peak_I = np . arange ( len ( I_detrend [ : , <NUM_LIT:0> ] ) ) <EOL> count = <NUM_LIT:0> <EOL> index_p = <NUM_LIT:0> <EOL> for i in range ( len ( I_detrend [ : , <NUM_LIT:0> ] ) ) : <EOL> if I_detrend [ i ] > <NUM_LIT> : <EOL> if i == <NUM_LIT:0> : <EOL> if I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i + <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = i <EOL> count = count + <NUM_LIT:1> <EOL> elif i == ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - <NUM_LIT:1> ) : <EOL> if I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i - <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> else : <EOL> if I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i + <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif I_detrend [ i , <NUM_LIT:0> ] < I_detrend [ i - <NUM_LIT:1> , <NUM_LIT:0> ] : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if i == <NUM_LIT:0> : <EOL> index_p = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> index_p = np . resize ( index_p , count + <NUM_LIT:1> ) <EOL> index_p [ len ( index_p ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_p = np . ones ( len ( I_detrend [ : , <NUM_LIT:0> ] ) , dtype = bool ) <EOL> mask_p [ index_p ] = False <EOL> value_I = I_detrend [ mask_p ] <EOL> loc_peak_I = loc_peak_I [ mask_p ] <EOL> count = <NUM_LIT:0> <EOL> for i in range ( len ( loc_peak_I ) - <NUM_LIT:1> ) : <EOL> if i == <NUM_LIT:0> : <EOL> if loc_peak_I [ i + <NUM_LIT:1> ] - loc_peak_I [ i ] < round ( <NUM_LIT:10> / scales [ <NUM_LIT:1> ] ) : <EOL> index = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if ( loc_peak_I [ i + <NUM_LIT:1> ] - loc_peak_I [ i ] ) < round ( <NUM_LIT:10> / scales [ <NUM_LIT:1> ] ) : <EOL> index = np . resize ( index , count + <NUM_LIT:1> ) <EOL> index [ len ( index ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif ( loc_peak_I [ i ] - loc_peak_I [ i - <NUM_LIT:1> ] ) < round ( <NUM_LIT:10> / scales [ <NUM_LIT:1> ] ) : <EOL> index = np . resize ( index , count + <NUM_LIT:1> ) <EOL> index [ len ( index ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_I = np . ones ( len ( value_I ) , dtype = bool ) <EOL> mask_I [ index ] = False <EOL> value_I = - value_I [ mask_I ] <EOL> loc_peak_I = loc_peak_I [ mask_I ] <EOL> from scipy . interpolate import UnivariateSpline <EOL> fit = UnivariateSpline ( loc_peak_I , value_I ) <EOL> P = fit ( np . arange ( len ( I_detrend ) ) ) <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . xlim ( <NUM_LIT:0> , len ( I_detrend ) - <NUM_LIT:1> ) <EOL> pl . plot ( loc_peak_I , value_I ) <EOL> pl . plot ( I_detrend ) <EOL> pl . plot ( P ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . xlabel ( '<STR_LIT>' ) <EOL> pl . ylabel ( '<STR_LIT>' ) <EOL> pl . show ( block = False ) <EOL> for i in range ( len ( I_detrend ) ) : <EOL> if P [ i ] > <NUM_LIT:0.1> : <EOL> I_detrend [ i , <NUM_LIT:0> ] = I_detrend [ i , <NUM_LIT:0> ] / abs ( P [ i ] ) <EOL> mean_distance_new = mean_distance <EOL> mean_ratio = np . zeros ( len ( mean_distance ) ) <EOL> L = np . round ( <NUM_LIT> * max ( mean_distance ) ) - np . round ( <NUM_LIT> * min ( mean_distance ) ) <EOL> corr_peak = np . nan ( np . zeros ( ( L , len ( mean_distance ) ) ) ) <EOL> for i_peak in range ( len ( mean_distance ) ) : <EOL> scale_min = np . round ( <NUM_LIT> * mean_distance_new [ i_peak ] ) - xmax_pattern - pixend <EOL> if scale_min < <NUM_LIT:0> : <EOL> scale_min = <NUM_LIT:0> <EOL> scale_max = np . round ( <NUM_LIT> * mean_distance_new [ i_peak ] ) - xmax_pattern - pixend <EOL> scale_peak = np . arange ( scale_min , scale_max + <NUM_LIT:1> ) <EOL> for i_scale in range ( len ( scale_peak ) ) : <EOL> template_resize_peak = np . concatenate ( [ template_truncated , np . zeros ( scale_peak [ i_scale ] ) , pattern ] ) <EOL> if len ( I_detrend [ : , <NUM_LIT:0> ] ) > len ( template_resize_peak ) : <EOL> template_resize_peak1 = np . concatenate ( ( template_resize_peak , np . zeros ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - len ( template_resize_peak ) ) ) ) <EOL> corr_template = scipy . signal . correlate ( I_detrend [ : , <NUM_LIT:0> ] , template_resize_peak ) <EOL> if len ( I_detrend [ : , <NUM_LIT:0> ] ) > len ( template_resize_peak ) : <EOL> val = np . dot ( I_detrend [ : , <NUM_LIT:0> ] , template_resize_peak1 . T ) <EOL> else : <EOL> I_detrend_2 = np . concatenate ( ( I_detrend [ : , <NUM_LIT:0> ] , np . zeros ( len ( template_resize_peak ) - len ( I_detrend [ : , <NUM_LIT:0> ] ) ) ) ) <EOL> val = np . dot ( I_detrend_2 , template_resize_peak . T ) <EOL> corr_peak [ i_scale , i_peak ] = val <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . xlim ( <NUM_LIT:0> , len ( I_detrend [ : , <NUM_LIT:0> ] ) ) <EOL> pl . plot ( I_detrend [ : , <NUM_LIT:0> ] ) <EOL> pl . plot ( template_resize_peak ) <EOL> pl . show ( block = False ) <EOL> pl . plot ( corr_peak [ : , i_peak ] , marker = '<STR_LIT:+>' , linestyle = '<STR_LIT:None>' , color = '<STR_LIT:r>' ) <EOL> pl . title ( '<STR_LIT>' ) <EOL> pl . show ( block = False ) <EOL> max_peak = np . amax ( corr_peak [ : , i_peak ] ) <EOL> index_scale_peak = np . where ( corr_peak [ : , i_peak ] == max_peak ) <EOL> good_scale_peak = scale_peak [ index_scale_peak ] [ <NUM_LIT:0> ] <EOL> Mcorr = Mcorr1 <EOL> Mcorr = np . resize ( Mcorr , i_peak + <NUM_LIT:2> ) <EOL> Mcorr [ i_peak + <NUM_LIT:1> ] = np . amax ( corr_peak [ : , <NUM_LIT:0> : ( i_peak + <NUM_LIT:1> ) ] ) <EOL> flag = <NUM_LIT:0> <EOL> if i_peak > <NUM_LIT:0> : <EOL> if ( Mcorr [ i_peak + <NUM_LIT:1> ] - Mcorr [ i_peak ] ) < <NUM_LIT> * np . mean ( Mcorr [ <NUM_LIT:1> : i_peak + <NUM_LIT:2> ] - Mcorr [ <NUM_LIT:0> : i_peak + <NUM_LIT:1> ] ) : <EOL> test = i_peak <EOL> template_resize_peak = np . concatenate ( ( template_truncated , np . zeros ( round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend ) , pattern ) ) <EOL> good_scale_peak = np . round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend <EOL> flag = <NUM_LIT:1> <EOL> if i_peak == <NUM_LIT:0> : <EOL> if ( Mcorr [ i_peak + <NUM_LIT:1> ] - Mcorr [ i_peak ] ) < <NUM_LIT> * Mcorr [ <NUM_LIT:0> ] : <EOL> template_resize_peak = np . concatenate ( ( template_truncated , np . zeros ( round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend ) , pattern ) ) <EOL> good_scale_peak = round ( mean_distance [ i_peak ] ) - xmax_pattern - pixend <EOL> flag = <NUM_LIT:1> <EOL> if flag == <NUM_LIT:0> : <EOL> template_resize_peak = np . concatenate ( ( template_truncated , np . zeros ( good_scale_peak ) , pattern ) ) <EOL> mean_distance_new [ i_peak ] = good_scale_peak + xmax_pattern + pixend <EOL> mean_ratio [ i_peak ] = np . mean ( mean_distance_new [ : , <NUM_LIT:0> : i_peak ] / mean_distance [ : , <NUM_LIT:0> : i_peak ] ) <EOL> template_truncated = template_resize_peak <EOL> if label . verbose == <NUM_LIT:1> : <EOL> pl . plot ( I_detrend [ : , <NUM_LIT:0> ] ) <EOL> pl . plot ( template_truncated ) <EOL> pl . xlim ( <NUM_LIT:0> , ( len ( I_detrend [ : , <NUM_LIT:0> ] ) - <NUM_LIT:1> ) ) <EOL> pl . show ( block = False ) <EOL> minpeakvalue = <NUM_LIT:0.5> <EOL> loc_disk = np . arange ( len ( template_truncated ) ) <EOL> count = <NUM_LIT:0> <EOL> index_disk = <NUM_LIT:0> <EOL> for i in range ( len ( template_truncated ) ) : <EOL> if template_truncated [ i ] >= minpeakvalue : <EOL> if i == <NUM_LIT:0> : <EOL> if template_truncated [ i ] < template_truncated [ i + <NUM_LIT:1> ] : <EOL> index_disk = i <EOL> count = count + <NUM_LIT:1> <EOL> elif i == ( len ( template_truncated ) - <NUM_LIT:1> ) : <EOL> if template_truncated [ i ] < template_truncated [ i - <NUM_LIT:1> ] : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> else : <EOL> if template_truncated [ i ] < template_truncated [ i + <NUM_LIT:1> ] : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> elif template_truncated [ i ] < template_truncated [ i - <NUM_LIT:1> ] : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> if i == <NUM_LIT:0> : <EOL> index_disk = i <EOL> count = count + <NUM_LIT:1> <EOL> else : <EOL> index_disk = np . resize ( index_disk , count + <NUM_LIT:1> ) <EOL> index_disk [ len ( index_disk ) - <NUM_LIT:1> ] = i <EOL> count = count + <NUM_LIT:1> <EOL> mask_disk = np . ones ( len ( template_truncated ) , dtype = bool ) <EOL> mask_disk [ index_disk ] = False <EOL> loc_disk = loc_disk [ mask_disk ] <EOL> X1 = np . where ( loc_disk > I_detrend . shape [ <NUM_LIT:0> ] ) <EOL> mask_disk1 = np . ones ( len ( loc_disk ) , dtype = bool ) <EOL> mask_disk1 [ X1 ] = False <EOL> loc_disk = loc_disk [ mask_disk1 ] <EOL> loc_disk = loc_disk + start_centerline_y - <NUM_LIT:1> <EOL> for i in range ( len ( loc_disk ) ) : <EOL> Index = np . array ( np . where ( y == loc_disk [ i ] ) ) . T <EOL> lim_plus = Index + <NUM_LIT:5> <EOL> lim_minus = Index - <NUM_LIT:5> <EOL> if lim_minus < <NUM_LIT:1> : <EOL> lim_minus = <NUM_LIT:1> <EOL> if lim_plus > len ( x ) : <EOL> lim_plus = len ( x ) <EOL> Vx = x [ lim_plus ] - x [ lim_minus ] <EOL> Vz = z [ lim_plus ] - z [ lim_minus ] <EOL> Vy = y [ lim_plus ] - y [ lim_minus ] <EOL> d = Vx * x1 [ Index ] + Vy * y [ Index ] + Vz * z [ Index ] <EOL> intersection = np . ones ( len ( x ) ) <EOL> for j in range ( len ( x ) ) : <EOL> intersection [ j ] = np . abs ( ( Vx * x [ j ] + Vy * y [ j ] + Vz * z [ j ] - d ) ) <EOL> min_intersection = np . amin ( intersection ) <EOL> index_intersection = np . where ( min_intersection == np . amin ( intersection ) ) <EOL> loc_disk [ i ] = y [ index_intersection ] <EOL> center_disk = np . array ( centerline ) <EOL> for i in range ( len ( loc_disk ) - <NUM_LIT:1> ) : <EOL> tmp = center_disk [ loc_disk [ i ] : loc_disk [ i + <NUM_LIT:1> ] ] <EOL> tmp [ np . where ( tmp == <NUM_LIT:1> ) ] = i + level_start <EOL> center_disk [ loc_disk [ i ] : loc_disk [ i + <NUM_LIT:1> ] ] = tmp <EOL> center_disk [ np . where ( center_disk == <NUM_LIT:1> ) ] = <NUM_LIT:0> <EOL> if level_start == <NUM_LIT:2> : <EOL> center_disk [ x [ <NUM_LIT:0> ] , round ( loc_disk [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:1> ] ) : loc_disk [ <NUM_LIT:0> ] , z [ <NUM_LIT:0> ] ] = <NUM_LIT:2> <EOL> center_disk [ x [ <NUM_LIT:0> ] , round ( loc_disk [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:1> ] ) : round ( loc_disk [ <NUM_LIT:0> ] - C1C2_distance [ <NUM_LIT:1> ] - <NUM_LIT:1> ) , z [ <NUM_LIT:0> ] ] = <NUM_LIT:1> <EOL> if orientation != '<STR_LIT>' : <EOL> a = orientation [ <NUM_LIT:0> ] <EOL> b = orientation [ <NUM_LIT:1> ] <EOL> c = orinetation [ <NUM_LIT:2> ] <EOL> if a == '<STR_LIT:A>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:P>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:S>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:I>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:R>' : a = '<STR_LIT>' <EOL> if a == '<STR_LIT:L>' : a = '<STR_LIT>' <EOL> if b == '<STR_LIT:A>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:P>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:S>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:I>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:R>' : b = '<STR_LIT>' <EOL> if b == '<STR_LIT:L>' : b = '<STR_LIT>' <EOL> if c == '<STR_LIT:A>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:P>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:S>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:I>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:R>' : c = '<STR_LIT>' <EOL> if c == '<STR_LIT:L>' : c = '<STR_LIT>' <EOL> command = fsloutput + '<STR_LIT>' + label . input_path + label . input_anat + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + label . input_path + label . input_anat + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> flag = <NUM_LIT:0> <EOL> if flag == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT>' , output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> else : <EOL> command = fsloutput + '<STR_LIT>' + output_centerline_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_centerline_vertebra <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + output_surface_vertebra + '<STR_LIT:U+0020>' + a + '<STR_LIT:U+0020>' + b + '<STR_LIT:U+0020>' + c + '<STR_LIT:U+0020>' + output_surface_vertebra <EOL> result = commands . getoutput ( command ) <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_centerline_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_surface_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> else : <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_centerline_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> if surface_do == <NUM_LIT:1> : <EOL> command = fsloutput + '<STR_LIT>' + input_anat + '<STR_LIT:U+0020>' + output_surface_vertebra + '<STR_LIT>' <EOL> result = commands . getoutput ( command ) <EOL> def j_detrend_new_v2 ( data , deg , func_type , robust ) : <EOL> mask = <NUM_LIT:0> <EOL> size = data . shape <EOL> if len ( size ) == <NUM_LIT:2> : <EOL> size = np . resize ( size , len ( size ) + <NUM_LIT:2> ) <EOL> size [ <NUM_LIT:2> ] = <NUM_LIT:1> <EOL> size [ <NUM_LIT:3> ] = <NUM_LIT:1> <EOL> if size [ <NUM_LIT:3> ] == <NUM_LIT:1> : <EOL> data2d = data . T <EOL> else : <EOL> data2d = ( np . reshape ( data , ( size [ <NUM_LIT:0> ] * size [ <NUM_LIT:1> ] * size [ <NUM_LIT:2> ] , size [ <NUM_LIT:3> ] ) ) ) . T <EOL> if mask != <NUM_LIT:1> : <EOL> mask1d = np . reshape ( mask , ( ( size [ <NUM_LIT:0> ] * size [ <NUM_LIT:1> ] * size [ <NUM_LIT:2> ] ) , <NUM_LIT:1> ) ) <EOL> else : <EOL> mask1d = np . ones ( ( ( size [ <NUM_LIT:0> ] * size [ <NUM_LIT:1> ] * size [ <NUM_LIT:2> ] ) , <NUM_LIT:1> ) ) <EOL> index_mask = np . where ( mask1d != <NUM_LIT:0> ) <EOL> nb_samples = data2d . shape [ <NUM_LIT:0> ] <EOL> nb_vectors = len ( index_mask ) <EOL> if func_type == '<STR_LIT>' : <EOL> D = ( np . arange ( - <NUM_LIT:1> , ( <NUM_LIT:1> + <NUM_LIT:2> / ( nb_samples - <NUM_LIT:1> ) ) , <NUM_LIT:2> / ( nb_samples - <NUM_LIT:1> ) ) ) . T <EOL> else : <EOL> N = nb_samples <EOL> K = deg <EOL> d = <NUM_LIT:0> <EOL> n = np . arange ( <NUM_LIT:0> , N ) <EOL> C = np . zeros ( ( N , K ) ) <EOL> C [ : , <NUM_LIT:0> ] = np . ones ( N ) / math . sqrt ( N ) <EOL> for k in range ( <NUM_LIT:1> , K ) : <EOL> C [ : , k ] = math . sqrt ( <NUM_LIT:2> / N ) * ( np . cos ( np . pi * ( <NUM_LIT:2> * n + <NUM_LIT:1> ) * ( k - <NUM_LIT:1> ) / ( <NUM_LIT:2> * N ) ) ) <EOL> D = C * ( math . sqrt ( nb_samples ) ) <EOL> data2d_d = np . zeros ( ( size [ <NUM_LIT:0> ] * size [ <NUM_LIT:1> ] * size [ <NUM_LIT:2> ] , size [ <NUM_LIT:3> ] ) ) <EOL> for i_vect in range ( <NUM_LIT:0> , nb_vectors ) : <EOL> data1d = data2d [ index_mask [ i_vect ] , : ] <EOL> if robust == <NUM_LIT:1> : <EOL> l = <NUM_LIT:1> <EOL> else : <EOL> l = np . dot ( np . dot ( np . linalg . pinv ( np . dot ( D . T , D ) ) , D . T ) , data1d ) <EOL> Dl = np . dot ( D , l ) <EOL> res_l = data1d - Dl <EOL> data1d_d = data1d - Dl <EOL> data2d_d [ index_mask [ i_vect ] , : ] = data1d_d <EOL> data_d = data2d_d <EOL> return data_d <EOL> def usage ( ) : <EOL> print '<STR_LIT:\n>' '<STR_LIT>' '<STR_LIT>' '<STR_LIT>' '<STR_LIT>' '<STR_LIT:\n>' '<STR_LIT>' '<STR_LIT>' '<STR_LIT>' '<STR_LIT:\n>' '<STR_LIT>' '<STR_LIT>' '<STR_LIT>' '<STR_LIT:\n>' '<STR_LIT>' '<STR_LIT>' <EOL> sys . exit ( <NUM_LIT:2> ) <EOL> def func ( x , a , b , c , d , e , f ) : <EOL> return a * x ** <NUM_LIT:5> + b * x ** <NUM_LIT:4> + c * x ** <NUM_LIT:3> + d * x ** <NUM_LIT:2> + e * x + f <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> main ( ) </s>
