Feature to convert a tracks file to a density map

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)
 

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()

Anima/math-tools/common_tools/shapes_to_density/animaShapesToDensity.cxx

@@ -0,0 +1,111 @@
+#include <animaShapesReader.h>
+#include <animaReadWriteFunctions.h>
+#include <itkVectorImage.h>
+#include <itkImageIOFactory.h>
+
+#include <vtkVector.h>
+
+#include <tclap/CmdLine.h>
+#include <fstream>
+
+template <typename T> 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 <std::string> inTrackArg("i","in-tracks","input track file (vtp, vtk)",true,"","input tracks",cmd);
+    TCLAP::ValueArg <std::string> geometryArg("g","geometry","Output image geometry",true,"","output geometry",cmd);
+    TCLAP::ValueArg <std::string> 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 <double,3> ImageType;
+
+    ImageType::Pointer geomImage = anima::readImage <ImageType> (geometryArg.getValue());
+
+    using OutImageType = itk::Image <double, 3>;
+
+//    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<vtkPolyData> 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 <OutImageType> (outArg.getValue(), outImage);
+
+    return EXIT_SUCCESS;
+}