• Remove this note when ready to publish

• Recommended Page Notes Name: YYYMNMDD <<Note contents>>

• OPTIONAL Includes :
• OPTIONAL Dynamic Tasklist :
  {dynamictasklist:ToDo DESCRIP YYMMDD}
• OPTIONAL Table of Contents: {toc:style=outline|maxLevel=2}

BrainWeb wget script

Those includes multi-modal MR brain images with different levels of noise. That will serve as a good test case of algorithm for different levels of noise.

/IPLlinux/hjohnson/HDNI/brainweb/anatomicalModelOfNormalBrain

• Evaluation of Performance Metrics for Bias Field Correction in MR Brain Images by Chua et al, 2009, Journal of MRI

: Montreal Simulated Brain Data Used for comparing two Bias Correction Methods in the Field at the time.

: Describes three different validation measures along the simulated data

/hjohnson/HDNI/20130225_BRAINSABC_Publication/BRAINSABC_brainweb

[eunyokim@wundt BRAINSABC_brainweb]$ ls
ExtendedAtlasDefinition.xml  Image_1mm_pn0_rf40/          Image_1mm_pn1_rf40/          
Image_1mm_pn3_rf40/          Image_1mm_pn5_rf40/          Image_1mm_pn7_rf40/          
Image_1mm_pn9_rf40/	     Image_1mm_pn0_rf0/           Image_1mm_pn1_rf0/           
Image_1mm_pn3_rf0/           Image_1mm_pn5_rf0/           Image_1mm_pn7_rf0/           
Image_1mm_pn9_rf0/           command.sh
Image_1mm_pn0_rf20/          Image_1mm_pn1_rf20/          Image_1mm_pn3_rf20/          
Image_1mm_pn5_rf20/          Image_1mm_pn7_rf20/          Image_1mm_pn9_rf20/

brainWebDir=/hjohnson/HDNI/brainweb/
for pn in pn0 pn1 pn3 pn5 pn7 pn9
do
    for rf in rf0 rf20 rf40
    do
        inputDir="$brainWebDir/"
        t1="$inputDir/T1_ai_msles2_1mm_${pn}_${rf}.nii"
        t2="$inputDir/T2_ai_msles2_1mm_${pn}_${rf}.nii"
        outputDir="/hjohnson/HDNI/20130225_BRAINSABC_Publication/BRAINSABC_brainweb/Image_1mm_${pn}_${rf}"
        mkdir  $outputDir
        cd $outputDir
        cmd="BRAINSABC  \
            --atlasDefinition /hjohnson/HDNI/20130225_BRAINSABC_Publication/BRAINSABC_brainweb//ExtendedAtlasDefinition.xml \
            --atlasToSubjectTransform $outputDir/atlas_to_subject.h5 \
            --atlasToSubjectTransformType SyN \
            --debuglevel 10 \
            --filterIteration 3 \
            --filterMethod GradientAnisotropicDiffusion \
            --gridSize 10,10,10 \
            --inputVolumeTypes T1,T2 \
            --inputVolumes $t1  \
            --inputVolumes $t2 \
            --interpolationMode Linear \
            --maxBiasDegree 4 \
            --maxIterations 3 \
            --outputDir ./ \
            --outputDirtyLabels $outputDir/volume_label_seg.nii.gz \
            --outputFormat NIFTI \
            --outputLabels brain_label_seg.nii.gz \
            --outputVolumes $outputDir/T1_ai_msles2_1mm_${pn}_${rf}_corrected.nii.gz \
            --outputVolumes $outputDir/T2_ai_msles2_1mm_${pn}_${rf}_corrected.nii.gz \
            --posteriorTemplate POSTERIOR_%s.nii.gz"
        echo "---------------------------------------------------------------------------"
        echo $cmd
        echo "---------------------------------------------------------------------------"
        $cmd |tee $outputDir/command.out
        wait
        cd /hjohnson/HDNI/20130225_BRAINSABC_Publication/BRAINSABC_brainweb/
    done
done

How do we want to use this experiment into paper?

• Is bias computed correctly? --> Looked at computed bias, hard to prove.
• Is label map consistent along different rf and pn? --> It is relatively consistent BUT it did worse when there is NO bias introduced. (If image is perfect, it does not do well)
• Then,, how?

Add description and script here

Interpretation 1:

Graph shows variation of [ mean intensity of WM regions from a) right to left, b) posterior to anterior, or c) inferior to superior ]

• • Each mean intensity normalized by its maximum since we only care about fluctuation.
  • Since each images involve botth Gaussian Noise (pn) and Bias (RF), mean intensity normalized by the mean intensity of image with only Gaussian noise (No Bias).
  • Therefore, I am quite sure that lower standard deviation meaning less bias.
  • In that sense, BRAINSABC with Syn Registration performs best, and Affine registration did a reasonable job.
  • Explanation for POOR performance on certain planes ( Saggital in both T1 and T2, and Coronal in T2): 

