minc-library.bpipe

BEASTMODEL_DIR="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2"
BEASTMODEL_NAME="mni_icbm152_t1_tal_nlin_sym_09c"
BEASTLIBRARY_DIR="${System.getenv().QUARANTINE_PATH}/resources/BEaST_libraries/combined"

//Model which has a facemask, must correspond to the same space as MNIMODEL
DEFACEMODEL="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c"

//We need a resample model so beast priors work (they're already in this space)
RESAMPLEMODEL="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c.mnc"

//Specify target linear registration model here options are
//sym or asym
//t1, t2, or pd
//Models older than icbm_2009c are also available, see resources directory
//Choice of model affects MNI_register, cutneck, and deface
REGISTRATIONMODEL="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c.mnc"
//REGISTRATIONMODEL="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_asym_09c_minc2/mni_icbm152_t1_tal_nlin_asym_09c.mnc"
//REGISTRATIONMODEL="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t2_tal_nlin_sym_09c.mnc"
//REGISTRATIONMODEL="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_asym_09c_minc2/mni_icbm152_t2_tal_nlin_asym_09c.mnc"

//Optional second registration model for two-model registration
//ADNI MNI MODEL
//REGISTRATIONMODEL="${System.getenv().QUARANTINE_PATH}/resources/mni_adni/mni_adni_t1w_tal_nlin_asym.mnc"
//PEDIATRIC MNI MODEL
//REGISTRATIONMODEL="${System.getenv().QUARANTINE_PATH}/resources/nihpd_sym_all_minc2/nihpd_sym_04.5-18.5_t1w.mnc"

//Mask for masked registration
REGISTRATIONBRAINMASK="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c_mask.mnc"
REGISTRATIONHEADMASK="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c_headmask.mnc"

//Mask of white matter for N4 improvement
REGISTRATIONWMMASK="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_wm_tal_nlin_sym_09c.mnc"
REGISTRATIONGMMASK="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_gm_tal_nlin_sym_09c.mnc"

REGISTRATIONOUTLINE="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c_outline.mnc"
REGISTRATIONANTIMASK="${System.getenv().QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c_antimask.mnc"

//Model to normalize intensity against
//NORMALIZEMODEL="/path/to/model.mnc"

convert = {

  if(file(input).name.endsWith('.mnc')) {
      transform('.mnc') to('.convert.mnc') {
        exec """
            mincconvert -2 -clobber -compress 9 $input.mnc $output.mnc
        """
      }
  }
  else
  if(file(input).name.endsWith('.nii.gz')) {
      transform('.nii.gz') to('.convert.mnc') {
      exec """
          nii2mnc -clobber $input.nii.gz $output.mnc
      """
      }
  }
  else
  if(file(input).name.endsWith('.nii')) {
      transform('.nii') to('.convert.mnc') {
      exec """
          nii2mnc -clobber $input.nii $output.mnc
      """
      }
  }
  branch.nativemnc = "${branch.name}".minus(".mnc").minus(".nii.gz").minus(".nii").concat(".convert")
}

n3correct = {
    //Runs the nu_correct (aka n3correct) command with optimal settings according to
    //http://www.ncbi.nlm.nih.gov/pubmed/19559796
    //Requires minc-toolkit
    exec "nu_correct -verbose -clobber -iter 200 -shrink 2 -stop 0.000001 -fwhm 0.05 -distance 30 $input.mnc $output.mnc"
    branch.nativemnc = "${branch.nativemnc}.n3correct"
}

nlm_denoise = {
    //Runs non-local-means filter on data
    //http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=4359947
    //This is probably only valid for very noise 3T, or most 1.5T data
    //Requires minc-toolkit
    exec "volume_denoise.pl $input.mnc $output.mnc --verbose --clobber"
    branch.nativemnc = "${branch.nativemnc}.nlm_denoise"
}

anlm_denoise = {
    //Runs the "improved" adaptive non local means filter
    //http://onlinelibrary.wiley.com/doi/10.1002/jmri.22003/full
    //This is probably only valid for very noise 3T, or most 1.5T data
    //Requires minc-toolkit
    exec "minc_anlm --rician --verbose $input.mnc $output.mnc"
    branch.nativemnc = "${branch.nativemnc}.anlm_denoise"
}

