Skip to content

Commit

Permalink
last anatomicuts changes (freesurfer#841)
Browse files Browse the repository at this point in the history
  • Loading branch information
@vsiless
vsiless authored Jan 7, 2021
1 parent c0a577b commit d0ebb96
Show file tree
Hide file tree
Showing 11 changed files with 782 additions and 326 deletions.
586 changes: 409 additions & 177 deletions anatomicuts/anatomiCutsUtils
View file

Large diffs are not rendered by default.

44 changes: 32 additions & 12 deletions anatomicuts/dmri_ac.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,9 @@ function forAll()
echo ${cluster_call}

cd ${SUBJECTS_DIR}

for s in */;
#for s in `ls -dir [a-zA-Z]*`;
do
subject=${s//[\/]/}
echo ${subject}
Expand Down Expand Up @@ -259,12 +261,12 @@ function preGA()
std=$4
bb=$5

Clean ${subject} ${targetSubject} ${lenght} ${std} ${bb}
Hungarian ${subject} ${targetSubject} ${lenght} ${std} ${bb}
Measures ${subject} ${targetSubject} ${lenght} ${std} ${bb}
ToAnat ${subject} ${targetSubject} ${lenght} ${std} ${bb}
#Clean ${subject} ${targetSubject} ${lenght} ${std} ${bb}
#Hungarian ${subject} ${targetSubject} ${lenght} ${std} ${bb}
#Measures ${subject} ${targetSubject} ${lenght} ${std} ${bb}
#ToAnat ${subject} ${targetSubject} ${lenght} ${std} ${bb}
SurfaceMeasures ${subject} ${targetSubject} ${lenght} ${std} ${bb}
ToTarget ${subject} ${targetSubject} ${lenght} ${std} ${bb}
#ToTarget ${subject} ${targetSubject} ${lenght} ${std} ${bb}

}
function GA()
Expand All @@ -278,8 +280,12 @@ function GA()
groupB=$7
groups=$8
thickness=$9
# mkdir -p ${ODMRI_DIR}/average/dmri.ac/${2}/${3}/
# rm ${ODMRI_DIR}/average/dmri.ac/${2}/${3}/*
#${FREESURFER_HOME}/bin/anatomiCutsUtils -f GAL -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups}
#${FREESURFER_HOME}/bin/anatomiCutsUtils -f GA -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups}
${FREESURFER_HOME}/bin/anatomiCutsUtils -f thicknessPerStructure -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} -t ${thickness}
${FREESURFER_HOME}/bin/anatomiCutsUtils -f diffAlongTheLine -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups}
#${FREESURFER_HOME}/bin/anatomiCutsUtils -f thicknessPerStructure -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} -t ${thickness}
#${FREESURFER_HOME}/bin/anatomiCutsUtils -f connectivityGraph -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} -t ${thickness}

#${FREESURFER_HOME}/bin/anatomiCutsUtils -f GA -m "DKI" -cf "/space/snoke/1/public/vivros/data/demos_fullID.csv" -cc 0:6 -cta 200 -ts ${targetSubject} -s ${ODMRI_DIR} -d " " -ga 3 -gb 1 -l ${lenght} -std ${std}
Expand Down Expand Up @@ -323,6 +329,7 @@ function Hungarian()
ci=${ODMRI_DIR}/${targetSubject}/dmri.ac/${lenght}/${std}
cj=${ODMRI_DIR}/${subject}/dmri.ac/${lenght}/${std}
mkdir -p ${cj}/match/
echo ${HungarianBin} -s1 ${si} -s2 ${sj} -h1 ${ci}/ -h2 ${cj}/ -o ${cj}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -labels -hungarian -c ${c} ${bb}
${HungarianBin} -s1 ${si} -s2 ${sj} -h1 ${ci}/ -h2 ${cj}/ -o ${cj}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -labels -hungarian -c ${c} ${bb}
done
}
Expand All @@ -341,11 +348,12 @@ function Measures()