Regarding that we only used polynormial approximation for bias field at maximum degree of 4, the simulated bias, which possibly has very higher degree of freedom, cannot be captured entirely. The maximal bias introduced in the direction of axial plane both in T1 and T2, and those are well processed.

To support of idea of enough correction of bias field even with those poor-looking graph, we might need to compare

1. 1. Brain label map from BABC showing it is very robust ( Little or no change across noise/bias level ) and/or
   2. something else???

Well, need to think more carefully.

Any comments/advice will be appreciated.

http://brainweb.bic.mni.mcgill.ca/brainweb/anatomic_normal.html

From the left to right, white matter FUZZY model, label map, and grey matter FUZZY model.

Make sure the fuzzy model is really fuzzy. " The voxel values in these volumes reflects the proportion of tissue present in that voxel, in the range [0,1]."

This attempt will compute statistics, mean and standard deviation, of ROI, white and grey matter from ground truth provided by BSD.

Three experiments are designed based on the paper by Chua et al, 2009:

1. original mask 
2. eroded mask by radius 1
3. dilated mask by radius 1

What the paper by Chua et al suggusted is that, indirect validation could be very unreliable, but the validty of indirect method could be increased by using very conservative (eroded) mask. 

1. Python code for processing:

   Expand

   Code Block ##{ def ErodedMask( maskFilename, radius, outputFilename ): import SimpleITK as sitk mask = sitk.ReadImage( maskFilename ) >0 eroded = sitk.BinaryErode( mask, 1 ) sitk.WriteImage( eroded, outputFilename ) ##} ##{ def DilatedMask( maskFilename, radius, outputFilename ): import SimpleITK as sitk mask = sitk.ReadImage( maskFilename ) >0 dilated = sitk.BinaryDilate( mask, 1 ) sitk.WriteImage( dilated, outputFilename ) ##} ##{ def generateStdCSVFileFromImageAndMask( imageFilename, maskFilename, outputCSVFilename ): import SimpleITK as sitk image = sitk.RescaleIntensity( sitk.ReadImage( imageFilename ), 0, 255 ) mask = sitk.ReadImage( maskFilename ) > 0 #mask = sitk.BinaryThreshold( sitk.ReadImage( mask ), 0.51, 1.0) #vDisplayCenterSlice([image,mask],["t1","mask"], 0.5) dim = image.GetDimension() size = image.GetSize() origin = image.GetOrigin() direction = image.GetDirection() print( "dim:::" + str( dim )) print( "size:::" + str( size )) print( "origin:::" + str( origin )) print( "direction:::" + str( direction ) ) DimName = {0: 'x', 1:'y', 2:'z' } summaryData = [] for iDim in range(0,dim): dimSliceSize = list( size ) dimSliceSize[ iDim ] = 1 for sliceNumber in range(0, size[ iDim ] ): startIndex = [0,0,0] startIndex[ iDim ] = sliceNumber myImageSlice = sitk.RegionOfInterest( image, tuple( dimSliceSize ), tuple( startIndex ) ) myMaskSlice = sitk.RegionOfInterest( mask, tuple( dimSliceSize ), tuple( startIndex ) ) stat = sitk.LabelStatisticsImageFilter() stat.Execute( myImageSlice, myMaskSlice>0 ) labelNo = max( stat.GetValidLabels()) mean = stat.GetMean( labelNo ) std = stat.GetSigma( labelNo ) if stat.GetCount( labelNo ) > 0: summaryData.append( tuple( [DimName[ iDim ], sliceNumber, mean, std] ) ) import csv with open( outputCSVFilename, "w") as the_file: csv.register_dialect("custom", delimiter=",", skipinitialspace=True) writer = csv.writer(the_file, dialect="custom") writer.writerow( tuple(["dimension","sliceNo","mean", "std"]) ) for tup in summaryData: writer.writerow(tup) ##} ## main part is omitted..