n4correct = {
    //Runs the improved n3correct aka N4correct
    //http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=5445030
    uses(GB:6,threads:4) {
    exec """
      iterativeN4_multispectral.sh $output.mnc $input.mnc
    ""","n4correct"
    }
    branch.nativemnc = "${branch.nativemnc}.n4correct"
}

n4correctsimple = {
    //Runs the improved n3correct aka N4correct
    //http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=5445030
    //Simple n4 done with a whole-scan mask for when internal otsu masking fails
    exec """
      tmpdir=`mktemp -d`

      minccalc -unsigned -byte -expression '1' $input.mnc \$tmpdir/mask.mnc
      ImageMath 3 \$tmpdir/weight.mnc ThresholdAtMean $input.mnc 1

      N4BiasFieldCorrection -d 3 -s 2 -b [200] -c [300x300x300x300,0.0] -i $input.mnc -o $output.mnc -x \$tmpdir/mask.mnc -w \$tmpdir/weight.mnc --histogram-sharpening [0.05,0.01,200]
    ""","n4correct"
    branch.nativemnc = "${branch.nativemnc}.n4correctsimple"
}


dumb_normalize = {
    //Performs a dumb normalization of the data, throwing away top and bottom 1%
    //Data is then scaled 0-65535 and stored in a sensible unsigned int format
    //Requires minc-toolkit
    exec """
      minccalc -clobber -expression "A[0]<(0.01*\$(mincstats -quiet -max $input.mnc)) || A[0]>(0.99*\$(mincstats -quiet -max $input.mnc)) ? 0 : A[0]" $input.mnc $output.tmp
      minccalc -clobber -expression "A[0]/\$(mincstats -quiet -max $output.tmp)*65535" $output.tmp $output.mnc
    """
}

nuyl_normalize = {
    //Performs histogram landmark normalization to a model using minc_nuyl
    //Requires NORMALIZEMODEL to be set to a reference mnc (atlas perhaps?)
    //See http://www.ncbi.nlm.nih.gov/pubmed/10784285
    exec """
      minc_nuyl --verbose --fix_zero_padding --source-mask ${NORMALIZEMODELMASK} --target_mask ${branch.nativemnc}.beastnormalize.mnc $input.mnc ${NORMALIZEMODEL} $output.mnc
    """
}

pol_normalize = {
    //Performs normalization against a model using volume_pol
    //Requires NORMALIZEMODEL to be set to a reference mnc (atlas perhaps?)
    //
    exec """
      volume_pol --verbose --clobber --order 1 --source_mask ${NORMALIZEMODELMASK} --target_mask ${branch.nativemnc}.beastnormalize.mnc $input.mnc ${NORMALIZEMODEL} $output.mnc
    """
}

cutneckapplyautocrop = {
    //Alternative implementation of cutneck which applies via autocrop in native space
    forward input.mnc

    //Transforms model mask to native space
    exec """
      itk_resample --clobber --labels --byte
        --like ${branch.nativemnc}.mnc
        --transform $input.xfm
        ${REGISTRATIONHEADMASK} ${output("${branch.nativemnc}.headmask.mnc")} &&
      minccalc  -clobber -expression 'A[0]*A[1]' ${branch.nativemnc}.mnc ${output("${branch.nativemnc}.headmask.mnc")} ${output("${branch.nativemnc}.cutneck.mnc")} &&
      autocrop -clobber -isoexpand 10mm -bbox ${output("${branch.nativemnc}.headmask.mnc")} ${output("${branch.nativemnc}.cutneck.mnc")} ${output("${branch.nativemnc}.cutneckapplyautocrop.mnc")}
    """
    branch.nativemnc = "${branch.nativemnc}.cutneckapplyautocrop"
}

