From f37a28523298ee54b6ab4c13d3cff00eb9bc1d0f Mon Sep 17 00:00:00 2001 From: Gregoire Ville Date: Thu, 6 Jun 2024 17:44:51 +0200 Subject: [PATCH 1/8] add of a feature to convert a tracks file to a density map --- Anima/math-tools/common_tools/CMakeLists.txt | 1 + .../shapes_to_density/CMakeLists.txt | 38 ++++++ .../animaShapesToDensity.cxx | 111 ++++++++++++++++++ 3 files changed, 150 insertions(+) create mode 100644 Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt create mode 100644 Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx diff --git a/Anima/math-tools/common_tools/CMakeLists.txt b/Anima/math-tools/common_tools/CMakeLists.txt index 4a0a3b040..5b64b5c1b 100644 --- a/Anima/math-tools/common_tools/CMakeLists.txt +++ b/Anima/math-tools/common_tools/CMakeLists.txt @@ -6,5 +6,6 @@ add_subdirectory(collapse_image) add_subdirectory(convert_image) add_subdirectory(convert_shape) add_subdirectory(merge_block_images) +add_subdirectory(shapes_to_density) add_subdirectory(vectorize_images) diff --git a/Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt b/Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt new file mode 100644 index 000000000..f0dc72832 --- /dev/null +++ b/Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt @@ -0,0 +1,38 @@ +if(BUILD_TOOLS) + +project(animaShapesToDensity) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +## ############################################################################# +## add executable +## ############################################################################# + +add_executable(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + ${VTK_LIBRARIES} + ${ITKIO_LIBRARIES} + AnimaDataIO + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx new file mode 100644 index 000000000..e93662993 --- /dev/null +++ b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx @@ -0,0 +1,111 @@ +#include +#include +#include +#include + +#include + +#include +#include + +template int sgn(T val) { + return (T(0) < val) - (val < T(0)); +} + +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + TCLAP::ValueArg inTrackArg("i","in-tracks","input track file (vtp, vtk)",true,"","input tracks",cmd); + TCLAP::ValueArg geometryArg("g","geometry","Output image geometry",true,"","output geometry",cmd); + TCLAP::ValueArg outArg("o","out","Output image representing fibers convert in binary image",true,"","output image",cmd); + + try + { + cmd.parse(argc,argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + typedef itk::Image ImageType; + + ImageType::Pointer geomImage = anima::readImage (geometryArg.getValue()); + + using OutImageType = itk::Image ; + +// OutImageType::DirectionType direction; +// direction.SetIdentity(); +// OutImageType::RegionType region; + OutImageType::Pointer outImage = OutImageType::New(); + + outImage->Initialize(); + outImage->SetRegions(geomImage->GetLargestPossibleRegion()); + outImage->SetSpacing(geomImage->GetSpacing()); + outImage->SetOrigin(geomImage->GetOrigin()); + outImage->SetDirection(geomImage->GetDirection()); + + outImage->Allocate(); + outImage->FillBuffer(0); + + anima::ShapesReader trackReader; + trackReader.SetFileName(inTrackArg.getValue()); + trackReader.Update(); + + vtkSmartPointer tracks = trackReader.GetOutput(); + + vtkIdType nbCells = tracks->GetNumberOfCells(); + + double ptVals[3], ptValsNext[3], ptValsTmp[3]; + OutImageType::IndexType index; + OutImageType::PointType point; + + for (int j = 0;j < nbCells;++j) + { + vtkCell *cell = tracks->GetCell(j); + vtkPoints *cellPts = cell->GetPoints(); + vtkIdType nbPts = cellPts->GetNumberOfPoints(); + + for (int i = 0;i < nbPts-1 ;++i) + { + cellPts->GetPoint(i,ptVals); + cellPts->GetPoint(i+1,ptValsNext); + + if(i == 0) + { + for (unsigned int k = 0;k < 3;++k) + point[k] = ptVals[k]; + + outImage->TransformPhysicalPointToIndex(point,index); + outImage->SetPixel(index, 1); + } + + double dist = 2; + double dx, dy, dz; + while(dist >= 1.5) + { + dx = ptValsNext[0] - ptVals[0]; + dy = ptValsNext[1] - ptVals[1]; + dz = ptValsNext[2] - ptVals[2]; + + dist = std::sqrt(dx*dx + dy*dy + dz*dz); + + ptVals[0] = ptVals[0] + sgn(dx); + ptVals[1] = ptVals[1] + sgn(dy); + ptVals[2] = ptVals[2] + sgn(dz); + + for (unsigned int k = 0;k < 3;++k) + point[k] = ptVals[k]; + + outImage->TransformPhysicalPointToIndex(point,index); + outImage->SetPixel(index, outImage->GetPixel(index) + 1.0); + } + } + } + + anima::writeImage (outArg.getValue(), outImage); + + return EXIT_SUCCESS; +} From 5d474be19eeb5f5946c8cb3cee08385b0b4d7720 Mon Sep 17 00:00:00 2001 From: Gregoire Ville Date: Wed, 10 Jul 2024 18:20:56 +0200 Subject: [PATCH 2/8] Improve the main algorithm and comment the code --- .../animaShapesToDensity.cxx | 108 +++++++++++------- 1 file changed, 65 insertions(+), 43 deletions(-) diff --git a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx index e93662993..976c68379 100644 --- a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx +++ b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx @@ -8,18 +8,21 @@ #include #include -template int sgn(T val) { - return (T(0) < val) - (val < T(0)); -} - +/** + * @brief Converts an input track file (vtp or vtk) into a density map. + * The density map is a 3D image giving for each voxel the number of streamlines from the input image passing through it. + * + */ int main(int argc, char **argv) { TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + //Define arguments TCLAP::ValueArg inTrackArg("i","in-tracks","input track file (vtp, vtk)",true,"","input tracks",cmd); TCLAP::ValueArg geometryArg("g","geometry","Output image geometry",true,"","output geometry",cmd); TCLAP::ValueArg outArg("o","out","Output image representing fibers convert in binary image",true,"","output image",cmd); + //Try to parse try { cmd.parse(argc,argv); @@ -30,82 +33,101 @@ int main(int argc, char **argv) return EXIT_FAILURE; } + //Define types typedef itk::Image ImageType; - - ImageType::Pointer geomImage = anima::readImage (geometryArg.getValue()); - using OutImageType = itk::Image ; -// OutImageType::DirectionType direction; -// direction.SetIdentity(); -// OutImageType::RegionType region; + //Read geometry image and allocate memory for output image + ImageType::Pointer geomImage = anima::readImage (geometryArg.getValue()); OutImageType::Pointer outImage = OutImageType::New(); + //Give to the output image the same properties as the geometry image outImage->Initialize(); outImage->SetRegions(geomImage->GetLargestPossibleRegion()); outImage->SetSpacing(geomImage->GetSpacing()); outImage->SetOrigin(geomImage->GetOrigin()); outImage->SetDirection(geomImage->GetDirection()); + //Allocate the correct size for the output image and intialize it with 0 in each voxel outImage->Allocate(); - outImage->FillBuffer(0); + outImage->FillBuffer(0.0); + //Read input track file anima::ShapesReader trackReader; trackReader.SetFileName(inTrackArg.getValue()); trackReader.Update(); + //Initialize a vtksmartpointer with the content of the track file, and store the total number of streamlines (1 cell = 1 streamline) vtkSmartPointer tracks = trackReader.GetOutput(); - vtkIdType nbCells = tracks->GetNumberOfCells(); - double ptVals[3], ptValsNext[3], ptValsTmp[3]; - OutImageType::IndexType index; - OutImageType::PointType point; + //Initialize data objects for the main algorithm + double ptVals[3], ptValsNext[3]; + OutImageType::IndexType index, indexNext; + OutImageType::PointType point, pointNext; + - for (int j = 0;j < nbCells;++j) + //Main algorithm + for (int j = 0; j < nbCells; ++j) { + //Get streamline j, its points and its number of points vtkCell *cell = tracks->GetCell(j); vtkPoints *cellPts = cell->GetPoints(); vtkIdType nbPts = cellPts->GetNumberOfPoints(); - for (int i = 0;i < nbPts-1 ;++i) + //Iterate through the streamline j's points + for (int i = 0; i < nbPts-1; ++i) { - cellPts->GetPoint(i,ptVals); - cellPts->GetPoint(i+1,ptValsNext); - - if(i == 0) + //Get the current point and the following one, write their coordinates in point and pointNext + cellPts->GetPoint(i, ptVals); + cellPts->GetPoint(i+1, ptValsNext); + for (unsigned int k = 0; k < 3; ++k) { - for (unsigned int k = 0;k < 3;++k) point[k] = ptVals[k]; - - outImage->TransformPhysicalPointToIndex(point,index); - outImage->SetPixel(index, 1); + pointNext[k] = ptValsNext[k]; } - double dist = 2; - double dx, dy, dz; - while(dist >= 1.5) + //Infer the points' index in the output image (namely voxel coordinates) from the points' real coordinates. Write the results in index and indexNext. + outImage->TransformPhysicalPointToIndex(point, index); + outImage->TransformPhysicalPointToIndex(pointNext, indexNext); + + //if index==indexNext, the current point and the following one are in the same voxel + //so we don't increment the number of streamlines passing through this voxel for now: we will do it in the next iteration + //that allows not to count one streamline several times in a voxel + if (index != indexNext) { - dx = ptValsNext[0] - ptVals[0]; - dy = ptValsNext[1] - ptVals[1]; - dz = ptValsNext[2] - ptVals[2]; - - dist = std::sqrt(dx*dx + dy*dy + dz*dz); - - ptVals[0] = ptVals[0] + sgn(dx); - ptVals[1] = ptVals[1] + sgn(dy); - ptVals[2] = ptVals[2] + sgn(dz); - - for (unsigned int k = 0;k < 3;++k) - point[k] = ptVals[k]; - - outImage->TransformPhysicalPointToIndex(point,index); + //increment the number of streamlines passing through the voxel with this index outImage->SetPixel(index, outImage->GetPixel(index) + 1.0); + + double dx = indexNext[0] - index[0]; //difference in x coordinate + double dy = indexNext[1] - index[1]; //difference in y coordinate + double dz = indexNext[2] - index[2]; //difference in z coordinate + double maxStep = std::max(std::max(std::abs(dx), std::abs(dy)), std::abs(dz)); //maximum between the three absolute differences + for (int l = 1; l < maxStep; l++) + { + //We enter this loop only if maxStep >=2. + //In that case, we create artificial index points, to try to modify all and only all the voxels through which the streamline passes between point and pointNext + //We repeat maxStep - 1 times and not maxStep times, because otherwise, index = indexNext, and we would increment indexNext twice + index[0] = index[0] + std::round(dx/maxStep); + index[1] = index[1] + std::round(dy/maxStep); + index[2] = index[2] + std::round(dz/maxStep); + outImage->SetPixel(index, outImage->GetPixel(index) + 1.0); //increment the number of streamlines passing through the voxel with this index + } } } + + //before moving to the next streamline, it remains to consider the last point of this streamline (its voxel has not been incremented yet) + cellPts->GetPoint(nbPts-1, ptVals); + for (unsigned int k = 0; k < 3; ++k) + { + point[k] = ptVals[k]; + } + outImage->TransformPhysicalPointToIndex(point, index); + outImage->SetPixel(index, outImage->GetPixel(index) + 1.0); } - anima::writeImage (outArg.getValue(), outImage); + //Write output image and exit + anima::writeImage (outArg.getValue(), outImage); return EXIT_SUCCESS; } From 86273fa6bc14c649ab5a493ecbb01d2362b5a84b Mon Sep 17 00:00:00 2001 From: Florent Date: Fri, 12 Jul 2024 14:39:50 +0200 Subject: [PATCH 3/8] Update animaShapesToDensity.cxx Improve doxydoc --- .../common_tools/shapes_to_density/animaShapesToDensity.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx index 976c68379..55a5b29e6 100644 --- a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx +++ b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx @@ -10,7 +10,7 @@ /** * @brief Converts an input track file (vtp or vtk) into a density map. - * The density map is a 3D image giving for each voxel the number of streamlines from the input image passing through it. + * @details The density map is a 3D image giving for each voxel the number of streamlines from the input image passing through it. * */ int main(int argc, char **argv) From 174018442b979b1f463dc562c88c406ea0f4a31f Mon Sep 17 00:00:00 2001 From: Gregoire Ville Date: Fri, 12 Jul 2024 18:10:18 +0200 Subject: [PATCH 4/8] Correction of an error in the main algorithm part --- .../shapes_to_density/animaShapesToDensity.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx index 55a5b29e6..899732895 100644 --- a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx +++ b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx @@ -63,7 +63,7 @@ int main(int argc, char **argv) //Initialize data objects for the main algorithm double ptVals[3], ptValsNext[3]; - OutImageType::IndexType index, indexNext; + OutImageType::IndexType index, indexNext, indexTemp; OutImageType::PointType point, pointNext; @@ -108,10 +108,10 @@ int main(int argc, char **argv) //We enter this loop only if maxStep >=2. //In that case, we create artificial index points, to try to modify all and only all the voxels through which the streamline passes between point and pointNext //We repeat maxStep - 1 times and not maxStep times, because otherwise, index = indexNext, and we would increment indexNext twice - index[0] = index[0] + std::round(dx/maxStep); - index[1] = index[1] + std::round(dy/maxStep); - index[2] = index[2] + std::round(dz/maxStep); - outImage->SetPixel(index, outImage->GetPixel(index) + 1.0); //increment the number of streamlines passing through the voxel with this index + indexTemp[0] = index[0] + std::round(l*dx/maxStep); + indexTemp[1] = index[1] + std::round(l*dy/maxStep); + indexTemp[2] = index[2] + std::round(l*dz/maxStep); + outImage->SetPixel(indexTemp, outImage->GetPixel(indexTemp) + 1.0); //increment the number of streamlines passing through the voxel with the index indexTemp } } } From 8cb68d0c300cb9ac6a44fc033325c2b56d7cbce3 Mon Sep 17 00:00:00 2001 From: Gregoire Ville Date: Tue, 29 Oct 2024 16:03:35 +0100 Subject: [PATCH 5/8] Fix syntax typos --- .../shapes_to_density/animaShapesToDensity.cxx | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx index 899732895..9604a88c4 100644 --- a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx +++ b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx @@ -20,7 +20,7 @@ int main(int argc, char **argv) //Define arguments TCLAP::ValueArg inTrackArg("i","in-tracks","input track file (vtp, vtk)",true,"","input tracks",cmd); TCLAP::ValueArg geometryArg("g","geometry","Output image geometry",true,"","output geometry",cmd); - TCLAP::ValueArg outArg("o","out","Output image representing fibers convert in binary image",true,"","output image",cmd); + TCLAP::ValueArg outArg("o","out","Output density map",true,"","output image",cmd); //Try to parse try @@ -106,12 +106,14 @@ int main(int argc, char **argv) for (int l = 1; l < maxStep; l++) { //We enter this loop only if maxStep >=2. - //In that case, we create artificial index points, to try to modify all and only all the voxels through which the streamline passes between point and pointNext + //In that case, we create artificial index points, + // to try to modify all and only all the voxels through which the streamline passes between point and pointNext //We repeat maxStep - 1 times and not maxStep times, because otherwise, index = indexNext, and we would increment indexNext twice indexTemp[0] = index[0] + std::round(l*dx/maxStep); indexTemp[1] = index[1] + std::round(l*dy/maxStep); indexTemp[2] = index[2] + std::round(l*dz/maxStep); - outImage->SetPixel(indexTemp, outImage->GetPixel(indexTemp) + 1.0); //increment the number of streamlines passing through the voxel with the index indexTemp + //increment the number of streamlines passing through the voxel with the index indexTemp + outImage->SetPixel(indexTemp, outImage->GetPixel(indexTemp) + 1.0); } } } From 53f7495ab28b1c1c939d20cd5a865cd25a60d136 Mon Sep 17 00:00:00 2001 From: Gregoire Ville Date: Wed, 20 Nov 2024 15:15:26 +0100 Subject: [PATCH 6/8] Merge actual animaFibersCounter.cxx and animaShapesToDensity.cxx inside a file animaFibersCounter_fusion.cxx + delete previous animaFibersCounter.cxx and shapes_to_density folder --- .../fibers_counter/animaFibersCounter.cxx | 103 -------- .../animaFibersCounter_fusion.cxx | 237 ++++++++++++++++++ Anima/math-tools/common_tools/CMakeLists.txt | 1 - .../shapes_to_density/CMakeLists.txt | 38 --- .../animaShapesToDensity.cxx | 135 ---------- 5 files changed, 237 insertions(+), 277 deletions(-) delete mode 100644 Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx create mode 100644 Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx delete mode 100644 Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt delete mode 100644 Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx diff --git a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx deleted file mode 100644 index 63a8a2a25..000000000 --- a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx +++ /dev/null @@ -1,103 +0,0 @@ -#include - -#include -#include - -#include -#include -#include -#include -#include - -#include - -int main(int argc, char **argv) -{ - TCLAP::CmdLine cmd("Filters fibers from a vtp file using a label image and specifying with several -t and -f which labels should be touched or are forbidden for each fiber. INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION); - - TCLAP::ValueArg inArg("i","input","input tracks file",true,"","input tracks",cmd); - TCLAP::ValueArg outArg("o","output","output mask image",true,"","output mask image",cmd); - TCLAP::ValueArg geomArg("g","geometry","Geometry image",true,"","geometry image",cmd); - - TCLAP::SwitchArg proportionArg("P","proportion","Output proportion of fibers going through each pixel",cmd,false); - - try - { - cmd.parse(argc,argv); - } - catch (TCLAP::ArgException& e) - { - std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; - return EXIT_FAILURE; - } - - typedef itk::Image OutputImageType; - OutputImageType::Pointer outputImage = anima::readImage (geomArg.getValue()); - outputImage->FillBuffer(0.0); - - anima::ShapesReader trackReader; - trackReader.SetFileName(inArg.getValue()); - trackReader.Update(); - - vtkSmartPointer tracks = trackReader.GetOutput(); - - vtkIdType nbCells = tracks->GetNumberOfCells(); - double incrementFactor = 1.0; - if (proportionArg.isSet()) - incrementFactor /= nbCells; - - double ptVals[3]; - OutputImageType::IndexType currentIndex; - OutputImageType::PointType currentPoint; - OutputImageType::RegionType region = outputImage->GetLargestPossibleRegion(); - - // Explores individual fibers - for (int j = 0;j < nbCells;++j) - { - vtkCell *cell = tracks->GetCell(j); - vtkPoints *cellPts = cell->GetPoints(); - vtkIdType nbPts = cellPts->GetNumberOfPoints(); - - // Explores points in fibers - for (int i = 0;i < nbPts;++i) - { - cellPts->GetPoint(i,ptVals); - - for (unsigned int k = 0;k < 3;++k) - currentPoint[k] = ptVals[k]; - - outputImage->TransformPhysicalPointToIndex(currentPoint,currentIndex); - bool insideBuffer = true; - for (unsigned int k = 0;k < 3;++k) - { - if ((currentIndex[k] < 0)||(currentIndex[k] >= region.GetSize(k))) - { - insideBuffer = false; - break; - } - } - - if (!insideBuffer) - continue; - - double countIndex = outputImage->GetPixel(currentIndex) + incrementFactor; - outputImage->SetPixel(currentIndex, countIndex); - } - } - - if (proportionArg.isSet()) - anima::writeImage (outArg.getValue(),outputImage); - else - { - using MaskImageType = itk::Image ; - using CastFilterType = itk::CastImageFilter ; - - CastFilterType::Pointer castFilter = CastFilterType::New(); - castFilter->SetInput(outputImage); - castFilter->Update(); - - anima::writeImage (outArg.getValue(), castFilter->GetOutput()); - } - - return EXIT_SUCCESS; -} diff --git a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx new file mode 100644 index 000000000..a6f9e416c --- /dev/null +++ b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx @@ -0,0 +1,237 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include + + + +/** + * @brief Increment voxels between two consecutive track points. + * @param index Index of the voxel of the current point + * @param indexNext Index of the voxel of the following point + * @param incrementTerm Term which will be added to each interpolate voxel (either 1, or 1/total_number_of_tracks if proportionArg is true) + * @param outImage Pointer to the output image, used to add incrementTerm to each interpolate voxel + * @param checkInterpolateIsInsideImage If true, we must check that interpolate voxels are inside image before adding incrementTerm + * @details Add interpolation voxels if two consecutive track points are separated by at least one voxel. Then increment the number or proprtion of tracks + * passing through these voxels. + */ +template +void incrementInterpolateVoxels(IndexType index, IndexType indexNext, double incrementTerm, itk::Image ::Pointer outImage, + bool checkInterpolataleIsInsideImage = false) +{ + IndexType indexTemp; //where the interpolate voxels' index will be stored + itk::Image ::RegionType region = outImage->GetLargestPossibleRegion(); //the whole image + bool insideImage = true; + + + double dx = indexNext[0] - index[0]; //difference in x coordinate + double dy = indexNext[1] - index[1]; //difference in y coordinate + double dz = indexNext[2] - index[2]; //difference in z coordinate + double maxStep = std::max(std::max(std::abs(dx), std::abs(dy)), std::abs(dz)); //maximum between the three absolute differences + for (int l = 1; l < maxStep; l++) + { + //We enter this loop only if maxStep >= 2. + //In that case, we create interpolation voxels, to try to modify all and only all the voxels through which the track passes between point and pointNext + //We repeat maxStep - 1 times and not maxStep times, because otherwise, index = indexNext, and we would increment indexNext twice + indexTemp[0] = index[0] + std::round(l*dx/maxStep); + indexTemp[1] = index[1] + std::round(l*dy/maxStep); + indexTemp[2] = index[2] + std::round(l*dz/maxStep); + + if (checkInterpolataleIsInsideImage) + { + unsigned int k = 0; + while (insideImage && k<3) + { + insideImage = (indexTemp[k] < 0) || (indexTemp[k] >= region.GetSize(k)); + ++k; + } + //insideImage is true if and only if the index of the interpolate voxel is within the image + } + + //increment the number or proportion of tracks passing through the voxel with the index indexTemp, if it is inside the image + if (insideImage) + { + outImage->SetPixel(indexTemp, outImage->GetPixel(indexTemp) + incrementTerm); + } + } +} + +/** + * @brief Convert an input track file (vtp or vtk) into a density map. + * @details The density map is a 3D image giving for each voxel the number or proportion of tracks from the input image passing through it. + */ +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + + //Define arguments + TCLAP::ValueArg inArg("i", "input", "Input tracks file (vtp, vtk)", true, "", "input tracks file (vtp, vtk)", cmd); + TCLAP::ValueArg outArg("o", "output", "Output density map", true, "", "output density map file", cmd); + TCLAP::ValueArg geomArg("g", "geometry", "Geometry image", true, "", "geometry image, file with ext .nii.gz or similar", cmd); + + TCLAP::SwitchArg proportionArg("P", "proportion", "Output proportion of tracks going through each voxel", cmd, false); + TCLAP::SwitchArg noInterpolationArg("N", "no-interpolatation", "No creation of interpolation voxels if two consecutive track points are separated by at least one voxel", cmd, false); + + + //Try to parse + try + { + cmd.parse(argc, argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + //Define types + using ImageType = itk::Image ; + using OutImageType = itk::Image ; + + //Read geometry image and allocate memory for output image + ImageType::Pointer geomImage = anima::readImage (geomArg.getValue()); + OutImageType::Pointer outImage = OutImageType::New(); + + //Give to the output image the same properties as the geometry image + outImage->Initialize(); + outImage->SetRegions(geomImage->GetLargestPossibleRegion()); + outImage->SetSpacing(geomImage->GetSpacing()); + outImage->SetOrigin(geomImage->GetOrigin()); + outImage->SetDirection(geomImage->GetDirection()); + + //Allocate the correct size for the output image and intialize it with 0 in each voxel + outImage->Allocate(); + outImage->FillBuffer(0.0); + + //Read input track file + anima::ShapesReader trackReader; + trackReader.SetFileName(inArg.getValue()); + trackReader.Update(); + + //Initialize a vtksmartpointer with the content of the track file, and store the total number of tracks (1 cell = 1 track) + vtkSmartPointer tracks = trackReader.GetOutput(); + vtkIdType nbCells = tracks->GetNumberOfCells(); + + //Set the increment term (1 if proportion arg is not set, 1/total_number_of_tracks else) + double incrementTerm = 1.0; + if (proportionArg.isSet()) + { + incrementTerm /= nbCells; + } + + //Set interpolateVoxels (true if and only if we need to add interpolate voxels) + bool interpolateVoxels = !noInterpolationArg.isSet(); + std::cout<GetCell(j); + vtkPoints *cellPts = cell->GetPoints(); + vtkIdType nbPts = cellPts->GetNumberOfPoints(); + + //Iterate through the track j's points + for (int i = 0; i < nbPts-1; ++i) + { + //Get the current point and the following one, write their coordinates in point and pointNext + cellPts->GetPoint(i, ptVals); + cellPts->GetPoint(i+1, ptValsNext); + for (unsigned int k = 0; k < 3; ++k) + { + point[k] = ptVals[k]; + pointNext[k] = ptValsNext[k]; + } + + //TransformPhysicalPointToIndex enables to move from point coordinates to voxel index inside the output image. It returns true if the point is inside the image and false otherwise (in theory it should not happen) + if (outImage->TransformPhysicalPointToIndex(point, index)) + { + if (!outImage->TransformPhysicalPointToIndex(pointNext, indexNext)) + { + //The current point is inside the image but not the following one + outImage->SetPixel(index, outImage->GetPixel(index) + incrementTerm); //increment the current point's voxel + if (interpolateVoxels) + { + //Increment voxels between index of the current point and index of the following one + //As pointNext is not inside the image, we may create interpolate voxels which are not neither. So we put set last argument to true to check that interpolateVoxels are inside the image within the function + incrementInterpolateVoxels(index, indexNext, incrementTerm, outImage, true); + } + } + + //Both current point and following one are inside the image (common case) + + //if index==indexNext, the current point and the following one are in the same voxel + //so we don't increment the number or proportion of tracks passing through this voxel for now: we will do it in the next iteration + //that allows not to count one track several times in a voxel + if (index != indexNext) + { + //increment the current point's voxel, and interpolation voxels if interpolateVoxels is true + outImage->SetPixel(index, outImage->GetPixel(index) + incrementTerm); + if (interpolateVoxels) + { + incrementInterpolateVoxels(index, indexNext, incrementTerm, outImage); + } + } + } + else + { + //Current point is not inside the image, but following point is + if (outImage->TransformPhysicalPointToIndex(pointNext, indexNext)) + { + if (interpolateVoxels) + { + //Increment interpolate voxels, but need to check if they are inside the image (so last argument set to true) + incrementInterpolateVoxels(index, indexNext, incrementTerm, outImage, true); + } + } + + //If neither current point, nor the following one are inside the image, we have nothing to do in this iteration. + + } + } + + //Before moving to the next track, it remains to consider the last point of this track (its voxel has not been incremented yet) + cellPts->GetPoint(nbPts-1, ptVals); + for (unsigned int k = 0; k < 3; ++k) + { + point[k] = ptVals[k]; + } + if (outImage->TransformPhysicalPointToIndex(point, index)) + { + outImage->SetPixel(index, outImage->GetPixel(index) + incrementTerm); + } + } + + + //Write output image and exit + if (proportionArg.isSet()) + { + anima::writeImage (outArg.getValue(), outImage); + } + else + { + //We should convert the output image type from double to int first + using UIntImageType = itk::Image ; + using CastFilterType = itk::CastImageFilter ; + + CastFilterType::Pointer castFilter = CastFilterType::New(); + castFilter->SetInput(outImage); + castFilter->Update(); + + anima::writeImage (outArg.getValue(), castFilter->GetOutput()); + } + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Anima/math-tools/common_tools/CMakeLists.txt b/Anima/math-tools/common_tools/CMakeLists.txt index 63f091afb..34de43f19 100644 --- a/Anima/math-tools/common_tools/CMakeLists.txt +++ b/Anima/math-tools/common_tools/CMakeLists.txt @@ -6,7 +6,6 @@ add_subdirectory(collapse_image) add_subdirectory(convert_image) add_subdirectory(convert_shape) add_subdirectory(merge_block_images) -add_subdirectory(shapes_to_density) add_subdirectory(vectorize_images) add_subdirectory(image_to_mask) diff --git a/Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt b/Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt deleted file mode 100644 index f0dc72832..000000000 --- a/Anima/math-tools/common_tools/shapes_to_density/CMakeLists.txt +++ /dev/null @@ -1,38 +0,0 @@ -if(BUILD_TOOLS) - -project(animaShapesToDensity) - -## ############################################################################# -## List Sources -## ############################################################################# - -list_source_files(${PROJECT_NAME} - ${CMAKE_CURRENT_SOURCE_DIR} - ) - -## ############################################################################# -## add executable -## ############################################################################# - -add_executable(${PROJECT_NAME} - ${${PROJECT_NAME}_CFILES} - ) - - -## ############################################################################# -## Link -## ############################################################################# - -target_link_libraries(${PROJECT_NAME} - ${VTK_LIBRARIES} - ${ITKIO_LIBRARIES} - AnimaDataIO - ) - -## ############################################################################# -## install -## ############################################################################# - -set_exe_install_rules(${PROJECT_NAME}) - -endif() diff --git a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx b/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx deleted file mode 100644 index 9604a88c4..000000000 --- a/Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx +++ /dev/null @@ -1,135 +0,0 @@ -#include -#include -#include -#include - -#include - -#include -#include - -/** - * @brief Converts an input track file (vtp or vtk) into a density map. - * @details The density map is a 3D image giving for each voxel the number of streamlines from the input image passing through it. - * - */ -int main(int argc, char **argv) -{ - TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); - - //Define arguments - TCLAP::ValueArg inTrackArg("i","in-tracks","input track file (vtp, vtk)",true,"","input tracks",cmd); - TCLAP::ValueArg geometryArg("g","geometry","Output image geometry",true,"","output geometry",cmd); - TCLAP::ValueArg outArg("o","out","Output density map",true,"","output image",cmd); - - //Try to parse - try - { - cmd.parse(argc,argv); - } - catch (TCLAP::ArgException& e) - { - std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; - return EXIT_FAILURE; - } - - //Define types - typedef itk::Image ImageType; - using OutImageType = itk::Image ; - - //Read geometry image and allocate memory for output image - ImageType::Pointer geomImage = anima::readImage (geometryArg.getValue()); - OutImageType::Pointer outImage = OutImageType::New(); - - //Give to the output image the same properties as the geometry image - outImage->Initialize(); - outImage->SetRegions(geomImage->GetLargestPossibleRegion()); - outImage->SetSpacing(geomImage->GetSpacing()); - outImage->SetOrigin(geomImage->GetOrigin()); - outImage->SetDirection(geomImage->GetDirection()); - - //Allocate the correct size for the output image and intialize it with 0 in each voxel - outImage->Allocate(); - outImage->FillBuffer(0.0); - - //Read input track file - anima::ShapesReader trackReader; - trackReader.SetFileName(inTrackArg.getValue()); - trackReader.Update(); - - //Initialize a vtksmartpointer with the content of the track file, and store the total number of streamlines (1 cell = 1 streamline) - vtkSmartPointer tracks = trackReader.GetOutput(); - vtkIdType nbCells = tracks->GetNumberOfCells(); - - //Initialize data objects for the main algorithm - double ptVals[3], ptValsNext[3]; - OutImageType::IndexType index, indexNext, indexTemp; - OutImageType::PointType point, pointNext; - - - //Main algorithm - for (int j = 0; j < nbCells; ++j) - { - //Get streamline j, its points and its number of points - vtkCell *cell = tracks->GetCell(j); - vtkPoints *cellPts = cell->GetPoints(); - vtkIdType nbPts = cellPts->GetNumberOfPoints(); - - //Iterate through the streamline j's points - for (int i = 0; i < nbPts-1; ++i) - { - //Get the current point and the following one, write their coordinates in point and pointNext - cellPts->GetPoint(i, ptVals); - cellPts->GetPoint(i+1, ptValsNext); - for (unsigned int k = 0; k < 3; ++k) - { - point[k] = ptVals[k]; - pointNext[k] = ptValsNext[k]; - } - - //Infer the points' index in the output image (namely voxel coordinates) from the points' real coordinates. Write the results in index and indexNext. - outImage->TransformPhysicalPointToIndex(point, index); - outImage->TransformPhysicalPointToIndex(pointNext, indexNext); - - //if index==indexNext, the current point and the following one are in the same voxel - //so we don't increment the number of streamlines passing through this voxel for now: we will do it in the next iteration - //that allows not to count one streamline several times in a voxel - if (index != indexNext) - { - //increment the number of streamlines passing through the voxel with this index - outImage->SetPixel(index, outImage->GetPixel(index) + 1.0); - - double dx = indexNext[0] - index[0]; //difference in x coordinate - double dy = indexNext[1] - index[1]; //difference in y coordinate - double dz = indexNext[2] - index[2]; //difference in z coordinate - double maxStep = std::max(std::max(std::abs(dx), std::abs(dy)), std::abs(dz)); //maximum between the three absolute differences - for (int l = 1; l < maxStep; l++) - { - //We enter this loop only if maxStep >=2. - //In that case, we create artificial index points, - // to try to modify all and only all the voxels through which the streamline passes between point and pointNext - //We repeat maxStep - 1 times and not maxStep times, because otherwise, index = indexNext, and we would increment indexNext twice - indexTemp[0] = index[0] + std::round(l*dx/maxStep); - indexTemp[1] = index[1] + std::round(l*dy/maxStep); - indexTemp[2] = index[2] + std::round(l*dz/maxStep); - //increment the number of streamlines passing through the voxel with the index indexTemp - outImage->SetPixel(indexTemp, outImage->GetPixel(indexTemp) + 1.0); - } - } - } - - //before moving to the next streamline, it remains to consider the last point of this streamline (its voxel has not been incremented yet) - cellPts->GetPoint(nbPts-1, ptVals); - for (unsigned int k = 0; k < 3; ++k) - { - point[k] = ptVals[k]; - } - outImage->TransformPhysicalPointToIndex(point, index); - outImage->SetPixel(index, outImage->GetPixel(index) + 1.0); - } - - - //Write output image and exit - anima::writeImage (outArg.getValue(), outImage); - return EXIT_SUCCESS; -} From f9fbc36fd452ae56f8bcdcb31d9c6ec16647dbca Mon Sep 17 00:00:00 2001 From: Gregoire Ville Date: Fri, 6 Dec 2024 10:42:22 +0100 Subject: [PATCH 7/8] Removal of a forgotten std::cout --- .../tractography/fibers_counter/animaFibersCounter_fusion.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx index a6f9e416c..42f668253 100644 --- a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx +++ b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx @@ -127,7 +127,6 @@ int main(int argc, char **argv) //Set interpolateVoxels (true if and only if we need to add interpolate voxels) bool interpolateVoxels = !noInterpolationArg.isSet(); - std::cout< Date: Tue, 10 Dec 2024 17:54:12 +0100 Subject: [PATCH 8/8] Rename animaFibersCounter_fusion.cxx and add of a reminder about orientation convention change. --- .../{animaFibersCounter_fusion.cxx => animaFibersCounter.cxx} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename Anima/diffusion/tractography/fibers_counter/{animaFibersCounter_fusion.cxx => animaFibersCounter.cxx} (95%) diff --git a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx similarity index 95% rename from Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx rename to Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx index 42f668253..334bccf5e 100644 --- a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter_fusion.cxx +++ b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx @@ -68,7 +68,7 @@ void incrementInterpolateVoxels(IndexType index, IndexType indexNext, double inc */ int main(int argc, char **argv) { - TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); + TCLAP::CmdLine cmd("This binary converts an input tracks file (vtp or vtk) into a density map, giving for each voxel the number or proportion of tracks from the input image passing through it. \n Warning : Anima uses LPS convention; if input file uses another convention (for example RAS) one must convert to LPS convention first. \n For example to move from RAS to LPS convention, one must rotate the input tracks file by 180° around z-axis. \n INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION); //Define arguments TCLAP::ValueArg inArg("i", "input", "Input tracks file (vtp, vtk)", true, "", "input tracks file (vtp, vtk)", cmd);