3
dt                 @   s  d Z ddlZddlmZ ddlmZmZmZ ddlZ	ddl
mZmZmZmZmZmZmZmZ ddlmZmZmZ ddlmZmZ dd	lmZmZ ejd
Zdd ZdddZ dddZ!G dd deZ"G dd deZ#G dd deZ$G dd deZ%G dd deZ&G dd deZ'dS )a/  
The rapidart module provides routines for artifact detection and region of
interest analysis.

These functions include:

  * ArtifactDetect: performs artifact detection on functional images

  * StimulusCorrelation: determines correlation between stimuli
    schedule and movement/intensity parameters
    N)deepcopy)loadfuncsNifti1Image   )BaseInterfacetraitsInputMultiPathOutputMultiPathTraitedSpecFileBaseInterfaceInputSpec	isdefined)ensure_list	save_jsonsplit_filename)find_indicesnormalize_mc_params)loggingconfigznipype.interfacec             C   s  |dkrddl m} || S t| |} dd }tjddddddddddddg}t| dk rvtj| |t| d f} t| f| _tjd	}| dd
 |dd
df< tjd	}|| d
 |dd
dd
f< tjd	}|| d	 j	 |d< tjd	}|| d |ddddf< tjd	}	tj
| dd |	dd
dd
f< tjd	}
| dd |
d< |dkrtj|tj|tj|tj|tj|	|
S tj|tj|tj|tj|tj|	|
S )zReturn affine matrix given a set of translation and rotation parameters

    params : np.array (upto 12 long) in native package format
    source : the package that generated the parameters
             supports SPM, AFNI, FSFAST, FSL, NIPY
    ZNIPYr   )to_matrix44c             S   s0   t jt j| t j| gt j|  t j| ggS )N)nparraycossin)x r   </tmp/pip-build-7vycvbft/nipype/nipype/algorithms/rapidart.py<lambda>5   s    z$_get_affine_matrix.<locals>.<lambda>      N      r         	   AFNIFSFASTr   r   r   r   r   r   r   r   )r)   r*   r   r   r   r   r   r   )r+   r,   )r&   r'   )Znipy.algorithms.registrationr   r   r   r   lenhstackshapeZeyeraveldiagdot)paramssourcer   ZrotfuncqTZRxZRyZRzSZShr   r   r   _get_affine_matrix&   s0    
"




"

,r8   c                s*    fddt  jd D }t|||S )a  Calculates the maximum overall displacement of the midpoints
    of the faces of a cube due to translation and rotation.

    Parameters
    ----------
    mc : motion parameter estimates
        [3 translation, 3 rotation (radians)]
    use_differences : boolean
    brain_pts : [4 x n_points] of coordinates

    Returns
    -------

    norm : at each time point
    displacement : euclidean distance (mm) of displacement at each coordinate

    c                s"   g | ]}t  |d d f qS )N)r8   ).0i)mcr4   r   r   
<listcomp>b   s    z_calc_norm.<locals>.<listcomp>r   )ranger/   _calc_norm_affine)r;   use_differencesr4   	brain_ptsaffinesr   )r;   r4   r   
_calc_normO   s    rB   c             C   s  |dkrLt jdddg}t jdddg}t jt j||ft jdf}d}n|}|j|jd  }t jt| |f}|dk	rt jt| t	|d f}xt
| D ]\}	}
t j|
|d	dddf j ||	ddf< |dk	rt jt jt jt j||	ddf d|jd f|d	dddf  d
d	d||	ddf< qW t jt| }|rt jt jd|ft j|dd	dfd	d}xt|jd	 D ]P}	t jt jt jt jt jt j||	ddf d
d|jd fd	d||	< qW n<d	dlm} t j||d	dd}t jt jt j|d
dd}||fS )a  Calculates the maximum overall displacement of the midpoints
    of the faces of a cube due to translation and rotation.

    Parameters
    ----------
    affines : list of [4 x 4] affine matrices
    use_differences : boolean
    brain_pts : [4 x n_points] of coordinates

    Returns
    -------

    norm : at each time point
    displacement : euclidean distance (mm) of displacement at each coordinate

    NF   K   n   -   r   r$   r"   r   r   )axis)nrG   )detrendZconstant)rG   typeiii)r   r$   )r   r1   vstackr.   onessizer/   zerosr-   int	enumerater2   r0   sqrtsumpowerZreshapeconcatenatediffr=   maxabsZscipy.signalrI   mean)rA   r?   r@   ZresposZresnegZall_ptsdisplacementZn_ptsZnewposr:   affineZnormdatarI   r   r   r   r>   f   sF    , $*r>   c            	   @   s@  e Zd ZeedddddZeedddddZejddd	d
