diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp index 4f2333d592..49a67a5f82 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp @@ -1,418 +1,417 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkParticleGrid.h" #include #include using namespace mitk; ParticleGrid::ParticleGrid(ItkFloatImageType* image, float particleLength, int cellCapacity) { // initialize counters m_NumParticles = 0; m_NumConnections = 0; m_NumCellOverflows = 0; m_ParticleLength = particleLength; // define isotropic grid from voxel spacing and particle length float cellSize = 2*m_ParticleLength; m_GridSize[0] = image->GetLargestPossibleRegion().GetSize()[0]*image->GetSpacing()[0]/cellSize +1; m_GridSize[1] = image->GetLargestPossibleRegion().GetSize()[1]*image->GetSpacing()[1]/cellSize +1; m_GridSize[2] = image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]/cellSize +1; m_GridScale[0] = 1/cellSize; m_GridScale[1] = 1/cellSize; m_GridScale[2] = 1/cellSize; m_CellCapacity = cellCapacity; // maximum number of particles per grid cell m_ContainerCapacity = 100000; // initial particle container capacity unsigned long numCells = m_GridSize[0]*m_GridSize[1]*m_GridSize[2]; // number of grid cells unsigned long gridSize = numCells*m_CellCapacity; if ( itk::NumericTraits::max() R) { if (m_NumParticles >= m_ContainerCapacity) { if (!ReallocateGrid()) return NULL; } int xint = int(R[0]*m_GridScale[0]); if (xint < 0) return NULL; if (xint >= m_GridSize[0]) return NULL; int yint = int(R[1]*m_GridScale[1]); if (yint < 0) return NULL; if (yint >= m_GridSize[1]) return NULL; int zint = int(R[2]*m_GridScale[2]); if (zint < 0) return NULL; if (zint >= m_GridSize[2]) return NULL; int idx = xint + m_GridSize[0]*(yint + m_GridSize[1]*zint); if (m_OccupationCount[idx] < m_CellCapacity) { Particle *p = &(m_Particles[m_NumParticles]); p->pos = R; p->mID = -1; p->pID = -1; m_NumParticles++; p->gridindex = m_CellCapacity*idx + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return p; } else { m_NumCellOverflows++; return NULL; } } bool ParticleGrid::TryUpdateGrid(int k) { Particle* p = &(m_Particles[k]); int xint = int(p->pos[0]*m_GridScale[0]); if (xint < 0) return false; if (xint >= m_GridSize[0]) return false; int yint = int(p->pos[1]*m_GridScale[1]); if (yint < 0) return false; if (yint >= m_GridSize[1]) return false; int zint = int(p->pos[2]*m_GridScale[2]); if (zint < 0) return false; if (zint >= m_GridSize[2]) return false; int idx = xint + m_GridSize[0]*(yint+ zint*m_GridSize[1]); int cellidx = p->gridindex/m_CellCapacity; if (idx != cellidx) // cell has changed { if (m_OccupationCount[idx] < m_CellCapacity) { // remove from old position in grid; int grdindex = p->gridindex; m_Grid[grdindex] = m_Grid[cellidx*m_CellCapacity + m_OccupationCount[cellidx]-1]; m_Grid[grdindex]->gridindex = grdindex; m_OccupationCount[cellidx]--; // insert at new position in grid p->gridindex = idx*m_CellCapacity + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return true; } else { m_NumCellOverflows++; return false; } } return true; } void ParticleGrid::RemoveParticle(int k) { Particle* p = &(m_Particles[k]); int gridIndex = p->gridindex; int cellIdx = gridIndex/m_CellCapacity; int idx = gridIndex%m_CellCapacity; // remove pending connections if (p->mID != -1) DestroyConnection(p,-1); if (p->pID != -1) DestroyConnection(p,+1); // remove from grid if (idx < m_OccupationCount[cellIdx]-1) { m_Grid[gridIndex] = m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1]; m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1] = NULL; m_Grid[gridIndex]->gridindex = gridIndex; } m_OccupationCount[cellIdx]--; // remove from container if (k < m_NumParticles-1) { Particle* last = &m_Particles[m_NumParticles-1]; // last particle // update connections of last particle because its index is changing if (last->mID!=-1) { if ( m_Particles[last->mID].mID == m_NumParticles-1 ) m_Particles[last->mID].mID = k; else if ( m_Particles[last->mID].pID == m_NumParticles-1 ) m_Particles[last->mID].pID = k; } if (last->pID!=-1) { if ( m_Particles[last->pID].mID == m_NumParticles-1 ) m_Particles[last->pID].mID = k; else if ( m_Particles[last->pID].pID == m_NumParticles-1 ) m_Particles[last->pID].pID = k; } m_Particles[k] = m_Particles[m_NumParticles-1]; // move very last particle to empty slot m_Particles[m_NumParticles-1].ID = m_NumParticles-1; // update ID of removed particle to match the index m_Particles[k].ID = k; // update ID of moved particle m_Grid[m_Particles[k].gridindex] = &m_Particles[k]; // update address of moved particle } m_NumParticles--; } void ParticleGrid::ComputeNeighbors(vnl_vector_fixed &R) { float xfrac = R[0]*m_GridScale[0]; float yfrac = R[1]*m_GridScale[1]; float zfrac = R[2]*m_GridScale[2]; int xint = int(xfrac); int yint = int(yfrac); int zint = int(zfrac); int dx = -1; if (xfrac-xint > 0.5) dx = 1; if (xint <= 0) { xint = 0; dx = 1; } if (xint >= m_GridSize[0]-1) { xint = m_GridSize[0]-1; dx = -1; } int dy = -1; if (yfrac-yint > 0.5) dy = 1; if (yint <= 0) {yint = 0; dy = 1; } if (yint >= m_GridSize[1]-1) {yint = m_GridSize[1]-1; dy = -1;} int dz = -1; if (zfrac-zint > 0.5) dz = 1; if (zint <= 0) {zint = 0; dz = 1; } if (zint >= m_GridSize[2]-1) {zint = m_GridSize[2]-1; dz = -1;} m_NeighbourTracker.cellidx[0] = xint + m_GridSize[0]*(yint+zint*m_GridSize[1]); m_NeighbourTracker.cellidx[1] = m_NeighbourTracker.cellidx[0] + dx; m_NeighbourTracker.cellidx[2] = m_NeighbourTracker.cellidx[1] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[3] = m_NeighbourTracker.cellidx[2] - dx; m_NeighbourTracker.cellidx[4] = m_NeighbourTracker.cellidx[0] + dz*m_GridSize[0]*m_GridSize[1]; m_NeighbourTracker.cellidx[5] = m_NeighbourTracker.cellidx[4] + dx; m_NeighbourTracker.cellidx[6] = m_NeighbourTracker.cellidx[5] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[7] = m_NeighbourTracker.cellidx[6] - dx; m_NeighbourTracker.cellidx_c[0] = m_CellCapacity*m_NeighbourTracker.cellidx[0]; m_NeighbourTracker.cellidx_c[1] = m_CellCapacity*m_NeighbourTracker.cellidx[1]; m_NeighbourTracker.cellidx_c[2] = m_CellCapacity*m_NeighbourTracker.cellidx[2]; m_NeighbourTracker.cellidx_c[3] = m_CellCapacity*m_NeighbourTracker.cellidx[3]; m_NeighbourTracker.cellidx_c[4] = m_CellCapacity*m_NeighbourTracker.cellidx[4]; m_NeighbourTracker.cellidx_c[5] = m_CellCapacity*m_NeighbourTracker.cellidx[5]; m_NeighbourTracker.cellidx_c[6] = m_CellCapacity*m_NeighbourTracker.cellidx[6]; m_NeighbourTracker.cellidx_c[7] = m_CellCapacity*m_NeighbourTracker.cellidx[7]; m_NeighbourTracker.cellcnt = 0; m_NeighbourTracker.pcnt = 0; } Particle* ParticleGrid::GetNextNeighbor() { if (m_NeighbourTracker.pcnt < m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]]) { return m_Grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt] + (m_NeighbourTracker.pcnt++)]; } else { for(;;) { m_NeighbourTracker.cellcnt++; if (m_NeighbourTracker.cellcnt >= 8) return 0; if (m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]] > 0) break; } m_NeighbourTracker.pcnt = 1; return m_Grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt]]; } } void ParticleGrid::CreateConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = P2->ID; else P1->pID = P2->ID; if (ep2 == -1) P2->mID = P1->ID; else P2->pID = P1->ID; m_NumConnections++; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = -1; else P1->pID = -1; if (ep2 == -1) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1) { Particle *P2 = 0; if (ep1 == 1) { P2 = &m_Particles[P1->pID]; P1->pID = -1; } else { P2 = &m_Particles[P1->mID]; P1->mID = -1; } if (P2->mID == P1->ID) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } bool ParticleGrid::CheckConsistency() { for (int i=0; iID != i) { std::cout << "Particle ID error!" << std::endl; return false; } if (p->mID!=-1) { Particle* p2 = &m_Particles[p->mID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } if (p->pID!=-1) { Particle* p2 = &m_Particles[p->pID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } } return true; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp new file mode 100755 index 0000000000..5fc4b37e71 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.cpp @@ -0,0 +1,377 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef __itkEvaluateDirectionImagesFilter_cpp +#define __itkEvaluateDirectionImagesFilter_cpp + +#include "itkEvaluateDirectionImagesFilter.h" +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +namespace itk { + +template< class PixelType > +EvaluateDirectionImagesFilter< PixelType > +::EvaluateDirectionImagesFilter(): + m_ImageSet(NULL), + m_ReferenceImageSet(NULL), + m_IgnoreMissingDirections(false), + m_Eps(0.0001) +{ + this->SetNumberOfOutputs(2); +} + +template< class PixelType > +void EvaluateDirectionImagesFilter< PixelType >::GenerateData() +{ + if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) + return; + + DirectionImageContainerType::Pointer set1 = DirectionImageContainerType::New(); + DirectionImageContainerType::Pointer set2 = DirectionImageContainerType::New(); + + for (int i=0; iSize(); i++) + { + typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); + duplicator->SetInputImage( m_ImageSet->GetElement(i) ); + duplicator->Update(); + set1->InsertElement(i, dynamic_cast(duplicator->GetOutput())); + } + for (int i=0; iSize(); i++) + { + typename itk::ImageDuplicator< DirectionImageType >::Pointer duplicator = itk::ImageDuplicator< DirectionImageType >::New(); + duplicator->SetInputImage( m_ReferenceImageSet->GetElement(i) ); + duplicator->Update(); + set2->InsertElement(i, dynamic_cast(duplicator->GetOutput())); + } + + m_ImageSet = set1; + m_ReferenceImageSet = set2; + + // angular error image + typename OutputImageType::Pointer outputImage = OutputImageType::New(); + outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + outputImage->Allocate(); + outputImage->FillBuffer(0.0); + this->SetNthOutput(0, outputImage); + + // length error image + outputImage = OutputImageType::New(); + outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + outputImage->Allocate(); + outputImage->FillBuffer(0.0); + this->SetNthOutput(1, outputImage); + + if (m_MaskImage.IsNull()) + { + m_MaskImage = UCharImageType::New(); + m_MaskImage->SetOrigin( outputImage->GetOrigin() ); + m_MaskImage->SetRegions( outputImage->GetLargestPossibleRegion() ); + m_MaskImage->SetSpacing( outputImage->GetSpacing() ); + m_MaskImage->SetDirection( outputImage->GetDirection() ); + m_MaskImage->Allocate(); + m_MaskImage->FillBuffer(1); + } + + m_MeanAngularError = 0.0; + m_MedianAngularError = 0; + m_MaxAngularError = 0.0; + m_MinAngularError = itk::NumericTraits::max(); + m_VarAngularError = 0.0; + m_AngularErrorVector.clear(); + + m_MeanLengthError = 0.0; + m_MedianLengthError = 0; + m_MaxLengthError = 0.0; + m_MinLengthError = itk::NumericTraits::max(); + m_VarLengthError = 0.0; + m_LengthErrorVector.clear(); + + if (m_ImageSet.IsNull() || m_ReferenceImageSet.IsNull()) + return; + + outputImage = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(0)); + typename OutputImageType::Pointer outputImage2 = static_cast< OutputImageType* >(this->ProcessObject::GetOutput(1)); + + ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); + ImageRegionIterator< OutputImageType > oit2(outputImage2, outputImage2->GetLargestPossibleRegion()); + ImageRegionIterator< UCharImageType > mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion()); + + int numImages = m_ImageSet->Size(); + int numReferenceImages = m_ReferenceImageSet->Size(); + + // fill missing directions with zeros + if (numImages>numReferenceImages) + { + DirectionType zeroDir; zeroDir.Fill(0.0); + for (int i=0; iSetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + img->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + img->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + img->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + img->Allocate(); + img->FillBuffer(zeroDir); + m_ReferenceImageSet->InsertElement(m_ReferenceImageSet->Size(), img); + } + numReferenceImages = numImages; + } + else if (numReferenceImages>numImages) + { + DirectionType zeroDir; zeroDir.Fill(0.0); + for (int i=0; iSetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + img->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + img->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + img->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + img->Allocate(); + img->FillBuffer(zeroDir); + m_ImageSet->InsertElement(m_ImageSet->Size(), img); + } + numImages = numReferenceImages; + } + int numDirections = numReferenceImages; + + // matrix containing the angular error between the directions + vnl_matrix< float > diffM; diffM.set_size(numDirections, numDirections); + + boost::progress_display disp(outputImage->GetLargestPossibleRegion().GetSize()[0]*outputImage->GetLargestPossibleRegion().GetSize()[1]*outputImage->GetLargestPossibleRegion().GetSize()[2]); + while( !oit.IsAtEnd() ) + { + ++disp; + if( mit.Get()!=1 ) + { + ++oit; + ++oit2; + ++mit; + continue; + } + typename OutputImageType::IndexType index = oit.GetIndex(); + + float maxAngularError = 1.0; + + diffM.fill(10); // initialize with invalid error value + + // get number of valid directions (length > 0) + int numRefDirs = 0; + int numTestDirs = 0; + for (int i=0; iGetElement(i)->GetPixel(index).GetVnlVector().magnitude() > m_Eps ) + numRefDirs++; + if (m_ImageSet->GetElement(i)->GetPixel(index).GetVnlVector().magnitude() > m_Eps ) + numTestDirs++; + } + + // i: index of reference direction + // j: index of test direction + for (int i=0; i refDir = m_ReferenceImageSet->GetElement(i)->GetPixel(index).GetVnlVector(); + + if (i testDir = m_ImageSet->GetElement(j)->GetPixel(index).GetVnlVector(); + if (j1.0) + diffM[i][j] = 1.0; + } + } + + float angularError = 0.0; + float lengthError = 0.0; + int counter = 0; + vnl_matrix< float > diffM_copy = diffM; + for (int k=0; k small error) + for (int i=0; ierror && diffM[i][j]<2) // found valid error entry + { + error = diffM[i][j]; + a = i; + b = j; + missingDir = false; + } + else if (diffM[i][j]<0 && error<0) // found missing direction + { + a = i; + b = j; + missingDir = true; + } + } + + if (a<0 || b<0 || m_IgnoreMissingDirections && missingDir) + continue; // no more directions found + + if (a>=numRefDirs && b>=numTestDirs) + { + MITK_INFO << "ERROR: missing test and reference direction. should not be possible. check code."; + continue; + } + + // remove processed directions from error matrix + diffM.set_row(a, 10.0); + diffM.set_column(b, 10.0); + + if (a>=numRefDirs) // missing reference direction (find next closest) + { + for (int i=0; ierror) + { + error = diffM_copy[i][b]; + a = i; + } + } + else if (b>=numTestDirs) // missing test direction (find next closest) + { + for (int i=0; ierror) + { + error = diffM_copy[a][i]; + b = i; + } + } + + float refLength = m_ReferenceImageSet->GetElement(a)->GetPixel(index).GetVnlVector().magnitude(); + float testLength = m_ImageSet->GetElement(b)->GetPixel(index).GetVnlVector().magnitude(); + + if (a>=numRefDirs || b>=numTestDirs || error<0) + error = 0; + + m_LengthErrorVector.push_back( fabs(refLength-testLength) ); + m_AngularErrorVector.push_back( acos(error)*180.0/M_PI ); + + m_MeanAngularError += m_AngularErrorVector.back(); + m_MeanLengthError += m_LengthErrorVector.back(); + + angularError += m_AngularErrorVector.back(); + lengthError += m_LengthErrorVector.back(); + counter++; + } + + if (counter>0) + { + lengthError /= counter; + angularError /= counter; + } + oit2.Set(lengthError); + oit.Set(angularError); + + ++oit; + ++oit2; + ++mit; + } + + std::sort( m_AngularErrorVector.begin(), m_AngularErrorVector.end() ); + m_MeanAngularError /= m_AngularErrorVector.size(); // mean + for (int i=0; im_MaxAngularError ) + m_MaxAngularError = m_AngularErrorVector.at(i); + if ( m_AngularErrorVector.at(i)1) + { + m_VarAngularError /= (m_AngularErrorVector.size()-1); // variance + + // median + if (m_AngularErrorVector.size()%2 == 0) + m_MedianAngularError = 0.5*( m_AngularErrorVector.at( m_AngularErrorVector.size()/2 ) + m_AngularErrorVector.at( m_AngularErrorVector.size()/2+1 ) ); + else + m_MedianAngularError = m_AngularErrorVector.at( (m_AngularErrorVector.size()+1)/2 ) ; + } + + std::sort( m_LengthErrorVector.begin(), m_LengthErrorVector.end() ); + m_MeanLengthError /= m_LengthErrorVector.size(); // mean + for (int i=0; im_MaxLengthError ) + m_MaxLengthError = m_LengthErrorVector.at(i); + if ( m_LengthErrorVector.at(i)1) + { + m_VarLengthError /= (m_LengthErrorVector.size()-1); // variance + + // median + if (m_LengthErrorVector.size()%2 == 0) + m_MedianLengthError = 0.5*( m_LengthErrorVector.at( m_LengthErrorVector.size()/2 ) + m_LengthErrorVector.at( m_LengthErrorVector.size()/2+1 ) ); + else + m_MedianLengthError = m_LengthErrorVector.at( (m_LengthErrorVector.size()+1)/2 ) ; + } +} + +} + +#endif // __itkEvaluateDirectionImagesFilter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h new file mode 100755 index 0000000000..964c17e5fc --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateDirectionImagesFilter.h @@ -0,0 +1,110 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +/*=================================================================== + +This file is based heavily on a corresponding ITK filter. + +===================================================================*/ +#ifndef __itkEvaluateDirectionImagesFilter_h_ +#define __itkEvaluateDirectionImagesFilter_h_ + +#include +#include +#include + +namespace itk{ +/** \class EvaluateDirectionImagesFilter + */ + +template< class PixelType > +class EvaluateDirectionImagesFilter : public ImageSource< Image< PixelType, 3 > > +{ + +public: + + typedef EvaluateDirectionImagesFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageSource< Image< PixelType, 3 > > Superclass; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + typedef typename Superclass::OutputImageType OutputImageType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(EvaluateDirectionImagesFilter, ImageToImageFilter) + + typedef Vector< float, 3 > DirectionType; + typedef Image< DirectionType, 3 > DirectionImageType; + typedef VectorContainer< int, DirectionImageType::Pointer > DirectionImageContainerType; + typedef Image< float, 3 > FloatImageType; + typedef Image< bool, 3 > BoolImageType; + typedef Image< unsigned char, 3 > UCharImageType; + + itkSetMacro( ImageSet , DirectionImageContainerType::Pointer) + itkSetMacro( ReferenceImageSet , DirectionImageContainerType::Pointer) + itkSetMacro( MaskImage , UCharImageType::Pointer) + itkSetMacro( IgnoreMissingDirections , bool) + + itkGetMacro( MeanAngularError, float) + itkGetMacro( MinAngularError, float) + itkGetMacro( MaxAngularError, float) + itkGetMacro( VarAngularError, float) + itkGetMacro( MedianAngularError, float) + + itkGetMacro( MeanLengthError, float) + itkGetMacro( MinLengthError, float) + itkGetMacro( MaxLengthError, float) + itkGetMacro( VarLengthError, float) + itkGetMacro( MedianLengthError, float) + +protected: + EvaluateDirectionImagesFilter(); + ~EvaluateDirectionImagesFilter() {} + + void GenerateData(); + + UCharImageType::Pointer m_MaskImage; + DirectionImageContainerType::Pointer m_ImageSet; + DirectionImageContainerType::Pointer m_ReferenceImageSet; + bool m_IgnoreMissingDirections; + double m_MeanAngularError; + double m_MedianAngularError; + double m_MaxAngularError; + double m_MinAngularError; + double m_VarAngularError; + std::vector< double > m_AngularErrorVector; + + double m_MeanLengthError; + double m_MedianLengthError; + double m_MaxLengthError; + double m_MinLengthError; + double m_VarLengthError; + std::vector< double > m_LengthErrorVector; + + double m_Eps; +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkEvaluateDirectionImagesFilter.cpp" +#endif + +#endif //__itkEvaluateDirectionImagesFilter_h_ + diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt index 860206a99f..8e7361bdee 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt @@ -1,47 +1,49 @@ OPTION(BUILD_FiberTrackingMiniApps "Build commandline tools for fiber tracking" OFF) IF(BUILD_FiberTrackingMiniApps) # include necessary modules here MitkExt QmitkExt MITK_CHECK_MODULE(_RESULT DiffusionCore FiberTracking ) IF(_RESULT) MESSAGE("Warning: FiberTrackingMiniApps is missing ${_RESULT}") ELSE(_RESULT) MITK_USE_MODULE( DiffusionCore FiberTracking ) # needed include directories INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ${ALL_INCLUDE_DIRECTORIES}) PROJECT( mitkFiberTrackingMiniApps ) # fill in the standalone executables here SET(FIBERTRACKINGMINIAPPS mitkFiberTrackingMiniApps ) # set additional files here SET(FIBERTRACKING_ADDITIONAL_FILES MiniAppManager.cpp GibbsTracking.cpp + StreamlineTracking.cpp + TractometerAngularErrorTool.cpp ) # create an executable foreach tool (only one at the moment) FOREACH(tool ${FIBERTRACKINGMINIAPPS}) ADD_EXECUTABLE( ${tool} ${tool}.cpp ${FIBERTRACKING_ADDITIONAL_FILES} ) TARGET_LINK_LIBRARIES( ${tool} ${ALL_LIBRARIES} ) ENDFOREACH(tool) ENDIF() MITK_INSTALL_TARGETS(EXECUTABLES mitkFiberTrackingMiniApps ) ENDIF(BUILD_FiberTrackingMiniApps) diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp index c527d8367d..f15d2bb328 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/GibbsTracking.cpp @@ -1,219 +1,217 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" #include #include #include #include #include #include #include #include #include #include #include +#include +#include "ctkCommandLineParser.h" template typename itk::ShCoefficientImageImporter< float, shOrder >::QballImageType::Pointer TemplatedConvertShCoeffs(mitk::Image* mitkImg, int toolkit) { typedef itk::ShCoefficientImageImporter< float, shOrder > FilterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImg); caster->Update(); typename FilterType::Pointer filter = FilterType::New(); switch (toolkit) { case 0: filter->SetToolkit(FilterType::FSL); break; case 1: filter->SetToolkit(FilterType::MRTRIX); break; default: filter->SetToolkit(FilterType::FSL); } filter->SetInputImage(caster->GetOutput()); filter->GenerateData(); return filter->GetQballImage(); } int GibbsTracking(int argc, char* argv[]) { - if (argc<6) - { - std::cout << std::endl; - std::cout << "Performes Gibbs Tracking" << std::endl; - std::cout << std::endl; - std::cout << "usage: " << argv[0] << " " << std::endl; - std::cout << std::endl; + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", ctkCommandLineParser::String, "input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image)", mitk::Any(), false); + parser.addArgument("parameters", "p", ctkCommandLineParser::String, "parameter file (.gtp)", mitk::Any(), false); + parser.addArgument("mask", "m", ctkCommandLineParser::String, "binary mask image"); + parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "sh coefficient convention (FSL, MRtrix) (default = FSL)"); + parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output fiber bundle (.fib)", mitk::Any(), false); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) return EXIT_FAILURE; - } + + string inFileName = mitk::any_cast(parsedArgs["input"]); + string paramFileName = mitk::any_cast(parsedArgs["parameters"]); + string outFileName = mitk::any_cast(parsedArgs["outFile"]); try { RegisterDiffusionCoreObjectFactory(); RegisterFiberTrackingObjectFactory(); - std::cout << "input: " << argv[1] << std::endl; - std::cout << "param: " << argv[2] << std::endl; - std::cout << "output: " << argv[3] << std::endl; - // instantiate gibbs tracker typedef itk::Vector OdfVectorType; typedef itk::Image ItkQballImageType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); // load input image const std::string s1="", s2=""; - std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[1], s1, s2, false ); + std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false ); // try to cast to qball image - QString inImageName(argv[1]); + QString inImageName(inFileName.c_str()); if( inImageName.endsWith(".qbi") ) { MITK_INFO << "Loading qball image ..."; mitk::QBallImage::Pointer mitkQballImage = dynamic_cast(infile.at(0).GetPointer()); ItkQballImageType::Pointer itk_qbi = ItkQballImageType::New(); mitk::CastToItkImage(mitkQballImage, itk_qbi); gibbsTracker->SetQBallImage(itk_qbi.GetPointer()); } else if( inImageName.endsWith(".dti") ) { MITK_INFO << "Loading tensor image ..."; typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.at(0).GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); gibbsTracker->SetTensorImage(itk_dti); } else if ( inImageName.endsWith(".nii") ) { MITK_INFO << "Loading sh-coefficient image ..."; mitk::Image::Pointer mitkImage = dynamic_cast(infile.at(0).GetPointer()); int nrCoeffs = mitkImage->GetLargestPossibleRegion().GetSize()[3]; int c=3, d=2-2*nrCoeffs; double D = c*c-4*d; int shOrder; if (D>0) { shOrder = (-c+sqrt(D))/2.0; if (shOrder<0) shOrder = (-c-sqrt(D))/2.0; } else if (D==0) shOrder = -c/2.0; MITK_INFO << "using SH-order " << shOrder; int toolkitConvention = 0; - if(argc==6) + + if (parsedArgs.count("shConvention")) { - QString convention(argv[5]); + QString convention(mitk::any_cast(parsedArgs["shConvention"]).c_str()); if ( convention.compare("MRtrix", Qt::CaseInsensitive)==0 ) { toolkitConvention = 1; MITK_INFO << "Using MRtrix style sh-coefficient convention"; } else MITK_INFO << "Using FSL style sh-coefficient convention"; } else MITK_INFO << "Using FSL style sh-coefficient convention"; switch (shOrder) { case 4: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<4>(mitkImage, toolkitConvention)); break; case 6: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<6>(mitkImage, toolkitConvention)); break; case 8: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<8>(mitkImage, toolkitConvention)); break; case 10: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<10>(mitkImage, toolkitConvention)); break; case 12: gibbsTracker->SetQBallImage(TemplatedConvertShCoeffs<12>(mitkImage, toolkitConvention)); break; default: MITK_INFO << "SH-order " << shOrder << " not supported"; } } else return EXIT_FAILURE; // global tracking - if(argc>=5) + if (parsedArgs.count("mask")) { typedef itk::Image MaskImgType; - std::cout << "mask: " << argv[4]<< std::endl; - infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[4], s1, s2, false ); - mitk::Image::Pointer mitkMaskImage = dynamic_cast(infile.at(0).GetPointer()); + mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::LoadImage(mitk::any_cast(parsedArgs["mask"])); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); gibbsTracker->SetMaskImage(itk_mask); } - else{ - std::cout << "perform tracking without mask" << std::endl; - } gibbsTracker->SetDuplicateImage(false); - gibbsTracker->SetLoadParameterFile( argv[2] ); + gibbsTracker->SetLoadParameterFile( paramFileName ); // gibbsTracker->SetLutPath( "" ); gibbsTracker->Update(); mitk::FiberBundleX::Pointer mitkFiberBundle = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); - std::string outfilename = argv[3]; - MITK_INFO << "searching writer"; mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) { if ( (*it)->CanWriteBaseDataType(mitkFiberBundle.GetPointer()) ) { - (*it)->SetFileName( outfilename.c_str() ); + (*it)->SetFileName( outFileName.c_str() ); (*it)->DoWrite( mitkFiberBundle.GetPointer() ); } } } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_INFO << "DONE"; return EXIT_SUCCESS; } RegisterFiberTrackingMiniApp(GibbsTracking); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.cpp index a27128840b..ff9b82b1ad 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.cpp +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/MiniAppManager.cpp @@ -1,91 +1,91 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "MiniAppManager.h" MiniAppManager* MiniAppManager::GetInstance() { static MiniAppManager instance; return &instance; } // Attention: Name of the miniApp must be the last argument!!! // it will be cut off from the rest of the arguments and then // the app will be run int MiniAppManager::RunMiniApp(int argc, char* argv[]) { try { std::string nameOfMiniApp; std::map< std::string, MiniAppFunction >::iterator it = m_Functions.begin(); if( argc < 2) { std::cout << "Please choose the mini app to execute: " << std::endl; for(int i=0; it != m_Functions.end(); ++i,++it) { std::cout << "(" << i << ")" << " " << it->first << std::endl; } std::cout << "Please select: "; int choose; std::cin >> choose; it = m_Functions.begin(); std::advance(it, choose); if( it != m_Functions.end() ) nameOfMiniApp = it->first; } else { - nameOfMiniApp = argv[argc-1]; - --argc; + nameOfMiniApp = argv[1]; + //--argc; } it = m_Functions.find(nameOfMiniApp); if(it == m_Functions.end()) { std::ostringstream s; s << "MiniApp (" << nameOfMiniApp << ") not found!"; throw std::invalid_argument(s.str().c_str()); } MITK_INFO << "Start " << nameOfMiniApp << " .."; MiniAppFunction func = it->second; return func( argc, argv ); } catch(std::exception& e) { MITK_ERROR << e.what(); } catch(...) { MITK_ERROR << "Unknown error occurred"; } return EXIT_FAILURE; } ///////////////////// // MiniAppFunction // ///////////////////// MiniAppManager::MiniAppFunction MiniAppManager::AddFunction(const std::string& name, MiniAppFunction func) { m_Functions.insert( std::pair(name, func) ); return func; } diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp new file mode 100755 index 0000000000..f95d8dfe8d --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp @@ -0,0 +1,174 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "ctkCommandLineParser.h" + +int StreamlineTracking(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", ctkCommandLineParser::String, "input tensor image (.dti)", mitk::Any(), false); + parser.addArgument("seed", "s", ctkCommandLineParser::String, "binary seed image"); + parser.addArgument("mask", "m", ctkCommandLineParser::String, "binary mask image"); + parser.addArgument("minFA", "t", ctkCommandLineParser::Float, "minimum fractional anisotropy threshold", 0.15, true); + parser.addArgument("minCurv", "c", ctkCommandLineParser::Float, "minimum curvature radius in mm (default = 0.5*minimum-spacing)"); + parser.addArgument("stepSize", "s", ctkCommandLineParser::Float, "stepsize in mm (default = 0.1*minimum-spacing)"); + parser.addArgument("tendf", "f", ctkCommandLineParser::Float, "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0, true); + parser.addArgument("tendg", "g", ctkCommandLineParser::Float, "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0, true); + parser.addArgument("numSeeds", "n", ctkCommandLineParser::Int, "Number of seeds per voxel.", 1, true); + parser.addArgument("minLength", "l", ctkCommandLineParser::Float, "minimum fiber length in mm", 20, true); + + parser.addArgument("interpolate", "a", ctkCommandLineParser::Bool, "Use linear interpolation", false, true); + parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output fiber bundle (.fib)", mitk::Any(), false); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + string dtiFileName = mitk::any_cast(parsedArgs["input"]); + string outFileName = mitk::any_cast(parsedArgs["outFile"]); + + float minFA = 0.15; + float minCurv = -1; + float stepSize = -1; + float tendf = 1; + float tendg = 0; + float minLength = 20; + int numSeeds = 1; + bool interpolate = false; + + if (parsedArgs.count("minCurv")) + minCurv = mitk::any_cast(parsedArgs["minCurv"]); + if (parsedArgs.count("minFA")) + minFA = mitk::any_cast(parsedArgs["minFA"]); + if (parsedArgs.count("stepSize")) + stepSize = mitk::any_cast(parsedArgs["stepSize"]); + if (parsedArgs.count("tendf")) + tendf = mitk::any_cast(parsedArgs["tendf"]); + if (parsedArgs.count("tendg")) + tendg = mitk::any_cast(parsedArgs["tendg"]); + if (parsedArgs.count("minLength")) + minLength = mitk::any_cast(parsedArgs["minLength"]); + if (parsedArgs.count("numSeeds")) + numSeeds = mitk::any_cast(parsedArgs["numSeeds"]); + + + if (parsedArgs.count("interpolate")) + interpolate = mitk::any_cast(parsedArgs["interpolate"]); + + + + try + { + RegisterDiffusionCoreObjectFactory(); + RegisterFiberTrackingObjectFactory(); + + // load input image + const std::string s1="", s2=""; + std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( dtiFileName, s1, s2, false ); + + MITK_INFO << "Loading tensor image ..."; + typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; + mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.at(0).GetPointer()); + ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); + mitk::CastToItkImage(mitkTensorImage, itk_dti); + + MITK_INFO << "Loading seed image ..."; + typedef itk::Image< unsigned char, 3 > ItkUCharImageType; + mitk::Image::Pointer mitkSeedImage = NULL; + if (parsedArgs.count("seed")) + mitkSeedImage = mitk::IOUtil::LoadImage(mitk::any_cast(parsedArgs["seed"])); + + MITK_INFO << "Loading mask image ..."; + mitk::Image::Pointer mitkMaskImage = NULL; + if (parsedArgs.count("mask")) + mitkMaskImage = mitk::IOUtil::LoadImage(mitk::any_cast(parsedArgs["mask"])); + + // instantiate tracker + typedef itk::StreamlineTrackingFilter< float > FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(itk_dti); + filter->SetSeedsPerVoxel(numSeeds); + filter->SetFaThreshold(minFA); + filter->SetMinCurvatureRadius(minCurv); + filter->SetStepSize(stepSize); + filter->SetF(tendf); + filter->SetG(tendg); + filter->SetInterpolate(interpolate); + filter->SetMinTractLength(minLength); + + if (mitkSeedImage.IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(mitkSeedImage, mask); + filter->SetSeedImage(mask); + } + + if (mitkMaskImage.IsNotNull()) + { + ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + mitk::CastToItkImage(mitkMaskImage, mask); + filter->SetMaskImage(mask); + } + + filter->Update(); + + vtkSmartPointer fiberBundle = filter->GetFiberPolyData(); + if ( fiberBundle->GetNumberOfLines()==0 ) + { + MITK_INFO << "No fibers reconstructed. Check parametrization."; + return EXIT_FAILURE; + } + mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberBundle); + + mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); + for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) + { + if ( (*it)->CanWriteBaseDataType(fib.GetPointer()) ) { + (*it)->SetFileName( outFileName.c_str() ); + (*it)->DoWrite( fib.GetPointer() ); + } + } + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + MITK_INFO << "DONE"; + return EXIT_SUCCESS; +} +RegisterFiberTrackingMiniApp(StreamlineTracking); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp new file mode 100755 index 0000000000..00996996dd --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp @@ -0,0 +1,247 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "MiniAppManager.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include "ctkCommandLineParser.h" +#include "ctkCommandLineParser.cpp" +#include +#include +#include +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +int TractometerAngularErrorTool(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", ctkCommandLineParser::String, "input tractogram (.fib, vtk ascii file format)", mitk::Any(), false); + parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "reference direction images", mitk::Any(), false); + parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", mitk::Any(), false); + parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); + parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); + parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results"); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + ctkCommandLineParser::StringContainerType referenceImages = mitk::any_cast(parsedArgs["reference"]); + + string fibFile = mitk::any_cast(parsedArgs["input"]); + + string maskImage(""); + if (parsedArgs.count("mask")) + maskImage = mitk::any_cast(parsedArgs["mask"]); + + float angularThreshold = 25; + if (parsedArgs.count("athresh")) + angularThreshold = mitk::any_cast(parsedArgs["athresh"]); + + string outRoot = mitk::any_cast(parsedArgs["out"]); + + bool verbose = false; + if (parsedArgs.count("verbose")) + verbose = mitk::any_cast(parsedArgs["verbose"]); + + try + { + RegisterDiffusionCoreObjectFactory(); + RegisterFiberTrackingObjectFactory(); + + typedef itk::Image ItkUcharImgType; + typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; + typedef itk::VectorContainer< int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; + typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; + + // load fiber bundle + mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::LoadDataNode(fibFile)->GetData()); + + // load reference directions + ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); + for (int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); + typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; + CasterType::Pointer caster = CasterType::New(); + caster->SetInput(img); + caster->Update(); + ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); + referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); + } + catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } + } + + // load/create mask image + ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); + if (maskImage.compare("")==0) + { + ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); + itkMaskImage->SetSpacing( dirImg->GetSpacing() ); + itkMaskImage->SetOrigin( dirImg->GetOrigin() ); + itkMaskImage->SetDirection( dirImg->GetDirection() ); + itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); + itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); + itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); + itkMaskImage->Allocate(); + itkMaskImage->FillBuffer(1); + } + else + { + mitk::Image::Pointer mitkMaskImage = dynamic_cast(mitk::IOUtil::LoadDataNode(maskImage)->GetData()); + CastToItkImage(mitkMaskImage, itkMaskImage); + } + + + // extract directions from fiber bundle + itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); + fOdfFilter->SetFiberBundle(inputTractogram); + fOdfFilter->SetMaskImage(itkMaskImage); + fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); + fOdfFilter->SetNormalizeVectors(true); + fOdfFilter->SetUseWorkingCopy(false); + fOdfFilter->Update(); + ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); + + if (verbose) + { + // write vector field + mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle(); + mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); + for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) + { + if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) { + QString outfilename(outRoot.c_str()); + outfilename += "_VECTOR_FIELD.fib"; + (*it)->SetFileName( outfilename.toStdString().c_str() ); + (*it)->DoWrite( directions.GetPointer() ); + } + } + + // write direction images + for (int i=0; iSize(); i++) + { + itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); + typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + QString outfilename(outRoot.c_str()); + outfilename += "_DIRECTION_"; + outfilename += QString::number(i); + outfilename += ".nrrd"; + MITK_INFO << "writing " << outfilename.toStdString(); + writer->SetFileName(outfilename.toStdString().c_str()); + writer->SetInput(itkImg); + writer->Update(); + } + + // write num direction image + { + ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); + typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; + WriterType::Pointer writer = WriterType::New(); + QString outfilename(outRoot.c_str()); + outfilename += "_NUM_DIRECTIONS.nrrd"; + MITK_INFO << "writing " << outfilename.toStdString(); + writer->SetFileName(outfilename.toStdString().c_str()); + writer->SetInput(numDirImage); + writer->Update(); + } + } + + // evaluate directions + EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); + evaluationFilter->SetImageSet(directionImageContainer); + evaluationFilter->SetReferenceImageSet(referenceImageContainer); + evaluationFilter->SetMaskImage(itkMaskImage); + evaluationFilter->Update(); + + if (verbose) + { + EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); + typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + QString outfilename(outRoot.c_str()); + outfilename += "_ERROR_IMAGE.nrrd"; + writer->SetFileName(outfilename.toStdString().c_str()); + writer->SetInput(angularErrorImage); + writer->Update(); + } + + QString logFile(outRoot.c_str()); logFile += "_ANGULAR_ERROR.csv"; + QFile file(logFile); + file.open(QIODevice::WriteOnly | QIODevice::Text); + QTextStream out(&file); + QString sens = QString("Mean:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMeanAngularError()); + sens += ";"; + + sens += QString("Median:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMedianAngularError()); + sens += ";"; + + sens += QString("Maximum:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMaxAngularError()); + sens += ";"; + + sens += QString("Minimum:"); + sens += ","; + sens += QString::number(evaluationFilter->GetMinAngularError()); + sens += ";"; + + sens += QString("STDEV:"); + sens += ","; + sens += QString::number(std::sqrt(evaluationFilter->GetVarAngularError())); + sens += ";"; + + out << sens; + + MITK_INFO << "DONE"; + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} +RegisterFiberTrackingMiniApp(TractometerAngularErrorTool); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.cpp new file mode 100755 index 0000000000..72da974c5b --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.cpp @@ -0,0 +1,727 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ +/*========================================================================= + + Library: CTK + + Copyright (c) Kitware Inc. + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0.txt + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + +=========================================================================*/ + +// STL includes +#include + + +// CTK includes +#include "ctkCommandLineParser.h" + +using namespace std; + +namespace +{ +// -------------------------------------------------------------------------- +class CommandLineParserArgumentDescription +{ +public: + + + CommandLineParserArgumentDescription( + const string& longArg, const string& longArgPrefix, + const string& shortArg, const string& shortArgPrefix, + ctkCommandLineParser::Type type, const string& argHelp, + const mitk::Any& defaultValue, bool ignoreRest, + bool deprecated, bool optional) + : LongArg(longArg), LongArgPrefix(longArgPrefix), + ShortArg(shortArg), ShortArgPrefix(shortArgPrefix), + ArgHelp(argHelp), IgnoreRest(ignoreRest), NumberOfParametersToProcess(0), + Deprecated(deprecated), DefaultValue(defaultValue), Value(type), ValueType(type), Optional(optional) + { + Value = defaultValue; + + switch (type) + { + case ctkCommandLineParser::String: + { + NumberOfParametersToProcess = 1; + } + break; + case ctkCommandLineParser::Bool: + { + NumberOfParametersToProcess = 0; + } + break; + case ctkCommandLineParser::StringList: + { + NumberOfParametersToProcess = -1; + } + break; + case ctkCommandLineParser::Int: + { + NumberOfParametersToProcess = 1; + } + break; + case ctkCommandLineParser::Float: + { + NumberOfParametersToProcess = 1; + } + break; + default: + MITK_INFO << "Type not supported: " << static_cast(type); + } + + } + + ~CommandLineParserArgumentDescription(){} + + bool addParameter(const string& value); + + string helpText(); + + string LongArg; + string LongArgPrefix; + string ShortArg; + string ShortArgPrefix; + string ArgHelp; + bool IgnoreRest; + int NumberOfParametersToProcess; + bool Deprecated; + bool Optional; + + mitk::Any DefaultValue; + mitk::Any Value; + ctkCommandLineParser::Type ValueType; +}; + +// -------------------------------------------------------------------------- +bool CommandLineParserArgumentDescription::addParameter(const string &value) +{ + switch (ValueType) + { + case ctkCommandLineParser::String: + { + Value = value; + } + break; + case ctkCommandLineParser::Bool: + { + if (value.compare("true")==0) + Value = true; + else + Value = false; + } + break; + case ctkCommandLineParser::StringList: + { + try + { + ctkCommandLineParser::StringContainerType list = mitk::any_cast(Value); + list.push_back(value); + Value = list; + } + catch(...) + { + ctkCommandLineParser::StringContainerType list; + list.push_back(value); + Value = list; + } + } + break; + case ctkCommandLineParser::Int: + { + stringstream ss(value); + int i; + ss >> i; + Value = i; + } + break; + case ctkCommandLineParser::Float: + { + stringstream ss(value); + float f; + ss >> f; + Value = f; + } + break; + default: + return false; + } + + return true; +} + +// -------------------------------------------------------------------------- +string CommandLineParserArgumentDescription::helpText() +{ + string text; + + string shortAndLongArg; + if (!this->ShortArg.empty()) + { + shortAndLongArg = " "; + shortAndLongArg += this->ShortArgPrefix; + shortAndLongArg += this->ShortArg; + } + + if (!this->LongArg.empty()) + { + if (this->ShortArg.empty()) + shortAndLongArg.append(" "); + else + shortAndLongArg.append(", "); + + shortAndLongArg += this->LongArgPrefix; + shortAndLongArg += this->LongArg; + } + + text = text + shortAndLongArg + ", " + this->ArgHelp; + + if (this->Optional) + text += " (optional)"; + + if (!this->DefaultValue.Empty()) + { + text = text + ", (default: " + this->DefaultValue.ToString() + ")"; + } + text += "\n"; + return text; +} + +} + +// -------------------------------------------------------------------------- +// ctkCommandLineParser::ctkInternal class + +// -------------------------------------------------------------------------- +class ctkCommandLineParser::ctkInternal +{ +public: + ctkInternal() + : Debug(false), FieldWidth(0), StrictMode(false) + {} + + ~ctkInternal() { } + + CommandLineParserArgumentDescription* argumentDescription(const string& argument); + + vector ArgumentDescriptionList; + map ArgNameToArgumentDescriptionMap; + map > GroupToArgumentDescriptionListMap; + + StringContainerType UnparsedArguments; + StringContainerType ProcessedArguments; + string ErrorString; + bool Debug; + int FieldWidth; + string LongPrefix; + string ShortPrefix; + string CurrentGroup; + string DisableQSettingsLongArg; + string DisableQSettingsShortArg; + bool StrictMode; +}; + +// -------------------------------------------------------------------------- +// ctkCommandLineParser::ctkInternal methods + +// -------------------------------------------------------------------------- +CommandLineParserArgumentDescription* +ctkCommandLineParser::ctkInternal::argumentDescription(const string& argument) +{ + string unprefixedArg = argument; + + if (!LongPrefix.empty() && argument.compare(0, LongPrefix.size(), LongPrefix)==0) + { + // Case when (ShortPrefix + UnPrefixedArgument) matches LongPrefix + if (argument == LongPrefix && !ShortPrefix.empty() && argument.compare(0, ShortPrefix.size(), ShortPrefix)==0) + { + unprefixedArg = argument.substr(ShortPrefix.size(),argument.size()); + } + else + { + unprefixedArg = argument.substr(LongPrefix.size(),argument.size()); + } + } + else if (!ShortPrefix.empty() && argument.compare(0, ShortPrefix.size(), ShortPrefix)==0) + { + unprefixedArg = argument.substr(ShortPrefix.size(),argument.size()); + } + else if (!LongPrefix.empty() && !ShortPrefix.empty()) + { + return 0; + } + + if (ArgNameToArgumentDescriptionMap.count(unprefixedArg)) + { + return this->ArgNameToArgumentDescriptionMap[unprefixedArg]; + } + return 0; +} + +// -------------------------------------------------------------------------- +// ctkCommandLineParser methods + +// -------------------------------------------------------------------------- +ctkCommandLineParser::ctkCommandLineParser() +{ + this->Internal = new ctkInternal(); +} + +// -------------------------------------------------------------------------- +ctkCommandLineParser::~ctkCommandLineParser() +{ + delete this->Internal; +} + +// -------------------------------------------------------------------------- +map ctkCommandLineParser::parseArguments(const StringContainerType& arguments, + bool* ok) +{ + // Reset + this->Internal->UnparsedArguments.clear(); + this->Internal->ProcessedArguments.clear(); + this->Internal->ErrorString.clear(); + // foreach (CommandLineParserArgumentDescription* desc, this->Internal->ArgumentDescriptionList) + for (int i=0; iArgumentDescriptionList.size(); i++) + { + CommandLineParserArgumentDescription* desc = Internal->ArgumentDescriptionList.at(i); + desc->Value = mitk::Any(desc->ValueType); + if (!desc->DefaultValue.Empty()) + { + desc->Value = desc->DefaultValue; + } + } + + bool error = false; + bool ignoreRest = false; + CommandLineParserArgumentDescription * currentArgDesc = 0; + vector parsedArgDescriptions; + for(int i = 1; i < arguments.size(); ++i) + { + string argument = arguments.at(i); + if (this->Internal->Debug) { MITK_DEBUG << "Processing" << argument; } + + // should argument be ignored ? + if (ignoreRest) + { + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: IgnoreRest flag was been set"; + } + this->Internal->UnparsedArguments.push_back(argument); + continue; + } + + // Skip if the argument does not start with the defined prefix + if (!(argument.compare(0, Internal->LongPrefix.size(), Internal->LongPrefix)==0 + || argument.compare(0, Internal->ShortPrefix.size(), Internal->ShortPrefix)==0)) + { + if (this->Internal->StrictMode) + { + this->Internal->ErrorString = "Unknown argument "; + this->Internal->ErrorString += argument; + error = true; + break; + } + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: It does not start with the defined prefix"; + } + this->Internal->UnparsedArguments.push_back(argument); + continue; + } + + // Skip if argument has already been parsed ... + bool alreadyProcessed = false; + for (int i=0; iProcessedArguments.size(); i++) + if (argument.compare(Internal->ProcessedArguments.at(i))==0) + { + alreadyProcessed = true; + break; + } + + if (alreadyProcessed) + { + if (this->Internal->StrictMode) + { + this->Internal->ErrorString = "Argument "; + this->Internal->ErrorString += argument; + this->Internal->ErrorString += " already processed !"; + error = true; + break; + } + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: Already processed !"; + } + continue; + } + + // Retrieve corresponding argument description + currentArgDesc = this->Internal->argumentDescription(argument); + + // Is there a corresponding argument description ? + if (currentArgDesc) + { + // If the argument is deprecated, print the help text but continue processing + if (currentArgDesc->Deprecated) + { + MITK_WARN << "Deprecated argument " << argument << ": " << currentArgDesc->ArgHelp; + } + else + { + parsedArgDescriptions.push_back(currentArgDesc); + } + + this->Internal->ProcessedArguments.push_back(currentArgDesc->ShortArg); + this->Internal->ProcessedArguments.push_back(currentArgDesc->LongArg); + int numberOfParametersToProcess = currentArgDesc->NumberOfParametersToProcess; + ignoreRest = currentArgDesc->IgnoreRest; + if (this->Internal->Debug && ignoreRest) + { + MITK_DEBUG << " IgnoreRest flag is True"; + } + + // Is the number of parameters associated with the argument being processed known ? + if (numberOfParametersToProcess == 0) + { + currentArgDesc->addParameter("true"); + } + else if (numberOfParametersToProcess > 0) + { + string missingParameterError = + "Argument %1 has %2 value(s) associated whereas exacly %3 are expected."; + for(int j=1; j <= numberOfParametersToProcess; ++j) + { + if (i + j >= arguments.size()) + { +// this->Internal->ErrorString = +// missingParameterError.arg(argument).arg(j-1).arg(numberOfParametersToProcess); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + string parameter = arguments.at(i + j); + if (this->Internal->Debug) + { + MITK_DEBUG << " Processing parameter" << j << ", value:" << parameter; + } + if (this->argumentAdded(parameter)) + { +// this->Internal->ErrorString = +// missingParameterError.arg(argument).arg(j-1).arg(numberOfParametersToProcess); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + if (!currentArgDesc->addParameter(parameter)) + { +// this->Internal->ErrorString = string( +// "Value(s) associated with argument %1 are incorrect. %2"). +// arg(argument).arg(currentArgDesc->ExactMatchFailedMessage); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + } + // Update main loop increment + i = i + numberOfParametersToProcess; + } + else if (numberOfParametersToProcess == -1) + { + if (this->Internal->Debug) + { + MITK_DEBUG << " Proccessing StringList ..."; + } + int j = 1; + while(j + i < arguments.size()) + { + if (this->argumentAdded(arguments.at(j + i))) + { + if (this->Internal->Debug) + { + MITK_DEBUG << " No more parameter for" << argument; + } + break; + } + string parameter = arguments.at(j + i); + + if (parameter.compare(0, Internal->LongPrefix.size(), Internal->LongPrefix)==0 + || parameter.compare(0, Internal->ShortPrefix.size(), Internal->ShortPrefix)==0) + { + j--; + break; + } + + if (this->Internal->Debug) + { + MITK_DEBUG << " Processing parameter" << j << ", value:" << parameter; + } + if (!currentArgDesc->addParameter(parameter)) + { +// this->Internal->ErrorString = string( +// "Value(s) associated with argument %1 are incorrect. %2"). +// arg(argument).arg(currentArgDesc->ExactMatchFailedMessage); +// if (this->Internal->Debug) { MITK_DEBUG << this->Internal->ErrorString; } + if (ok) { *ok = false; } + return map(); + } + j++; + } + // Update main loop increment + i = i + j; + } + } + else + { + if (this->Internal->StrictMode) + { + this->Internal->ErrorString = "Unknown argument "; + this->Internal->ErrorString += argument; + error = true; + break; + } + if (this->Internal->Debug) + { + MITK_DEBUG << " Skipping: Unknown argument"; + } + this->Internal->UnparsedArguments.push_back(argument); + } + } + + if (ok) + { + *ok = !error; + } + + map parsedArguments; + + int obligatoryArgs = 0; + vector::iterator it; + for(it = Internal->ArgumentDescriptionList.begin(); it != Internal->ArgumentDescriptionList.end();++it) + { + CommandLineParserArgumentDescription* desc = *it; + + if(!desc->Optional) + obligatoryArgs++; + } + + int parsedObligatoryArgs = 0; + for(it = parsedArgDescriptions.begin(); it != parsedArgDescriptions.end();++it) + { + CommandLineParserArgumentDescription* desc = *it; + + string key; + if (!desc->LongArg.empty()) + { + key = desc->LongArg; + } + else + { + key = desc->ShortArg; + } + + if(!desc->Optional) + parsedObligatoryArgs++; + + std::pair elem; elem.first = key; elem.second = desc->Value; + parsedArguments.insert(elem); + } + + if (obligatoryArgs>parsedObligatoryArgs) + { + parsedArguments.clear(); + cout << helpText(); + } + + return parsedArguments; +} + +// ------------------------------------------------------------------------- +map ctkCommandLineParser::parseArguments(int argc, char** argv, bool* ok) +{ + StringContainerType arguments; + + // Create a StringContainerType of arguments + for(int i = 0; i < argc; ++i) + arguments.push_back(argv[i]); + + return this->parseArguments(arguments, ok); +} + +// ------------------------------------------------------------------------- +string ctkCommandLineParser::errorString() const +{ + return this->Internal->ErrorString; +} + +// ------------------------------------------------------------------------- +const ctkCommandLineParser::StringContainerType& ctkCommandLineParser::unparsedArguments() const +{ + return this->Internal->UnparsedArguments; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::addArgument(const string& longarg, const string& shortarg, + Type type, const string& argHelp, + const mitk::Any& defaultValue, bool optional, bool ignoreRest, + bool deprecated) +{ + if (longarg.empty() && shortarg.empty()) { return; } + + /* Make sure it's not already added */ + bool added = this->Internal->ArgNameToArgumentDescriptionMap.count(longarg); + if (added) { return; } + + added = this->Internal->ArgNameToArgumentDescriptionMap.count(shortarg); + if (added) { return; } + + CommandLineParserArgumentDescription* argDesc = + new CommandLineParserArgumentDescription(longarg, this->Internal->LongPrefix, + shortarg, this->Internal->ShortPrefix, type, + argHelp, defaultValue, ignoreRest, deprecated, optional); + + int argWidth = 0; + if (!longarg.empty()) + { + this->Internal->ArgNameToArgumentDescriptionMap[longarg] = argDesc; + argWidth += longarg.size() + this->Internal->LongPrefix.size(); + } + if (!shortarg.empty()) + { + this->Internal->ArgNameToArgumentDescriptionMap[shortarg] = argDesc; + argWidth += shortarg.size() + this->Internal->ShortPrefix.size() + 2; + } + argWidth += 5; + + // Set the field width for the arguments + if (argWidth > this->Internal->FieldWidth) + { + this->Internal->FieldWidth = argWidth; + } + + this->Internal->ArgumentDescriptionList.push_back(argDesc); + this->Internal->GroupToArgumentDescriptionListMap[this->Internal->CurrentGroup].push_back(argDesc); +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::addDeprecatedArgument( + const string& longarg, const string& shortarg, const string& argHelp) +{ + addArgument(longarg, shortarg, StringList, argHelp, mitk::Any(), false, true); +} + +// -------------------------------------------------------------------------- +int ctkCommandLineParser::fieldWidth() const +{ + return this->Internal->FieldWidth; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::beginGroup(const string& description) +{ + this->Internal->CurrentGroup = description; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::endGroup() +{ + this->Internal->CurrentGroup.clear(); +} + +// -------------------------------------------------------------------------- +string ctkCommandLineParser::helpText() const +{ + string text; + vector deprecatedArgs; + + // Loop over grouped argument descriptions + map >::iterator it; + for(it = Internal->GroupToArgumentDescriptionListMap.begin(); it != Internal->GroupToArgumentDescriptionListMap.end();++it) + { + if (!(*it).first.empty()) + { + text = text + "\n" + (*it).first + "\n"; + } + + vector::iterator it2; + for(it2 = (*it).second.begin(); it2 != (*it).second.end(); ++it2) + { + CommandLineParserArgumentDescription* argDesc = *it2; + if (argDesc->Deprecated) + { + deprecatedArgs.push_back(argDesc); + } + else + { + text += argDesc->helpText(); + } + } + } + + if (!deprecatedArgs.empty()) + { + text += "\nDeprecated arguments:\n"; + vector::iterator it2; + for(it2 = deprecatedArgs.begin(); it2 != deprecatedArgs.end(); ++it2) + { + CommandLineParserArgumentDescription* argDesc = *it2; + text += argDesc->helpText(); + } + } + + return text; +} + +// -------------------------------------------------------------------------- +bool ctkCommandLineParser::argumentAdded(const string& argument) const +{ + return this->Internal->ArgNameToArgumentDescriptionMap.count(argument); +} + +// -------------------------------------------------------------------------- +bool ctkCommandLineParser::argumentParsed(const string& argument) const +{ + for (int i=0; iProcessedArguments.size(); i++) + if (argument.compare(Internal->ProcessedArguments.at(i))==0) + return true; + return false; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::setArgumentPrefix(const string& longPrefix, const string& shortPrefix) +{ + this->Internal->LongPrefix = longPrefix; + this->Internal->ShortPrefix = shortPrefix; +} + +// -------------------------------------------------------------------------- +void ctkCommandLineParser::setStrictModeEnabled(bool strictMode) +{ + this->Internal->StrictMode = strictMode; +} + diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.h b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.h new file mode 100755 index 0000000000..e7eb221379 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/ctkCommandLineParser.h @@ -0,0 +1,422 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ +/*========================================================================= + + Library: CTK + + Copyright (c) Kitware Inc. + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0.txt + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + +=========================================================================*/ + +#ifndef __ctkCommandLineParser_h +#define __ctkCommandLineParser_h + +#include +#include + + +/** + * \ingroup Core + * + * The CTK command line parser. + * + * Use this class to add information about the command line arguments + * your program understands and to easily parse them from a given list + * of strings. + * + * This parser provides the following features: + * + *
    + *
  • Add arguments by supplying a long name and/or a short name. + * Arguments are validated using a regular expression. They can have + * a default value and a help string.
    • + *
  • Deprecated arguments.
    • + *
  • Custom regular expressions for argument validation.
    • + *
  • Set different argument name prefixes for native platform look and feel.
    • + *
  • QSettings support. Default values for arguments can be read from + * a QSettings object.
    • + *
  • Create a help text for the command line arguments with support for + * grouping arguments.
    • + *