linear_bestlinreg = {
    //Two target bestlinreg
    exec """
      tmpdir=`mktemp -d` &&
      bestlinreg_g -noverbose -clobber -nmi -lsq12 -target_mask ${REGISTRATIONBRAINMASK} $input.mnc ${REGISTRATIONMODEL} ${output.xfm.prefix}_inverse.xfm &&
      xfminvert -clobber ${output.xfm.prefix}_inverse.xfm $output.xfm &&
      bestlinreg_g -noverbose -clobber -nmi -lsq12 -source_mask ${REGISTRATIONANTIMASK} ${REGISTRATIONMODEL} $input.mnc ${output.xfm.prefix}_ICV.xfm &&
      itk_resample --clobber --order 3
        --like ${RESAMPLEMODEL}
        --transform ${output.xfm.prefix}_inverse.xfm
        $input.mnc \$tmpdir/transformed.mnc &&
      mincmath -clobber -clamp -const2 0 \$(mincstats -quiet -max \$tmpdir/transformed.mnc) \$tmpdir/transformed.mnc $output.mnc &&
      rm -r \$tmpdir
    """
}

linear_antsRegistration = {
    exec """
      tmpdir=`mktemp -d` &&
      antsRegistration --dimensionality 3 --verbose --minc \
        --output [\$tmpdir/reg] \
        --use-histogram-matching 1 \
--initial-moving-transform [ ${REGISTRATIONMODEL},$input.mnc,1 ] \
--transform Translation[ 0.1 ] \
        --metric Mattes[ ${REGISTRATIONMODEL},$input.mnc,1,32,None ] \
        --convergence [ 500x500x500x500x500x500x500x500,1e-6,10 ] \
        --shrink-factors 7x7x7x7x7x7x7x7 \
        --smoothing-sigmas 6.78129076418x6.35574237559x5.93006674681x5.50423435717x5.07820577132x4.65192708599x4.22532260674x3.79828256043mm \
        --masks [ NOMASK,NOMASK ] \
--transform Rigid[ 0.1 ] \
        --metric Mattes[ ${REGISTRATIONMODEL},$input.mnc,1,32,None ] \
        --convergence [ 500x500x500x500x500x500x500x500,1e-6,10 ] \
        --shrink-factors 7x7x7x7x7x7x6x5 \
        --smoothing-sigmas 5.07820577132x4.65192708599x4.22532260674x3.79828256043x3.37064139994x2.94213702015x2.51232776601x2.08040503813mm \
        --masks [ NOMASK,NOMASK ] \
--transform Similarity[ 0.1 ] \
        --metric Mattes[ ${REGISTRATIONMODEL},$input.mnc,1,32,None ] \
        --convergence [ 500x500x500x500x500x450x150,1e-6,10 ] \
        --shrink-factors 7x7x6x5x4x3x2 \
        --smoothing-sigmas 3.37064139994x2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037mm \
        --masks [ NOMASK,NOMASK ] \
--transform Similarity[ 0.1 ] \
        --metric Mattes[ ${REGISTRATIONMODEL},$input.mnc,1,32,None ] \
        --convergence [ 500x500x500x500x500x450x150,1e-6,10 ] \
        --shrink-factors 7x7x6x5x4x3x2 \
        --smoothing-sigmas 3.37064139994x2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037mm \
        --masks [ ${REGISTRATIONBRAINMASK},NOMASK ] \
--transform Affine[ 0.1 ] \
        --metric Mattes[ ${REGISTRATIONMODEL},$input.mnc,1,64,None ] \
        --convergence [ 500x450x150x50,1e-6,10 ] \
        --shrink-factors 4x3x2x1 \
        --smoothing-sigmas 1.64470459404x1.20112240879x0.735534255037x0.0mm \
        --masks [ ${REGISTRATIONBRAINMASK},NOMASK ] &&
      xfminvert \$tmpdir/reg0_GenericAffine.xfm ${output.xfm.prefix}_inverse.xfm &&
      mv -f \$tmpdir/reg0_GenericAffine.xfm $output.xfm &&
      bestlinreg_g -noverbose -clobber -nmi -lsq12 -source_mask ${REGISTRATIONANTIMASK} ${REGISTRATIONMODEL} $input.mnc ${output.xfm.prefix}_ICV.xfm &&
      itk_resample --clobber --order 3
        --like ${RESAMPLEMODEL}
        --transform ${output.xfm.prefix}_inverse.xfm
        $input.mnc \$tmpdir/transformed.mnc &&
      mincmath -clobber -clamp -const2 0 \$(mincstats -quiet -max \$tmpdir/transformed.mnc) \$tmpdir/transformed.mnc $output.mnc &&
      rm -r \$tmpdir
    """
}