ddddZ	ej
ddgdddddZejdddgddZejddgdddZejddgddZejddgddZejdddZejddddddZedd d!Zejd"d#Zejddd$d%Zejdd&dd'Zejd(d)d*d+d,dd'Zejdd-dd'Zejd.d/dd'Zd0S )1ArtifactDetectInputSpecT)existsz(Names of realigned functional data files)desc	mandatoryzJNames of realignment parameters corresponding to the functional data files)r^   r]   SPMZFSLr&   ZNiPyr'   zSource of movement parametersFr   zUse differences between successive motion (first element) and intensity parameter (second element) estimates in order to determine outliers.  (default is [True, False]))Zminlenmaxlen
usedefaultr]   norm_thresholdzIUses a composite of the motion parameters in order to determine outliers.)ra   requiresr]   rotation_thresholdtranslation_thresholdzVThreshold to use to detect motion-related outliers when composite motion is being used)xorr^   r]   zAThreshold (in radians) to use to detect rotation-related outliers)r^   rf   r]   z?Threshold (in mm) to use to detect translation-related outlierszHIntensity Z-threshold use to detection images that deviate from the mean
spm_globalfilethreshaR  Type of mask that should be used to mask the functional data. *spm_global* uses an spm_global like calculation to determine the brain mask. *file* specifies a brain mask file (should be an image file consisting of 0s and 1s). *thresh* specifies a threshold to use. By default all voxels are used,unless one of these mask types are definedz,Mask file to be used if mask_type is 'file'.)r\   r]   z3Mask threshold to be used if mask_type is 'thresh'.)r]   z2Intersect the masks when computed from spm_global.)ra   r]   zsave plots containing outliers)r]   ra   ZpngsvgepsZpdfzfile type of the outlier plotzyuse the brain mask to determine bounding boxfor composite norm (worksfor SPM and Nipy - currentlyinaccurate for FSL, AFNIg       @z4use this threshold when mask type equal's spm_globalN)__name__
__module____qualname__r	   r   realigned_filesrealignment_parametersr   Enumparameter_sourceZListBoolr?   Booluse_normZFloatrb   rd   re   zintensity_threshold	mask_type	mask_filemask_thresholdintersect_mask	save_plot	plot_typebound_by_brainmaskglobal_thresholdr   r   r   r   r[      s   
r[   c               @   sr   e Zd ZeeddddZeeddddZeeddZeeddddZeeddZ	eed	dZ
eed
dZdS )ArtifactDetectOutputSpecT)r\   zfOne file for each functional run containing a list of 0-based indices corresponding to outlier volumes)r]   zeOne file for each functional run containing the global intensity values determined from the brainmaskz>One file for each functional run containing the composite normzOne file for each functional run containing information about the different types of artifacts and if design info is provided then details of stimulus correlated motion and a listing or artifacts by event type.zGOne image file for each functional run containing the detected outliersz]One image file for each functional run containing the mask used for global signal calculationzSOne image file for each functional run containing the voxel displacement timeseriesN)rl   rm   rn   r
   r   outlier_filesintensity_files
norm_filesstatistic_files
plot_files
mask_filesdisplacement_filesr   r   r   r   r~   )  s(   
	r~   c                   sR   e Zd ZdZeZeZ fddZdd Z	dd Z
dd	 ZdddZdd Z  ZS )ArtifactDetecta,  Detects outliers in a functional imaging series

    Uses intensity and motion parameters to infer outliers. If `use_norm` is
    True, it computes the movement of the center of each face a cuboid centered
    around the head and returns the maximal movement across the centers. If you
    wish to use individual thresholds instead, import `Undefined` from
    `nipype.interfaces.base` and set `....inputs.use_norm = Undefined`


    Examples
    --------

    >>> ad = ArtifactDetect()
    >>> ad.inputs.realigned_files = 'functional.nii'
    >>> ad.inputs.realignment_parameters = 'functional.par'
    >>> ad.inputs.parameter_source = 'FSL'
    >>> ad.inputs.norm_threshold = 1
    >>> ad.inputs.use_differences = [True, False]
    >>> ad.inputs.zintensity_threshold = 3
    >>> ad.run()  # doctest: +SKIP
    c                s   t t| jf | d S )N)superr   __init__)selfinputs)	__class__r   r   r   |  s    zArtifactDetect.__init__c             C   s  t |ttfr|}nt |tr(|d }ntdt|\}}}tjj|djd|df}tjj|djd|df}tjj|djd|df}	tjj|djd	|df}