+ * + * Here is an example how to use this class inside a main function: + * + * \code + * #include + * #include + * #include + * + * int main(int argc, char** argv) + * { + * QCoreApplication app(argc, argv); + * // This is used by QSettings + * QCoreApplication::setOrganizationName("MyOrg"); + * QCoreApplication::setApplicationName("MyApp"); + * + * ctkCommandLineParser parser; + * // Use Unix-style argument names + * parser.setArgumentPrefix("--", "-"); + * // Enable QSettings support + * parser.enableSettings("disable-settings"); + * + * // Add command line argument names + * parser.addArgument("disable-settings", "", mitk::Any::Bool, "Do not use QSettings"); + * parser.addArgument("help", "h", mitk::Any::Bool, "Show this help text"); + * parser.addArgument("search-paths", "s", mitk::Any::StringList, "A list of paths to search"); + * + * // Parse the command line arguments + * bool ok = false; + * map parsedArgs = parser.parseArguments(QCoreApplication::arguments(), &ok); + * if (!ok) + * { + * QTextStream(stderr, QIODevice::WriteOnly) << "Error parsing arguments: " + * << parser.errorString() << "\n"; + * return EXIT_FAILURE; + * } + * + * // Show a help message + * if (parsedArgs.contains("help") || parsedArgs.contains("h")) + * { + * QTextStream(stdout, QIODevice::WriteOnly) << parser.helpText(); + * return EXIT_SUCCESS; + * } + * + * // Do something + * + * return EXIT_SUCCESS; + * } + * \endcode + */ + +using namespace std; + +class ctkCommandLineParser +{ + +public: + + enum Type { + String = 0, + Bool = 1, + StringList = 2, + Int = 3, + Float = 4 + }; + + typedef std::vector< std::string > StringContainerType; + + /** + * Constructs a parser instance. + * + * If QSettings support is enabled by a call to enableSettings() + * a default constructed QSettings instance will be used when parsing + * the command line arguments. Make sure to call QCoreApplication::setOrganizationName() + * and QCoreApplication::setApplicationName() before using default + * constructed QSettings objects. + * + * @param newParent The QObject parent. + */ + ctkCommandLineParser(); + + ~ctkCommandLineParser(); + + /** + * Parse a given list of command line arguments. + * + * This method parses a list of string elements considering the known arguments + * added by calls to addArgument(). If any one of the argument + * values does not match the corresponding regular expression, + * ok is set to false and an empty map object is returned. + * + * The keys in the returned map object correspond to the long argument string, + * if it is not empty. Otherwise, the short argument string is used as key. The + * mitk::Any values can safely be converted to the type specified in the + * addArgument() method call. + * + * @param arguments A StringContainerType containing command line arguments. Usually + * given by QCoreApplication::arguments(). + * @param ok A pointer to a boolean variable. Will be set to true + * if all regular expressions matched, false otherwise. + * @return A map object mapping the long argument (if empty, the short one) + * to a mitk::Any containing the value. + */ + + map parseArguments(const StringContainerType &arguments, bool* ok = 0); + + /** + * Convenient method allowing to parse a given list of command line arguments. + * @see parseArguments(const StringContainerType &, bool*) + */ + map parseArguments(int argc, char** argv, bool* ok = 0); + + /** + * Returns a detailed error description if a call to parseArguments() + * failed. + * + * @return The error description, empty if no error occured. + * @see parseArguments(const StringContainerType&, bool*) + */ + string errorString() const; + + /** + * This method returns all unparsed arguments, i.e. all arguments + * for which no long or short name has been registered via a call + * to addArgument(). + * + * @see addArgument() + * + * @return A list containing unparsed arguments. + */ + const StringContainerType& unparsedArguments() const; + + /** + * Checks if the given argument has been added via a call + * to addArgument(). + * + * @see addArgument() + * + * @param argument The argument to be checked. + * @return true if the argument was added, false + * otherwise. + */ + bool argumentAdded(const string& argument) const; + + /** + * Checks if the given argument has been parsed successfully by a previous + * call to parseArguments(). + * + * @param argument The argument to be checked. + * @return true if the argument was parsed, false + * otherwise. + */ + bool argumentParsed(const string& argument) const; + + /** + * Adds a command line argument. An argument can have a long name + * (like --long-argument-name), a short name (like -l), or both. The type + * of the argument can be specified by using the type parameter. + * The following types are supported: + * + * + * + * + * + * + * + * + *
Type# of parametersDefault regular exprExample
mitk::Any::String1.*--test-string StringParameter
mitk::Any::Bool0does not apply--enable-something
mitk::Any::StringList-1.*--test-list string1 string2
mitk::Any::Int1-?[0-9]+--test-int -5
+ * + * The regular expressions are used to validate the parameters of command line + * arguments. You can restrict the valid set of parameters by calling + * setExactMatchRegularExpression() for your argument. + * + * Optionally, a help string and a default value can be provided for the argument. If + * the mitk::Any type of the default value does not match type, an + * exception is thrown. Arguments with default values are always returned by + * parseArguments(). + * + * You can also declare an argument deprecated, by setting deprecated + * to true. Alternatively you can add a deprecated argument by calling + * addDeprecatedArgument(). + * + * If the long or short argument has already been added, or if both are empty strings, + * the method call has no effect. + * + * @param longarg The long argument name. + * @param shortarg The short argument name. + * @param type The argument type (see the list above for supported types). + * @param argHelp A help string describing the argument. + * @param defaultValue A default value for the argument. + * @param ignoreRest All arguments after the current one will be ignored. + * @param deprecated Declares the argument deprecated. + * + * @see setExactMatchRegularExpression() + * @see addDeprecatedArgument() + * @throws std::logic_error If the mitk::Any type of defaultValue + * does not match type, a std::logic_error is thrown. + */ + void addArgument(const string& longarg, const string& shortarg, + Type type, const string& argHelp = string(), + const mitk::Any& defaultValue = mitk::Any(), bool optional=true, + bool ignoreRest = false, bool deprecated = false); + + /** + * Adds a deprecated command line argument. If a deprecated argument is provided + * on the command line, argHelp is displayed in the console and + * processing continues with the next argument. + * + * Deprecated arguments are grouped separately at the end of the help text + * returned by helpText(). + * + * @param longarg The long argument name. + * @param shortarg The short argument name. + * @param argHelp A help string describing alternatives to the deprecated argument. + */ + void addDeprecatedArgument(const string& longarg, const string& shortarg, + const string& argHelp); + + /** + * Sets a custom regular expression for validating argument parameters. The method + * errorString() can be used the get the last error description. + * + * @param argument The previously added long or short argument name. + * @param expression A regular expression which the arugment parameters must match. + * @param exactMatchFailedMessage An error message explaining why the parameter did + * not match. + * + * @return true if the argument was found and the regular expression was set, + * false otherwise. + * + * @see errorString() + */ + bool setExactMatchRegularExpression(const string& argument, const string& expression, + const string& exactMatchFailedMessage); + + /** + * The field width for the argument names without the help text. + * + * @return The argument names field width in the help text. + */ + int fieldWidth() const; + + /** + * Creates a help text containing properly formatted argument names and help strings + * provided by calls to addArgument(). The arguments can be grouped by + * using beginGroup() and endGroup(). + * + * @param charPad The padding character. + * @return The formatted help text. + */ + string helpText() const; + + /** + * Sets the argument prefix for long and short argument names. This can be used + * to create native command line arguments without changing the calls to + * addArgument(). For example on Unix-based systems, long argument + * names start with "--" and short names with "-", while on Windows argument names + * always start with "/". + * + * Note that all methods in ctkCommandLineParser which take an argument name + * expect the name as it was supplied to addArgument. + * + * Example usage: + * + * \code + * ctkCommandLineParser parser; + * parser.setArgumentPrefix("--", "-"); + * parser.addArgument("long-argument", "l", mitk::Any::String); + * StringContainerType args; + * args << "program name" << "--long-argument Hi"; + * parser.parseArguments(args); + * \endcode + * + * @param longPrefix The prefix for long argument names. + * @param shortPrefix The prefix for short argument names. + */ + void setArgumentPrefix(const string& longPrefix, const string& shortPrefix); + + /** + * Begins a new group for documenting arguments. All newly added arguments via + * addArgument() will be put in the new group. You can close the + * current group by calling endGroup() or be opening a new group. + * + * Note that groups cannot be nested and all arguments which do not belong to + * a group will be listed at the top of the text created by helpText(). + * + * @param description The description of the group + */ + void beginGroup(const string& description); + + /** + * Ends the current group. + * + * @see beginGroup(const string&) + */ + void endGroup(); + + /** + * Enables QSettings support in ctkCommandLineParser. If an argument name is found + * in the QSettings instance with a valid mitk::Any, the value is considered as + * a default value and overwrites default values registered with + * addArgument(). User supplied values on the command line overwrite + * values in the QSettings instance, except for arguments with multiple parameters + * which are merged with QSettings values. Call mergeSettings(false) + * to disable merging. + * + * See ctkCommandLineParser(QSettings*) for information about how to + * supply a QSettings instance. + * + * Additionally, a long and short argument name can be specified which will disable + * QSettings support if supplied on the command line. The argument name must be + * registered as a regular argument via addArgument(). + * + * @param disableLongArg Long argument name. + * @param disableShortArg Short argument name. + * + * @see ctkCommandLineParser(QSettings*) + */ + void enableSettings(const string& disableLongArg = "", + const string& disableShortArg = ""); + + /** + * Controlls the merging behavior of user values and QSettings values. + * + * If merging is on (the default), user supplied values for an argument + * which can take more than one parameter are merged with values stored + * in the QSettings instance. If merging is off, the user values overwrite + * the QSettings values. + * + * @param merge true enables QSettings merging, false + * disables it. + */ + void mergeSettings(bool merge); + + /** + * Can be used to check if QSettings support has been enabled by a call to + * enableSettings(). + * + * @return true if QSettings support is enabled, false + * otherwise. + */ + bool settingsEnabled() const; + + + /** + * Can be used to teach the parser to stop parsing the arguments and return False when + * an unknown argument is encountered. By default StrictMode is disabled. + * + * @see parseArguments(const StringContainerType &, bool*) + */ + void setStrictModeEnabled(bool strictMode); + +private: + class ctkInternal; + ctkInternal * Internal; +}; + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 5f6ac6b882..912c88e510 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,100 +1,105 @@ set(CPP_FILES -# DataStructures -> FiberBundleX + MiniApps/ctkCommandLineParser.h + + # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures IODataStructures/mitkFiberTrackingObjectFactory.cpp # Rendering Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp #Rendering/mitkPlanarFigureMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp ) set(H_FILES + MiniApps/ctkCommandLineParser.h + # Rendering Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h #Rendering/mitkPlanarFigureMapper3D.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h IODataStructures/mitkFiberTrackingObjectFactory.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h Algorithms/itkFibersFromPlanarFiguresFilter.h Algorithms/itkTractsToDWIImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkKspaceImageFilter.h Algorithms/itkDftImageFilter.h Algorithms/itkAddArtifactsToDwiImageFilter.h Algorithms/itkFieldmapGeneratorFilter.h + Algorithms/itkEvaluateDirectionImagesFilter.h # (old) Tractography Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/itkStreamlineTrackingFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h # Signal Models SignalModels/mitkDiffusionSignalModel.h SignalModels/mitkTensorModel.h SignalModels/mitkBallModel.h SignalModels/mitkDotModel.h SignalModels/mitkAstroStickModel.h SignalModels/mitkStickModel.h SignalModels/mitkDiffusionNoiseModel.h SignalModels/mitkRicianNoiseModel.h SignalModels/mitkKspaceArtifact.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin # Shaders Shaders/mitkShaderFiberClipping.xml ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui index 418098a551..ceca723203 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkZeppelinModelParametersWidgetControls.ui @@ -1,124 +1,124 @@ QmitkZeppelinModelParametersWidgetControls 0 0 410 116 0 0 Form 0 Fractional anisotropy of resulting tensor. <html><head/><body><p>-</p></body></html> Diffusivity along second and third eigenvector. 5 1.000000000000000 0.000100000000000 0.000250000000000 T2 relaxation time of this compartment (in milliseconds). 10000 - 100 + 200 Diffusivity along largest eigenvector. 5 1.000000000000000 0.000100000000000 0.001000000000000 <html><head/><body><p><span style=" font-style:italic;">d</span><span style=" vertical-align:sub;">⟂</span>:</p></body></html> <html><head/><body><p><span style=" font-style:italic;">d</span><span style=" vertical-align:sub;">||</span>:</p></body></html> <html><head/><body><p><span style=" font-style:italic;">T2</span> - Relaxation:</p></body></html> <html><head/><body><p><span style=" font-style:italic;">FA</span>:</p></body></html> m_T2box m_D1box m_D2box