if [[ -e ${SUBJECTS_DIR}/${subject}/dmri/DKI/dki_RD.nii.gz ]] ;
then
string="${stats_ac_bin} -i ${anatomicuts}/ -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 7 FA ${diff}/DKI/dki_FA.nii MD ${diff}/DKI/dki_MD.nii RD ${diff}/DKI/dki_RD.nii AD ${diff}/DKI/dki_AD.nii MK ${diff}/DKI/dki_MK.nii RK ${diff}/DKI/dki_RK.nii AK ${diff}/DKI/dki_AK.nii -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
string="${stats_ac_bin} -i ${anatomicuts}/ -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 7 FA ${diff}/DKI/dki_FA.nii.gz MD ${diff}/DKI/dki_MD.nii.gz RD ${diff}/DKI/dki_RD.nii.gz AD ${diff}/DKI/dki_AD.nii.gz MK ${diff}/DKI/dki_MK.nii.gz RK ${diff}/DKI/dki_RK.nii.gz AK ${diff}/DKI/dki_AK.nii.gz -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
else
string="${stats_ac_bin} -i ${anatomicuts}/ -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 4 FA ${diff}/DTI/dti_FA.nii.gz MD ${diff}/DTI/dti_MD.nii RD ${diff}/DTI/dti_RD.nii AD ${diff}/DTI/dti_AD.nii -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
string="${stats_ac_bin} -i ${anatomicuts}/ -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 4 FA ${diff}/DTI/dti_FA.nii.gz.gz MD ${diff}/DTI/dti_MD.nii.gz RD ${diff}/DTI/dti_RD.nii.gz AD ${diff}/DTI/dti_AD.nii.gz -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
fi

echo ${string}
${string}
done
}
Expand All @@ -363,8 +371,9 @@ function SurfaceMeasures()
annot=${SUBJECTS_DIR}/${subject}/label/
mri=${SUBJECTS_DIR}/${subject}/mri/

diff=${ODMRI_DIR}/${subject}/dmri/