tjj|djd
|d| j	j
f}tjj|djd||f}tjj|djd||f}|||	|
|||fS )a  Generate output files based on motion filenames

        Parameters
        ----------

        motionfile: file/string
            Filename for motion parameter file
        output_dir: string
            output directory in which the files will be generated
        r   zUnknown type of file zart.z_outliers.txtzglobal_intensity.z.txtzstats.znorm.zplot..zdisp.zmask.)
isinstancestrbyteslist	Exceptionr   ospathjoinr   r{   )r   
motionfile
output_dirinfile_filenameextartifactfileintensityfile	statsfilenormfileplotfiledisplacementfilemaskfiler   r   r   _get_output_filenames  s.    

z$ArtifactDetect._get_output_filenamesc             C   sX  | j  j }g |d< g |d< g |d< g |d< t| jjrX| jjrXg |d< | jjrXg |d< t| jjrt| jjrtg |d< xtt| jj	D ]\}}| j
|tj \}}}}}}	}
|d j|| |d j|| |d j|| |d j||
 t| jjo| jjr,|d j|| | jjr,|d j||	 t| jjr| jjr|d j|| qW |S )Nr   r   r   r   r   r   r   )_outputsgetr   r   rt   r|   rz   rP   r   ro   r   r   getcwdinsert)r   outputsr:   fZoutlierfiler   r   r   r   r   r   r   r   r   _list_outputs  s0    	
zArtifactDetect._list_outputsc             C   s   dd l }|jtjdd dd lj}|j| |j|j |j	 g |j
dt|d g t|r|jtj|d d d f djtj|j |j	 gt|dfjd |jd |j| d S )	Nr   	executionmatplotlib_backendr   r   rzScans - 0-based)r   r   )
matplotlibuser   r   matplotlib.pyplotpyplotZplotZylimminrV   Zxlimr-   r   Ztiler6   ZxlabelZylabel)r   Zwaveoutliersnamer   pltr   r   r   _plot_outliers_with_wave  s    

 
z'ArtifactDetect._plot_outliers_with_waveNc       5      C   s~  ddl m} |stj }t|ttfr0t|}n<t|trlt	|dkrTt|d }ndd |D }t
j|}|j\}}	}
}|jtjd}|j}tj|df}| jj}|dkr>tjd | jj}|rtj||	|
ftd}xJt|D ]>}|d	d	d	d	d	d	|f }|tj|| jj k}|| }qW x@t|D ]4}|d	d	d	d	d	d	|f }tj|| ||< q0W t	t|tj||	|
fd
 k rd}tj|df}|s:tjd tj||	|
