3 dZ@sLddlZddljZddlZddlZddlZ ddl m Z ddl m Z ddlmZmZmZmZmZmZmZmZe jdZd+d d Zd d Zd dZddZddZd,ddZddZGdddeZ GdddeZ!GdddeZ"GdddeZ#Gdd d eZ$Gd!d"d"eZ%d#d$Z&Gd%d&d&eZ'Gd'd(d(eZ(Gd)d*d*eZ)dS)-N)logging)split_filename) BaseInterfaceBaseInterfaceInputSpectraitsFile TraitedSpec DirectoryOutputMultiPath isdefinedznipype.interfaceFcCsdtj|}|jddkr,|r(tjdgSdStjtj|dddjdd}|rZtj|Stj|S)a Euclidean length of track line Parameters ---------- xyz : array-like shape (N,3) array representing x,y,z of N points in a track along : bool, optional If True, return array giving cumulative length along track, otherwise (default) return scalar giving total length. Returns ------- L : scalar or array shape (N-1,) scalar in case of `along` == False, giving total length, array if `along` == True, giving cumulative lengths. Examples -------- >>> xyz = np.array([[1,1,1],[2,3,4],[0,0,0]]) >>> expected_lens = np.sqrt([1+2**2+3**2, 2**2+3**2+4**2]) >>> length(xyz) == expected_lens.sum() True >>> len_along = length(xyz, along=True) >>> np.allclose(len_along, expected_lens.cumsum()) True >>> length([]) 0 >>> length([[1, 2, 3]]) 0 >>> length([], along=True) array([0]) rr)axis)npZasarrayshapearraysqrtZdiffsumZcumsum)xyzZalongdistsr=/tmp/pip-build-7vycvbft/nipype/nipype/interfaces/cmtk/cmtk.pylengths"    rc Cst|}g}xtd|D]~}t||dft|d}t||dft|d}t||dft|d}||||fdks|j||||fqWttj|j}|S)Nrrr) lenrangeintfloatappendlistdictfromkeyskeys) ZpointsmmroiData voxelSizeZn_points rois_crossedjxyzrrrget_rois_crossedIsr*cCstj||ftjd}xd|D]\}xVt|D]J\}}x@t|D]4\}}||kr:||ks:||d|dfd7<q:Wq(WqW||j}|S)N)dtyper)rzerosuint enumerateT)n_roislist_of_roi_crossed_listsconnectivity_matrixr%Zidx_iZroi_iZidx_jZroi_jrrrget_connectivity_matrixYs ( r3cCst|}d }g}g}xt|D]x\}} tttd||} | |krb| ddkrb| }td|t| d||} t| dkr|jt| |j|qWt ||} |t|} t j d| | d||t j d|| d| d|t j d| |fS) z-Create the intersection arrays for each fiberrdrz%4.0f%%zpFound %i (%f percent out of %i fibers) fibers that start or terminate in a voxel which is not labeled. (orphans)gY@zValid fibers: %i (%f percent)z4Returning the intersecting point connectivity matrix) rr.rroundrprintr*rrr3ifloggerinfo)Z streamlinesr#r$r0Zn_fibZpcfinal_fiber_idsr1iZfiberZpcNr%r2disrrrcreate_allpoints_cmates0      r=cCst|}tj|ddf}tj|ddf}xt|D]x\}}|d}|dddf||dddf<|dddf||dddf<|dddf||dddf<|dddf||dddf<t||ddft|d||ddf<t||ddft|d||ddf<t||ddft|d||ddf<t||ddft|d||ddf<t||ddft|d||ddf<t||ddft|d||ddf<q4Wtjd||fS) aCreate the endpoints arrays for each fiber. Parameters ---------- fib : array-like the fibers data voxelSize : tuple 3-tuple containing the voxel size of the ROI image Returns ------- endpoints : ndarray of size [#fibers, 2, 3] containing for each fiber the index of its first and last point in the voxelSize volume endpointsmm : ndarray of size [#fibers, 2, 3] endpoints in millimeter coordinates rrrNrzReturning the endpoint matrixr5r5)rrr,r.rrr8r9)fibr$n endpoints endpointsmmr;fifrrrcreate_endpoints_arrays"(((((, rDcFCsddlj}i}tjdtj|d} tj|d} tjd|tjj|d\} } t | |d<tj |} t j | j }| jj}t| |\}}tjd | t j| |tjd | t j| |t | }tjd |t j|d f}g}g}t|\}}}|d kr tj|}n |dkr tj|}n td|t |j}d|jt|jdkr\|j}ntj}x|jddD]p\}}|jt|f|tt jt j t j!|t|dkdd}t|d|d |d g|jt|d<qrW|rvtjdt"| |||\}} tj|d}!t#| | |!| |d<t j$|}|j}"tj%t j$|}#tj&|#dd}#|"j'dd|#j(ddDd}$xt)|j*dD]}%ydt|||%ddf||%ddf||%dd ff}&t|||%ddf||%ddf||%dd ff}'Wn$t+k rtj,d|%PYnX|&dks0|'dkrH|$d7}$dJ||%df<q|&|ks\|'|krtj,dtj,dtj,d |&|'tj,d!q|'|&kr|&}(|'}&|(}'|&||%df<|'||%df<|j-|&|'g|j-|%|j.|&|'rd"|j/|&|'kr|j/|&|'d"j-|%n|j0|&|'|%gd#qWg})|r6| }*n|}*x$|*D]}+|)j-t1| |+dq@Wt j2|)},t j2|td$}-tjd%|$|$d&||tjd'||$d(|$d&|tj}.|.j3||.j}/|.j}0|.j}1xZ|j(ddD]H\}}2}|j4||2i}3d"|krt |d"|3d)<t j |-dddft|k|-dddft|2k@d}+t5t j|,|+|3d*<t5t j6|,|+|3d+<t5t j7|,|+|3d,<n d|3d)<d|3d*<d|3d+<d|3d,<||2ks|j0||2f|3d"|kr|.j0||2|3d)d-|/j0||2|3d*d-|0j0||2|3d+d-|1j0||2|3d,d-qWtjd.|tj8|tj|tj9|.td$}4d)|4i}5tj9|/t j:d$}6d/|6i}7tj9|0t j:d$}8d0|8i}9tj9|1t j:d$}:d,|:i};|rt|\}}}tj|d1|}tjd6|>|j;|>|=tj|d7|}?tjd8|?|j;|?|7tj|d9|}@tjd:|@|j;|@|9tj|d;|}Atjd<|A|j;|A|;tj|d=}Btjd>|Bt j|B|,tj|d?}Ctjd@|Ct j|Ct j2|t j szcmat..css"|]\}}}|||dfVqdS)weightNr).0uvdrrr szcmat..zNAN INDEXERROR EXCEPTION OCCURED FOR FIBER %s. PLEASE CHECK ENDPOINT GENERATIONzGStart or endpoint of fiber terminate in a voxel which is labeled higherz6than is expected by the parcellation node information.zStart ROI: %i, End ROI: %izThis needs bugfixing!fiblist)rO)r+zpFound %i (%f percent out of %i fibers) fibers that start or terminate in a voxel which is not labeled. (orphans)gY@zValid fibers: %i (%f%%)r4Znumber_of_fibersZfiber_length_meanZfiber_length_medianZfiber_length_std)rIzWriting network as %sZmean_fiber_lengthZmedian_fiber_lengthZ_intersectionsz"Writing intersection network as %sz.matzWriting matlab matrix as %s intersectionsz!Writing intersection matrix as %sZ_mean_fiber_lengthz-Writing matlab mean fiber length matrix as %sZ_median_fiber_lengthz/Writing matlab median fiber length matrix as %sZ_fiber_length_stdz2Writing matlab fiber length deviation matrix as %sz_final_fiberslength.npyz&Storing final fiber length array as %sz_filtered_fiberslabel.npyz-Storing all fiber labels (with orphans) as %sz_final_fiberslabels.npyz-Storing final fiber labels (no orphans) as %sz6Filtering tractography - keeping only no orphan fibersz_streamline_final.trkZendpoint_n_fibZendpoints_percentZintersections_percentz_statistics.matz'Saving matrix creation statistics as %sr5)=Zscipy.ioior8r9opabspathnbtrackvisreadrloadr asanyarraydataobjheaderZ get_zoomsrDsaver,rnxZ read_gpickle read_graphml TypeErrornodesrcopyGraphadd_nodertuplemeanwhereflipudr= save_fibersZmatrixZfrom_numpy_matrixZ relabel_nodesZadd_weighted_edges_fromedgesrr IndexErrorerrorrZhas_edgeZedgeZadd_edgerrZadd_nodes_fromZ remove_edgerZmedianZstd write_gpickleZto_numpy_matrixZfloat64ZsavematZint32)FZ track_fileroi_fileresolution_network_fileZ matrix_nameZmatrix_mat_name endpoint_namerPsiostatsZen_fnameZ en_fnamemmr>hdrZroir#Z roiVoxelSizer@rAr?Z fiberlabelsZfinal_fiberlabelsZfinal_fibers_idxpathnameextgpZnROIsGrKrMrZintersection_matrixr:Zfinalfibers_fnameIHr<r;ZstartROIZendROItmpZfinalfiberlengthZfinal_fibers_indicesidxZfinal_fiberlength_arrayZfinal_fiberlabels_arrayZnumfibZfibmeanZ fibmedianZfibdevrLZdiZ numfib_mlabZ numfib_dictZ fibmean_mlabZ fibmean_dictZfibmedian_mlabZfibmedian_dictZ fibdev_mlabZ fibdev_dictZintersection_matrix_nameZintersect_dictZintersection_matrix_mat_namemean_fiber_length_matrix_namemedian_fiber_length_matrix_namefiber_length_std_matrix_nameZfiberlengths_fnameZfiberlabels_fnameZfiberlabels_noorphans_fnameZout_stats_filerrrcmatst                    0   .2       "                           r~cCsX|j}g}x|D]}|j||qWt|}||d<tjd|tjj||||S)z9Stores a new trackvis file fname using only given indicesZn_countz%Writing final non-orphan fibers as %s)r`rrr8r9rTrUwrite)ZoldhdrZoldfibfnameindicesZhdrnewZ outstreamsr;Z n_fib_outrrrrgs  rgc@seZdZeddddZeddddZeddddZejddddZ edd d Z ed dd dZ edd d Z edd d Z eddd Zeddd Zeddd ZdS)CreateMatrixInputSpecTzFreesurfer aparc+aseg file)exists mandatorydesczTrackvis tract filez2Parcellation files from Connectome Mapping ToolkitFzuCounts all of the fiber-region traversals in the connectivity matrix (requires significantly more computational time)) usedefaultrz*NetworkX graph describing the connectivity)genfilerz cmatrix.matz)Matlab matrix describing the connectivityzBMatlab matrix describing the mean fiber lengths between each node.zMMatlab matrix describing the deviation in fiber lengths connecting each node.zIMatlab connectivity matrix if all region/fiber intersections are counted.z&Name for the generated endpoint arraysN)__name__ __module__ __qualname__r rl tract_filermrBoolcount_region_intersectionsout_matrix_fileout_matrix_mat_file%out_mean_fiber_length_matrix_mat_file'out_median_fiber_length_matrix_mat_file$out_fiber_length_std_matrix_mat_file out_intersection_matrix_mat_fileout_endpoint_array_namerrrrrs8  rc@seZdZedddZedddZeedddZeedddZedddZ edddZ edddZ edddZ ed ddZ ed ddZed ddZed ddZed ddZedddZedddZedddZeedddZedddZdS)CreateMatrixOutputSpecz*NetworkX graph describing the connectivityT)rrz:All of the gpickled network files output by this interfacez5All of the MATLAB .mat files output by this interfacez)Matlab matrix describing the connectivityzBMatlab matrix describing the mean fiber lengths between each node.zDMatlab matrix describing the median fiber lengths between each node.zMMatlab matrix describing the deviation in fiber lengths connecting each node.z2Saved Numpy array with the endpoints of each fiberzCSaved Numpy array with the endpoints of each fiber (in millimeters)z0Saved Numpy array with the lengths of each fiberz0Saved Numpy array with the labels for each fiberz;Saved Numpy array with the labels for each non-orphan fiberz[TrackVis file containing only those fibers originate in one and terminate in another regionz=TrackVis file containing all fibers which connect two regionszDSaved Matlab .mat file with the number of fibers saved at each stageN)rrrr matrix_fileintersection_matrix_filer matrix_filesmatlab_matrix_filesmatrix_mat_fileintersection_matrix_mat_file!mean_fiber_length_matrix_mat_file#median_fiber_length_matrix_mat_file fiber_length_std_matrix_mat_file endpoint_fileendpoint_file_mmfiber_length_filefiber_label_filefiber_labels_noorphansfiltered_tractography&filtered_tractography_by_intersectionsfiltered_tractographies stats_filerrrrrs^         rc@s0eZdZdZeZeZddZddZ ddZ dS) CreateMatrixac Performs connectivity mapping and outputs the result as a NetworkX graph and a Matlab matrix Example ------- >>> import nipype.interfaces.cmtk as cmtk >>> conmap = cmtk.CreateMatrix() >>> conmap.roi_file = 'fsLUT_aparc+aseg.nii' >>> conmap.tract_file = 'fibers.trk' >>> conmap.run() # doctest: +SKIP cCsNt|jjr.t|jj\}}}tj|d}n |jd}tj|jj}t|\}}}|dkshd}||}t|jjrtj|jj}ntj|jd} t|jj rtj|jj } ntj|jd} t|jj rtj|jj } ntj|jd} t|jj st|jj \}}}tj|}ntj|jj }t |jj |jj|jj||||jj|S)Nz.pckz.matz_mean_fiber_length.matz_median_fiber_length.matz_fiber_length_std.mat)r inputsrrrRrS_gen_outfilenamerrrrrrr~rlrmr)selfruntimerrrs_rrrtrr{rr|rr}rnrrr_run_interfaceYsJ             zCreateMatrix._run_interfacec Cs.|jj}t|jjrHt|jj\}}}tj|d}tj|d}n tj|jd}tj|jd}||d<||d<tj|jj }t|\}}}|dksd}||}||d<t|jj rtj|jj |d<ntj|jd|d<t|jj rtj|jj |d <ntj|jd |d <t|jj r>> import nipype.interfaces.cmtk as cmtk >>> rg = cmtk.ROIGen() >>> rg.inputs.aparc_aseg_file = 'aparc+aseg.nii' >>> rg.inputs.use_freesurfer_LUT = True >>> rg.inputs.freesurfer_dir = '/usr/local/freesurfer' >>> rg.run() # doctest: +SKIP The label dictionary is written to disk using Pickle. Resulting data can be loaded using: >>> file = open("FreeSurferColorLUT_adapted_aparc+aseg_out.pck", "r") >>> file = open("fsLUT_aparc+aseg.pck", "r") >>> labelDict = pickle.load(file) # doctest: +SKIP >>> labelDict # doctest: +SKIP c TCs|jj}t|\}}}tjd||tj|}tj|j }tj |}tj |} tjd| d} |jj r|jj d|_tjd|jd} nT|jj rt|jjrtj|jj|_t|j\} } }tjd| || } nd} d } t|jjrtj|jj}ntj| d |d }t|jjr.tj|jj}ntj| d |d }| r*tjd tj|jtj|jdddddddgddddd}tj |}|| krtjdtjd tjd!|i}xVtd|D]H}||d||d||d||d||dg|||d<qWtjd"tj|i}dd#gdd$gdd%gdd&gdd'gd(d)gd*d+gd,d-gd.d/gd0d1gd2d3gd4d5gd6d7gd8d9gd:d;gdd?gd@dAgdBdCgdDdEgdFdGgdHdIgdJdKgdLdMgdNdOgdPdQgdRdSgdTdUgdVdWgdXdYgdZd[gd\d]gd^d_gd`dagdbdcgdddegdfdggdhdigdjdkgdldmgdndogdpdqgdrdsgdtdugdvdwgdxdygdzd{gd|d}gdcd~gdedgdgdgdidgdmdgdodgddgddgddgdkdgddgddgddgddgddgddgddgddgddgddgddgddgddgddgddgddgddgdd0gdd2gdd4gdd6gddPgdd>gdd@gddi?i@iAiBiCiDiEi FiGiHiIiJi Ki LMNOPQRS)r+zGrey matter mask createdz Number of grey matter labels: %sZ originalIDz Label %s not in provided mappinglabelscolorsazSaving ROI File to %sz-Saving Dictionary File to %s in Pickle formatw)rrrrvrr)rrrrrr)!rrrr8r9rTrWrrXrYuniquesizerrrr rRrSrrZloadtxtrjrr,rr-Z Nifti1ImageZaffinerZr[openpickledump) rrrZ aparcpathZ aparcnameZaparcextZ niiAPARCimgZ niiAPARCdataZ niiDataLabelsZ numDataLabels write_dictprefixlutpathlutnamelutextrlrZ LUTlabelsRGBAZ numLUTLabelsZ LUTlabelDictrZmapDictZMAPPINGZniiGMmaZgreyMaskLabelsZ numGMLabelsZ labelDictZ GMlabelDictr roi_imagerCrrrrHsP                               zROIGen._run_interfacecCsx|jj}t|jjr,tj|jj|d<ntj|jd|d<t|jjr`tj|jj|d<ntj|jd|d<|S)NrlZniirZpck) _outputsrr rrrRrSrr)rrrrrr s   zROIGen._list_outputscCsft|jj\}}}|jjr d}n2|jj rNt|jjrNt|jj\}}}|}nd}|d|d|S)Nrrr.)rrrrr r)rrtrrsr)r*r+r,rrrrszROIGen._gen_outfilenameN) rrrrrrrrrrrrrrrr/sD rc Cstj}tj|}tj|}tj|j}x~|jddD]n\}}|j t |f|t tj tj tj|t |dkdd} t | d| d| d g|jt |d<q6Wtj|||S) NT)rFrGr)rrrrE)r\rar]rTrWrrXrYr_rbrrcrdrerfrk) rlrm out_filenamervrur.r#rKrMrrrr create_nodes#s   $. r2c@s6eZdZeddddZeddddZeddddZdS) CreateNodesInputSpecTzRegion of interest file)rrrz1Parcellation file from Connectome Mapping Toolkitznodenetwork.pckz/Output gpickled network with the nodes defined.)rrN)rrrr rlrmr1rrrrr35sr3c@seZdZeddZdS)CreateNodesOutputSpecz/Output gpickled network with the nodes defined.)rN)rrrr node_networkrrrrr4Csr4c@s(eZdZdZeZeZddZddZ dS) CreateNodesam Generates a NetworkX graph containing nodes at the centroid of each region in the input ROI file. Node data is added from the resolution network file. Example ------- >>> import nipype.interfaces.cmtk as cmtk >>> mknode = cmtk.CreateNodes() >>> mknode.inputs.roi_file = 'ROI_scale500.nii.gz' >>> mknode.run() # doctest: +SKIP cCs<tjdt|jj|jj|jjtjdtj|jj|S)NzCreating nodes...zSaving node network to %s) r8r9r2rrlrmr1rRrS)rrrrrrXs  zCreateNodes._run_interfacecCs"|jj}tj|jj|d<|S)Nr5)r/rrRrSrr1)rrrrrrbs zCreateNodes._list_outputsN) rrrrr3rr4rrrrrrrr6Gs   r6)F)F)*r&os.pathrrrRnumpyrZnibabelrTZnetworkxr\rZutils.filemaniprbaserrrr r r r r getLoggerr8rr*r3r=rDr~rgrrrrrrr2r3r4r6rrrrs:   (  - !9  (EMu