resample_to_lsq6_space = {
    //Must be run as the last stage, resamples all the useful files (brains and masks) to lsq6 space with the target model
    output.dir="lsq6"

    //Clever descale command, extract the scale and shear values from lsq12 registration, invert, and apply to lsq12, result is a lsq6 which is pretty nice
    exec """
      tmpdir=`mktemp -d` &&
      param2xfm \$(xfm2param ${input.xfm.prefix}_inverse.xfm | grep -E 'scale|shear') \$tmpdir/${branch.nativemnc}.scale.xfm &&
      xfminvert \$tmpdir/${branch.nativemnc}.scale.xfm \$tmpdir/${branch.nativemnc}.unscale.xfm &&
      xfmconcat ${input.xfm.prefix}_inverse.xfm \$tmpdir/${branch.nativemnc}.unscale.xfm ${output("${branch.nativemnc}.lsq6.xfm")} &&
      mincresample -clobber -tfm_input_sampling
        -transform ${output("${branch.nativemnc}.lsq6.xfm")}
        ${branch.nativemnc}.mnc ${output("${branch.nativemnc}.lsq6.mnc")} &&
      mincresample -clobber -tfm_input_sampling -byte -unsigned -keep_real_range
        -transform ${output("${branch.nativemnc}.lsq6.xfm")}
        ${branch.nativemnc}.beastmask.mnc ${output("${branch.nativemnc}.beastmask.lsq6.mnc")} &&
      mincresample -clobber -tfm_input_sampling
        -transform ${output("${branch.nativemnc}.lsq6.xfm")}
        ${branch.nativemnc}.beastextract.mnc ${output("${branch.nativemnc}.beastextract.lsq6.mnc")}
    """
}

beastnormalize = {
    //Linearly scales the intensities to the range [0;100] using 0.1%-99.9%
    //of the voxels in the intensity histogram
    //No registration or n3 since this is done in prior steps
    branch.modelmnc = "$input.mnc.prefix"
    exec """
      tmpdir=`mktemp -d`;
      minc_anlm $input.mnc \$tmpdir/${branch.nativemnc}.denoise.mnc &&
      volume_pol --order 1 --min 0 --max 100 --noclamp  \$tmpdir/${branch.nativemnc}.denoise.mnc ${RESAMPLEMODEL} --source_mask ${REGISTRATIONBRAINMASK} --target_mask ${REGISTRATIONBRAINMASK} --clobber $output.mnc
    """
}

beastmask = {
    //Generates a brain mask for subject brain in model space
    uses(GB:6,threads:4) {
    exec "mincbeast -clobber -verbose -fill -median -same_res -flip -conf ${BEASTLIBRARY_DIR}/default.1mm.conf ${BEASTLIBRARY_DIR} $input.mnc $output.mnc"
    }
}