|f}x|t|D ]p}|d	d	d	d	d	d	|f }|tj|| jj k}||d	d	d	d	d	d	|f< tj|| tj| ||< qW n|dkrt| jj }|jtjd}|j}|dk}xt|D ]4}|d	d	d	d	d	d	|f }tj|| ||< qzW n|dkrxzt|D ]@}|d	d	d	d	d	d	|f }|| jj!k}tj|| ||< qW n,tj||	|
f}tj||dkd	d	f d}|j"|dd}| jj#d rztj$tjd/tj%|dddfdd}|tj&| tj'| }tt(|| jj)k}tj*|}t+|}| j,||\}}}}}}} t-|j.tj/|}!|!j0|  | jj1rd	}"| jj2r^tj3|}#tj4|#d tj4|#d |#d ffj5}$tj6|tj7|$tj|$jd dffj5}"t8|| jj#d | jj9|"d\}%}&t|%| jj:k}'t|%dk }(|&d	k	rtj||	|
|ftj;d})x:t|D ].}*|&|*d	d	f |)|#d |#d |#d |*f< qW t-|)|}+|+j0| n| jj#d r@tj$tjd0tj%|dddfdd}|d	d	ddf },|d	d	ddf }-ttj<t(|,| jj=kddk}'ttj<t(|-| jj>kddk}(tj?tj@|tj@|'|(}.tjA||.ddd tjA||ddd | jj1rtjA||%ddd tB| jjCr| jjCrdd	lD}/|/jEtFjGdd dd	lHjI}0|0jJ }1tB| jj1rj| jj1rj|0jKd n
|0jKd | jL||d tB| jj1r| jj1r|0jKd  | jL|%tj@|'|(d! nNd"}2| jj#d rd#}2|0jKd$ | jL|,|'d%|2  |0jKd& | jL|-|(d'|2  |0jM| |0jN|1 tj@|'|(}3||d(t	tjO||3t	tjP||3t	tjP|3|d)d*d+| jj#d itj&|ddjQ tjR|ddjQ tjS|ddjQ tj'|ddjQ d,gid-d+| jj#d itj&|ddjQ tjR|ddjQ tjS|ddjQ tj'|ddjQ d,gig}4| jj1rp|4jTdd.tj&|%ddjQ tjR|%ddjQ tjS|%ddjQ tj'|%ddjQ d,i tU||4 d	S )1z5
        Core routine for detecting outliers
        r   )signalr   c             S   s   g | ]}t |qS r   )r   )r9   r   r   r   r   r<     s    z8ArtifactDetect._detect_outliers_core.<locals>.<listcomp>)Zdtyperg   zart: using spm globalN
   Fznot intersect_mask is Truerh   g      ?ri   )rG   )rH   rG   r   )r@   r$   r"   s   %d )fmt	delimiters   %.2fs   %.4fr   r      i7  Z	Intensity   z	Norm (mm)r   rU   i8  zTranslation (mm)i9  zRotation (rad))Zmotion_fileZfunctional_file)Zcommon_outliersZintensity_outliersmotion_outliersZmotionzusing differences)rX   r   rV   stdZ	intensityZmotion_norm)r   r   )r   r$   )VZscipyr   r   r   r   r   r   r   r   r-   r   Zconcat_imagesr/   Z	get_fdatar   Zfloat32rZ   rN   r   rv   ifloggerdebugry   rL   boolr=   Znanmeanr}   r   prodinfoZnansumrw   rx   rI   r?   rT   rU   rX   r   rW   ru   loadtxtr   r   r   ZastypeZuint8to_filenamert   r|   ZnonzerorK   r6   r2   r.   rB   rr   rb   Zfloat64rR   re   rd   uniqueZunion1dZsavetxtr   rz   r   r   r   r   r   r   figureZsubplotr   ZsavefigcloseZintersect1dZ	setdiff1dtolistr   rV   r   r   )5r   Zimgfiler   Zrunidxcwdr   ZnimZimagesr   yzZ
timepointsdatarZ   gZmasktypery   maskt0ZvolZmask_tmpZmaskimggzZiidxmc_inr;   r   r   r   r   r   r   r   Zmask_imgr@   Zvoxel_coordsZcoordsZnormvalrY   ZtidxZridxZdmapr:   ZdimgZtravalZrotvalr   r   r   ZfigrU   r   statsr   r   r   _detect_outliers_core  s$   




"
$

$





"$

.
 








z$ArtifactDetect._detect_outliers_corec             C   sL   t | jj}t | jj}x.t|D ]"\}}| j||| |tj d q"W |S )zExecute this module.)r   )r   r   ro   rp   rP   r   r   r   )r   runtimeZfuncfilelistmotparamlistr:   Zimgfr   r   r   _run_interface  s
    zArtifactDetect._run_interface)N)rl   rm   rn   __doc__r[   
input_specr~   output_specr   r   r   r   r   r   __classcell__r   r   )r   r   r   b  s   )"
 Ur   c               @   sP   e Zd ZeedddddZeedddddZeddddZej	dddZ
d	S )
StimCorrInputSpecT)r\   zJNames of realignment parameters corresponding to the functional data files)r^   r]   z(Name of file containing intensity valuesz,SPM mat file (use pre-estimate SPM.mat file))r\   r^   r]   z9state if the design matrix contains concatenated sessionsN)rl   rm   rn   r	   r   rp   intensity_valuesspm_mat_filer   rs   concatenated_designr   r   r   r   r     s   r   c               @   s   e Zd ZeeddddZdS )StimCorrOutputSpecT)r\   z+List of files containing correlation values)r]   N)rl   rm   rn   r
   r   stimcorr_filesr   r   r   r   r     s   r   c               @   sD   e Zd ZdZeZeZdd ZdddZ	dddZ