2. Usage (help)

   Expand

   Code Block [eunyokim@wundt BRAINSABC_brainweb]$ python BiasCorrectionResultComparisonCMD.py --help usage: BiasCorrectionResultComparisonCMD.py [-h] [--inputErodeFilename INPUTERODEFILENAME] [--inputErodeRadius INPUTERODERADIUS] [--outputErodeFilename OUTPUTERODEDFILENAME] [--inputDilateFilename INPUTDILATEFILENAME] [--inputDilateRadius INPUTDILATERADIUS] [--outputDilateFilename OUTPUTDILATEDFILENAME] [--inputImageFilename INPUTIMAGEFILENAME] [--inputMaskFilename INPUTMASKFILENAME] [--outputCSVFilename OUTPUTCSVFILENAME] [--listOfPosteriors LISTOFPOSTERIORS] [--outputCombinedPosteriorFilename OUTPUTCOMBINEDPOSTERIORFILENAME] **************************** 10-cross validation analysis optional arguments: -h, --help show this help message and exit erodeArg: --inputErodeFilename INPUTERODEFILENAME inputFilename input mask file name to be eroded --inputErodeRadius INPUTERODERADIUS inputRadius input radius to be eroded --outputErodeFilename OUTPUTERODEDFILENAME outputErodedFilename dilateArg: --inputDilateFilename INPUTDILATEFILENAME inputFilename input mask file name to be dilated --inputDilateRadius INPUTDILATERADIUS inputRadius input radius to be dilated --outputDilateFilename OUTPUTDILATEDFILENAME outputDilatedFilename csvFileArg: --inputImageFilename INPUTIMAGEFILENAME inputImageFilename --inputMaskFilename INPUTMASKFILENAME inputMaskFilename --outputCSVFilename OUTPUTCSVFILENAME outputCSVFilename posteriorCombineArg: --listOfPosteriors LISTOFPOSTERIORS List of posteriors to be combined --outputCombinedPosteriorFilename OUTPUTCOMBINEDPOSTERIORFILENAME output combined posterior filename from given in list of posteriors

3. Shell script to compute mean and standard deviation from masks:

   Expand

   Code Block inputDir="/hjohnson/HDNI/20130125_BRAINSABC_Publication/BRAINSABC_brainweb" #phantom_1.0mm_normal_wht.nii #phantom_1.0mm_normal_gry.nii outputDir="/hjohnson/HDNI/20130125_BRAINSABC_Publication/BRAINSABC_brainweb/anatomicalModel" for roi in phantom_1.0mm_normal_wht phantom_1.0mm_normal_gry.nii do refMask="/hjohnson/HDNI/brainweb/anatomicalModelOfNormalBrain/$roi.nii" echo "Process $refMask" echo "====== ====== ====== ====== ====== ====== ====== ====== ====== " ## erode ErodeMask="$outputDir/Erode/${roi}_Erode_R1.nii.gz" mkdir -p "$outputDir/Erode/" python BiasCorrectionResultComparisonCMD.py --inputErodeFilename $refMask --inputErodeRadius 1 --outputErodeFilename $ErodeMask ## dilate DilateMask="$outputDir/Dilate/${roi}_Dilate_R1.nii.gz" mkdir -p "$outputDir/Dilate/" python BiasCorrectionResultComparisonCMD.py --inputDilateFilename $refMask --inputDilateRadius 1 --outputDilateFilename $DilateMask # for erode original and dilated for option in Erode Dilate "" do if [ $option == "" ]; then echo "set mask to original" mask=$refMask currOutputDir="$outputDir/Orig/" fi if [ $option == Erode ]; then echo "set mask to Erodel" mask=$ErodeMask currOutputDir="$outputDir/Erode/" fi if [ $option == Dilate ]; then echo "set mask to Dilate" mask=$DilateMask currOutputDir="$outputDir/Dilate/" fi # processing BABC outputs for currDir in $inputDir/Image* do echo "Processing $currDir..." currDirName=(`basename $currDir`) python BiasCorrectionResultComparisonCMD.py --inputImageFilename ${currDir}/t1_average_BRAINSABC.nii.gz --inputMaskFilename $mask --outputCSVFilename ${currOutputDir}/t1_${currDirName}.csv python BiasCorrectionResultComparisonCMD.py --inputImageFilename ${currDir}/t2_average_BRAINSABC.nii.gz --inputMaskFilename $mask --outputCSVFilename ${currOutputDir}/t2_${currDirName}.csv done # processing raw outputs brainWebDir=/hjohnson/HDNI/brainweb/ for pn in pn0 pn1 pn3 pn5 pn7 pn9 do for rf in rf0 rf20 rf40 do echo "processing raw data of ${pn} ${rf}..." currDir="$brainWebDir/" t1="$currDir/T1_ai_msles2_1mm_${pn}_${rf}.nii" python BiasCorrectionResultComparisonCMD.py --inputImageFilename $t1 --inputMaskFilename $mask --outputCSVFilename $currOutputDir/t1_ai_msles2_1mm_${pn}_${rf}_${roi}_raw.csv t2="$currDir/T2_ai_msles2_1mm_${pn}_${rf}.nii" python BiasCorrectionResultComparisonCMD.py --inputImageFilename $t2 --inputMaskFilename $mask --outputCSVFilename $currOutputDir/t2_ai_msles2_1mm_${pn}_${rf}_${roi}_raw.csv done done done done

4. Above trial looks really ugly. Should think about the validity of computation.

Quick Validation by Sight: Red: Ours, Blue: :Ground Truth from Montreal Group

Our definition includes Globus and Thalamus. As we see, they match very good other than those globus and thalamus area. As it is, it is Dice Index matric is compatabe to Atropos based on what they published.