diff --git a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx index 63a8a2a25..334bccf5e 100644 --- a/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx +++ b/Anima/diffusion/tractography/fibers_counter/animaFibersCounter.cxx @@ -1,29 +1,88 @@ +#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; -#include -#include -#include -#include -#include -#include -#include + 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 + } -#include + //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("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::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); - 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); + //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 fibers going through each pixel",cmd,false); + 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); + cmd.parse(argc, argv); } catch (TCLAP::ArgException& e) { @@ -31,73 +90,147 @@ int main(int argc, char **argv) return EXIT_FAILURE; } - typedef itk::Image OutputImageType; - OutputImageType::Pointer outputImage = anima::readImage (geomArg.getValue()); - outputImage->FillBuffer(0.0); + //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(); - vtkSmartPointer tracks = trackReader.GetOutput(); - + //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(); - double incrementFactor = 1.0; + + //Set the increment term (1 if proportion arg is not set, 1/total_number_of_tracks else) + double incrementTerm = 1.0; if (proportionArg.isSet()) - incrementFactor /= nbCells; + { + incrementTerm /= nbCells; + } + + //Set interpolateVoxels (true if and only if we need to add interpolate voxels) + bool interpolateVoxels = !noInterpolationArg.isSet(); - double ptVals[3]; - OutputImageType::IndexType currentIndex; - OutputImageType::PointType currentPoint; - OutputImageType::RegionType region = outputImage->GetLargestPossibleRegion(); + //Initialize data objects for the main algorithm + double ptVals[3], ptValsNext[3]; + OutImageType::IndexType index, indexNext, indexTemp; + OutImageType::PointType point, pointNext; - // Explores individual fibers - for (int j = 0;j < nbCells;++j) + + //Main algorithm + for (int j = 0; j < nbCells; ++j) { + //Get track j, its points and its number of points vtkCell *cell = tracks->GetCell(j); vtkPoints *cellPts = cell->GetPoints(); vtkIdType nbPts = cellPts->GetNumberOfPoints(); - // Explores points in fibers - for (int i = 0;i < nbPts;++i) + //Iterate through the track j's points + for (int i = 0; i < nbPts-1; ++i) { - cellPts->GetPoint(i,ptVals); - - for (unsigned int k = 0;k < 3;++k) - currentPoint[k] = ptVals[k]; + //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]; + } - outputImage->TransformPhysicalPointToIndex(currentPoint,currentIndex); - bool insideBuffer = true; - for (unsigned int k = 0;k < 3;++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 ((currentIndex[k] < 0)||(currentIndex[k] >= region.GetSize(k))) + 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) { - insideBuffer = false; - break; + //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. - if (!insideBuffer) - continue; + } + } - double countIndex = outputImage->GetPixel(currentIndex) + incrementFactor; - outputImage->SetPixel(currentIndex, countIndex); + //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(),outputImage); + { + anima::writeImage (outArg.getValue(), outImage); + } else { - using MaskImageType = itk::Image ; - using CastFilterType = itk::CastImageFilter ; + //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(outputImage); + castFilter->SetInput(outImage); castFilter->Update(); - anima::writeImage (outArg.getValue(), castFilter->GetOutput()); + anima::writeImage (outArg.getValue(), castFilter->GetOutput()); } return EXIT_SUCCESS; -} +} \ No newline at end of file