d	d
 Zdd ZdS )StimulusCorrelationa   Determines if stimuli are correlated with motion or intensity
    parameters.

    Currently this class supports an SPM generated design matrix and requires
    intensity parameters. This implies that one must run
    :ref:`ArtifactDetect <nipype.algorithms.rapidart.ArtifactDetect>`
    and :ref:`Level1Design <nipype.interfaces.spm.model.Level1Design>` prior to
    running this or provide an SPM.mat file and intensity parameters through
    some other means.

    Examples
    --------

    >>> sc = StimulusCorrelation()
    >>> sc.inputs.realignment_parameters = 'functional.par'
    >>> sc.inputs.intensity_values = 'functional.rms'
    >>> sc.inputs.spm_mat_file = 'SPM.mat'
    >>> sc.inputs.concatenated_design = False
    >>> sc.run() # doctest: +SKIP

    c             C   s>   t jj|\}}t jj|\}}t jj|djd|df}|S )a  Generate output files based on motion filenames

        Parameters
        ----------
        motionfile: file/string
            Filename for motion parameter file
        output_dir: string
            output directory in which the files will be generated
        r   zqa.z_stimcorr.txt)r   r   splitsplitextr   )r   r   r   r   r   corrfiler   r   r   r     s    
z)StimulusCorrelation._get_output_filenamesNc             C   s8  |st j }tj|}tj|}|jd df|_|jd }|jd }tjtj||f|f}	tj|	dd}
| j||}t|d}|j	d |j	d|  xVt
|D ]J}|j	d|  x,|
||tj| f D ]}|j	d|  qW |j	d	 qW |j	d
|  x,t
|D ] }|j	d||
|df f  qW |j  dS )zD
        Core routine for determining stimulus correlation

        r   r   )ZrowvarwzStats for:
zStimulus correlated motion:
%s
zSCM.%d:z %.2f
z"Stimulus correlated intensity:
%s
zSCI.%d: %.2f
Nr(   )r   r   r   r   r/   r.   Zcorrcoefr   openwriter=   aranger   )r   r   r   designmatrixr   r   Zg_inZdcolZmccolZconcat_matrixcmr   rh   r:   vr   r   r   _stimcorr_core   s,    





 z"StimulusCorrelation._stimcorr_corec             C   s   |d d d j d d j}|d d d jd | jd }|dkrh|d d d jd | jd d }|d d d jd | jd ttt| d }|j	|j
 ddj	|j
 dd}|S )z
        Parameters
        ----------
        spmmat: scipy matlab object
            full SPM.mat file loaded into a scipy object
        sessidx: int
            index to session that needs to be extracted.
        r_   r   Nr   )rG   )xXXZSessUrowcolr   r=   r-   Ztaker   )r   spmmatsessidxrowsr   r   colsZ	outmatrixr   r   r   _get_spm_submatrix  s    	 $4 z&StimulusCorrelation._get_spm_submatrixc             C   s   ddl j}| jj}| jj}|j| jjdd}g }xtt|D ]z}|}d}	| jj	rd}t
j|| }
t
j|t
j|
jd  }	|j|
jd  | j|||	}| j|| || |tj  q>W |S )zExecute this module.r   NF)Zstruct_as_record)Zscipy.ioior   rp   r   Zloadmatr   r=   r-   r   r   r   rR   r   r/   appendr   r   r   r   )r   r   sior   Zintensityfilesr   Znrowsr:   r   r   r   Zmatrixr   r   r   r   /  s     
 z"StimulusCorrelation._run_interfacec             C   sR   | j  j }g }x0t| jjD ] \}}|j|| j|tj  qW |rN||d< |S )Nr   )	r   r   rP   r   rp   r   r   r   r   )r   r   filesr:   r   r   r   r   r   C  s    z!StimulusCorrelation._list_outputs)N)N)rl   rm   rn   r   r   r   r   r   r   r   r   r   r   r   r   r   r   r     s   

r   )N)N)(r   r   copyr   Znibabelr   r   r   numpyr   Zinterfaces.baser   r   r	   r
   r   r   r   r   Zutils.filemanipr   r   r   Z
utils.miscr   r   r   r   r   	getLoggerr   r8   rB   r>   r[   r~   r   r   r   r   r   r   r   r   <module>   s*   (

)

C 9  Y