${FREESURFER_HOME}/bin/dmri_extractSurfaceMeasurements -i ${anatomicuts}/*trk -sl ${surf}/lh.pial -tl ${surf}/lh.thickness -cl ${surf}/lh.curv.pial -sr ${surf}/rh.pial -tr ${surf}/rh.thickness -cr ${surf}/rh.curv.pial -rid ${mri}/brain.nii.gz -ria ${mri}/brain.nii.gz -al ${annot}/lh.aparc.annot -ar ${annot}/rh.aparc.annot -o ${anatomicutsdiff}/measures/ -p ${anatomicutsdiff}/match/${targetSubject}_${subject}_c200_hungarian.csv
${FREESURFER_HOME}/bin/dmri_extractSurfaceMeasurements -i ${anatomicuts}/*.trk -sl ${surf}/lh.pial -tl ${surf}/lh.thickness -cl ${surf}/lh.curv.pial -sr ${surf}/rh.pial -tr ${surf}/rh.thickness -cr ${surf}/rh.curv.pial -rid ${mri}/brain.nii.gz -ria ${mri}/brain.nii.gz -al ${annot}/lh.aparc.annot -ar ${annot}/rh.aparc.annot -o ${anatomicutsdiff}/measures/ -p ${anatomicutsdiff}/match/${targetSubject}_${subject}_c200_hungarian.csv -fa 7 FA ${diff}/DKI/dki_FA.nii.gz MD ${diff}/DKI/dki_MD.nii.gz RD ${diff}/DKI/dki_RD.nii.gz AD ${diff}/DKI/dki_AD.nii.gz MK ${diff}/DKI/dki_MK.nii.gz RK ${diff}/DKI/dki_RK.nii.gz AK ${diff}/DKI/dki_AK.nii.gz

}
function ToAnat()
Expand Down Expand Up @@ -392,6 +401,7 @@ function ToAnat()
mkdir -p ${ODMRI_DIR}/${subject}/dmri.ac/${lenght}/${std}/${common}/

common_clustering=${ODMRI_DIR}/${subject}/dmri.ac/${lenght}/${std}/${common}/
mkdir ${diff}/xfms/

flirt -in ${diff}/DTI/dti_FA.nii.gz -ref ${mri}/brain.nii.gz -omat ${diff}/xfms/fa2brain.mat

Expand Down Expand Up @@ -464,7 +474,13 @@ function average()
targetSubject=$1
lenght=$2
std=$3

labels_file=$4
labels_cols=$5
groupA=$6
groupB=$7
groups=$8
thickness=$9
#
mkdir -p ${ODMRI_DIR}/average/dmri.ac/${lenght}/${std}/images
#correspondences="["
#imagesFolder="["
Expand All @@ -480,23 +496,27 @@ function average()
if [ ${#correspondences} -ge 3 ]; then
correspondences=${correspondences},
imagesFolder=${imagesFolder},
subjectsNames=${subjectsNames},
fi
correspondences=${correspondences}${ODMRI_DIR}/${s}/dmri.ac/${lenght}/${std}/match/${s2}_${s}_c200_hungarian.csv
imagesFolder=${imagesFolder}${ODMRI_DIR}/${s}/dmri.ac/${lenght}/${std}/to${targetSubject}/images/
subjectsNames=${subjectsNames}${s}

fi
done
correspondences=${correspondences}
imagesFolder=${imagesFolder}
subjectsNames=${subjectsNames}

echo $correspondences
echo $imagesFolder
clusterIndeces=[$( echo `seq 0 1 199` | sed 's/ /,/g' )]
echo $clusterIndeces
mkdir -p ${outputFolder}
#correspondences, imagesFolder, outputFolder, clusterIndeces
cd /space/erebus/2/users/vsiless/code/freesurfer/anatomicuts/
cd ${SUBJECTS_DIR}
#pbsubmit -n 1 -c "python3 -c \"import anatomiCutsUtils; anatomiCutsUtils.averageCorrespondingClusters($correspondences, $imagesFolder, $outputFolder,$clusterIndeces) \" "
anatomiCutsUtils -f averageCorrespondingClusters -co ${correspondences} -if ${imagesFolder} -of ${outputFolder} -in ${clusterIndeces}
${FREESURFER_HOME}/bin/anatomiCutsUtils -f averageCorrespondingClusters -co ${correspondences} -if ${imagesFolder} -of ${outputFolder} -in ${clusterIndeces} -sn ${subjectsNames} -cf "${labels_file}" -cc ${labels_cols} -d "," -pt ${groups}
#python3 -c "import anatomiCutsUtils; anatomiCutsUtils.averageCorrespondingClusters(${correspondences}, ${imagesFolder}, ${outputFolder},${clusterIndeces}) "

}
Expand Down
156 changes: 105 additions & 51 deletions anatomicuts/dmri_extractSurfaceMeasurements.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,7 +226,7 @@ int main(int narg, char* arg[])

ClusterToolsType::Pointer clusterTools = ClusterToolsType::New();
clusterTools->GetPolyDatas(TRKFiles, &polydatas, ref_Image.at(0));
meshes = clusterTools->PolydataToMesh(polydatas);
meshes = clusterTools->FixSampleClusters(polydatas,10);

//Loading the surface for each hemisphere and metric
//Left Curvature
Expand Down Expand Up @@ -322,8 +322,8 @@ int main(int narg, char* arg[])
getline ( file, value, ',' );
long long v2 = atoll(value.c_str());
correspondences.push_back(v2);
std::cout << " v2 " << v2 << std::endl;
}
std::cout << " Files " << meshes->size() << std::endl;
// Cycling through the TRK files
for(int i = 0; i < meshes->size(); i++)
{
Expand All @@ -343,9 +343,12 @@ int main(int narg, char* arg[])
if (FA_FOUND)
{
for (int a = 0; a < image_fileNames.size(); a++)
oFile << ", mean" << image_fileNames.at(a) << ", stde" << image_fileNames.at(a);
oFile << ",mean" << image_fileNames.at(a) << ",stde" << image_fileNames.at(a);
for (int a = 0; a < image_fileNames.size(); a++)
for (int b =0; b<10; b++)
oFile << ","<< b << "_"<< image_fileNames.at(a);
}
oFile << endl;
oFile << ","<<endl;

// Initialization of a new stream for every TRK files

Expand All @@ -357,22 +360,25 @@ int main(int narg, char* arg[])
for (; inputCellIt != input->GetCells()->End(); ++inputCellIt, ++counter)
{
vector<float> meanFA;
vector<vector<float>> allFA;
vector<float> stdeFA;

// If there are image files, then find the mean and stde of FA
if (FA_FOUND)
{
{

for (int p = 0; p < volumes.size(); p++)
{
allFA.push_back(vector<float>());
// Creating streamline variable and finding first point
CellType::PointIdIterator it = inputCellIt.Value()->PointIdsBegin();
input->GetPoint(*it, &firstPt);

input->GetPoint(*inputCellIt.Value()->PointIdsEnd(), &lastPt);
float val = (firstPt[0] -lastPt[0])*(firstPt[1]-lastPt[1])*(firstPt[2]-lastPt[2]);

// Getting the FA value at all points
vector<float> FA_values;
ImageType::IndexType index;
if (volumes.at(p)->TransformPhysicalPointToIndex(firstPt, index))
FA_values.push_back(volumes.at(p)->GetPixel(index));

// Cycling through the points in the stream
for (; it != inputCellIt.Value()->PointIdsEnd(); it++)
Expand All @@ -381,7 +387,13 @@ int main(int narg, char* arg[])

// If FA options is used, then add the FA values to the vector
if (volumes.at(p)->TransformPhysicalPointToIndex(lastPt, index))
{
FA_values.push_back(volumes.at(p)->GetPixel(index));
if( val>0)
allFA[p].push_back(volumes.at(p)->GetPixel(index));
else
allFA[p].insert(allFA[p].begin(),volumes.at(p)->GetPixel(index));
}
}

// Calculating the MeanFA and stdeFA
Expand Down Expand Up @@ -437,81 +449,123 @@ int main(int narg, char* arg[])


// Outputting values to the file
std::cout << ID1 << " " <<ID2<< std::endl;
oFile << "StreamLine" << counter << ", ";
int structure;
if (ID1 == Left_ID1)
if (ID1 >=0)
{
CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID1].annotation, &structure);
std::cout << ID1 << " " << surfCL->vertices[ID1].annotation << " " <<structure<< std::endl;
if (ID1 == Left_ID1)
{
CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID1].annotation, &structure);
std::cout <<"L" << ID1 << " " << surfCL->vertices[ID1].annotation << " " <<structure<< std::endl;
}
else
{
CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID1].annotation, &structure);
std::cout << "R" << ID1 << " " << surfCR->vertices[ID1].annotation << " " <<structure<< std::endl;
}
oFile << structure << ",";
}
else
{
CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID1].annotation, &structure);
std::cout << ID1 << " " << surfCR->vertices[ID1].annotation << " " <<structure<< std::endl;
oFile << ",";
}
oFile << structure << ",";
if (ID2 == Left_ID2)
if(ID2>=0)
{
CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID2].annotation, &structure);
std::cout << ID2 << " " << surfCL->vertices[ID2].annotation << " " <<structure<< std::endl;
if (ID2 == Left_ID2)
{
CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID2].annotation, &structure);
std::cout << "L"<< ID2 << " " << surfCL->vertices[ID2].annotation << " " <<structure<< std::endl;
}
else
{
CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID2].annotation, &structure);
std::cout <<"R"<< ID2 << " " << surfCR->vertices[ID2].annotation << " " <<structure<< std::endl;
}

oFile << structure << ",";
}
else
{
CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID2].annotation, &structure);
std::cout << ID2 << " " << surfCR->vertices[ID2].annotation << " " <<structure<< std::endl;
oFile << ",";
}

oFile << structure << ",";

if (ID1 == Left_ID1)
if (ID1 >=0)
{
oFile << surfCL->vertices[ID1].curv << ",";
values[0]+= surfCL->vertices[ID1].curv ;
if (ID1 == Left_ID1)
{
oFile << surfCL->vertices[ID1].curv << ",";
values[0]+= surfCL->vertices[ID1].curv ;
}
else
{
oFile << surfCR->vertices[ID1].curv << ",";
values[0]+= surfCR->vertices[ID1].curv ;
}
}
else
{
oFile << surfCR->vertices[ID1].curv << ",";
values[0]+= surfCR->vertices[ID1].curv ;
oFile << ",";
}

if (ID2 == Left_ID2)
if (ID2 >=0)
{
oFile << surfCL->vertices[ID2].curv << ",";
values[1]+= surfCL->vertices[ID2].curv ;

if (ID2 == Left_ID2)
{
oFile << surfCL->vertices[ID2].curv << ",";
values[1]+= surfCL->vertices[ID2].curv ;

}
else
{
oFile << surfCR->vertices[ID2].curv << ",";
values[1]+= surfCR->vertices[ID2].curv ;
}
}
else
{
oFile << surfCR->vertices[ID2].curv << ",";
values[1]+= surfCR->vertices[ID2].curv ;
}

if (ID1 == Left_ID1)
oFile << ",";
}

if (ID1 >=0)
{
oFile << surfTL->vertices[ID1].curv << ",";
values[2]+= surfTL->vertices[ID1].curv ;
if (ID1 == Left_ID1)
{
oFile << surfTL->vertices[ID1].curv << ",";
values[2]+= surfTL->vertices[ID1].curv ;
}
else
{
oFile << surfTR->vertices[ID1].curv << ",";
values[2]+= surfTR->vertices[ID1].curv ;
}
}
else
{
oFile << surfTR->vertices[ID1].curv << ",";
values[2]+= surfTR->vertices[ID1].curv ;
oFile << ",";
}


if (ID2 == Left_ID2)
if (ID2 >=0)
{
oFile << surfTL->vertices[ID2].curv;
values[3]+= surfTL->vertices[ID2].curv ;
if (ID2 == Left_ID2)
{
oFile << surfTL->vertices[ID2].curv <<",";
values[3]+= surfTL->vertices[ID2].curv ;
}
else
{
oFile << surfTR->vertices[ID2].curv <<",";
values[3]+= surfTR->vertices[ID2].curv ;
}
}
else
{
oFile << surfTR->vertices[ID2].curv;
values[3]+= surfTR->vertices[ID2].curv ;
}
oFile << ",";
}
if (FA_FOUND)
{
for (int m = 0; m < stdeFA.size(); m++)
oFile << "," << meanFA.at(m) << "," << stdeFA.at(m);
for (int m = 0; m < volumes.size(); m++)
oFile << meanFA.at(m) << "," << stdeFA.at(m)<<",";
for (int m = 0; m < volumes.size(); m++)
for (int b = 0; b < allFA[m].size(); b++)
oFile << allFA[m][b]<<",";
}

oFile << endl;
Expand Down
2 changes: 1 addition & 1 deletion anatomicuts/dmri_neighboringRegions.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,6 @@ using namespace std;
int main(int argc, char *argv[])
{
// TODO

std::cout << " chau " << std::endl;
return 0;
}
5 changes: 5 additions & 0 deletions freeview/LayerSurface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -809,7 +809,12 @@ void LayerSurface::InitializeActors()

vtkSmartPointer<vtkPolyDataMapper> mapper3 = vtkSmartPointer<vtkPolyDataMapper>::New();
m_vertexPoly2D[i] = vtkSmartPointer<vtkPolyData>::New();
#if VTK_MAJOR_VERSION > 5
mapper3->SetInputData(m_vertexPoly2D[i]);
#else
mapper3->SetInput(m_vertexPoly2D[i]);
#endif

mapper3->ScalarVisibilityOff();
m_vertexActor2D[i]->SetMapper(mapper3);
m_vertexActor2D[i]->SetProperty( m_vertexActor2D[i]->MakeProperty() );
Expand Down
Loading

0 comments on commit d0ebb96

Please sign in to comment.