beastextract = {
    //Applies the mask from beastmask to the original input file to extract the brain

    //Apply brain mask to original registered brain
    exec """
      minccalc -verbose -clobber -expression 'A[0]*A[1]' ${branch.modelmnc}.mnc $input.mnc ${output("${branch.modelmnc}.beastextract.mnc")}
    """

    forward input.mnc

    //Transform brain mask to native space
    exec """
      tmpdir=`mktemp -d` &&
      itk_resample --clobber --labels --byte
        --like ${branch.nativemnc}.mnc
        --transform $input.xfm
        $input.mnc ${output("${branch.nativemnc}.beastmask.mnc")} &&
      minccalc -unsigned -verbose -clobber -expression 'A[0]*A[1]' ${branch.nativemnc}.mnc ${output("${branch.nativemnc}.beastmask.mnc")} \$tmpdir/extracted.mnc &&
      autocrop -isoexpand 10mm -bbox \$tmpdir/extracted.mnc \$tmpdir/extracted.mnc ${output("${branch.nativemnc}.beastextract.mnc")} &&
      rm -r \$tmpdir
    """

    //Generate a brainmatter mask for volume estimations
    exec """
      tmpdir=`mktemp -d` &&
      ThresholdImage 3 ${branch.nativemnc}.mnc \$tmpdir/otsu.mnc Otsu 1 &&
      ImageMath 3 ${output("${branch.nativemnc}.brainvolume.mnc")} m \$tmpdir/otsu.mnc ${output("${branch.nativemnc}.beastmask.mnc")} &&
      rm -r \$tmpdir
    """

}

beast = segment {
    //Runs the beast brain extraction toolchain
    //Requires models defined in $BEASTMODEL_DIR, $BEASTMODEL_NAME and a library in $BEASTLIBRARY_DIR
    //Requires minc-toolkit
    //http://www.ncbi.nlm.nih.gov/pubmed/21945694
    beastnormalize + beastmask + beastextract
}

defaceapply = {
    //Applies the inverse face mask to the MRI volume
    //Expects the brain to be in MNI space
    exec """
      minccalc -verbose -clobber -expression 'A[0]*(1 - A[1])' $input.mnc ${DEFACEMODEL}_face_mask_res.mnc $output.mnc
    """

    exec """
      itk_resample --clobber --labels --byte
        --like ${branch.nativemnc}.mnc
        --transform $input.xfm
        ${DEFACEMODEL}_face_mask_res.mnc ${output("${branch.nativemnc}.facemask.mnc")} &&
        minccalc -verbose -clobber -expression 'A[0]*(1 - A[1])' ${branch.nativemnc}.mnc ${output("${branch.nativemnc}.facemask.mnc")} ${output("${branch.nativemnc}.defaceapply.mnc")}
    """
}

deface = segment {
    //Linear registers brain to MNI space
    //Applies inverse facemask to MRI data
    linear_bestlinreg + defaceapply
}

cutneck = segment {
    //Linearly registers brain to MNI space
    //Applies neckmask to MRI data
    linear_bestlinreg + cutneckapplyautocrop
}

clean_and_center = {
    //Fixes a bunch of possible broken bits in MINC files
    //Requires minc-toolkit-extras
    exec "clean_and_center_minc.pl $input.mnc $output.mnc"
    branch.nativemnc = "${branch.nativemnc}.clean_and_center"
}

QC = {
    //Generates QC images
    output.dir="QC"
    exec """
      tmpdir=\$(mktemp -d) &&

      create_verify_image \$tmpdir/t.rgb
        -width 1920 -autocols 6 -autocol_planes t
        -row ${branch.modelmnc}.beastnormalize.mnc color:gray
        volume_overlay:$input.mnc:0.4:red volume_overlay:${REGISTRATIONOUTLINE}:1:blue &&

      create_verify_image \$tmpdir/s.rgb
        -width 1920 -autocols 6 -autocol_planes s
        -row ${branch.modelmnc}.beastnormalize.mnc color:gray
        volume_overlay:$input.mnc:0.4:red volume_overlay:${REGISTRATIONOUTLINE}:1:blue &&

      create_verify_image \$tmpdir/c.rgb
        -width 1920 -autocols 6 -autocol_planes c
        -row ${branch.modelmnc}.beastnormalize.mnc color:gray
        volume_overlay:$input.mnc:0.4:red volume_overlay:${REGISTRATIONOUTLINE}:1:blue &&

        convert -append \$tmpdir/*.rgb $output.jpg &&

        rm -r \$tmpdir
    """

}

preprocess = segment {
	//Default best-practices preprocessing pipeline to run on all data
	convert + n4correct + linear_antsRegistration + cutneckapplyautocrop + beast + QC + resample_to_lsq6_space
}