diff --git a/Algorithms/CMakeLists.txt b/Algorithms/CMakeLists.txt index 8c124b3e7aaf9b8ea8bff6b664bae1969cf7c4e0..e7ff595d12817268afd0b84972f3986ec5c3ceba 100644 --- a/Algorithms/CMakeLists.txt +++ b/Algorithms/CMakeLists.txt @@ -12,14 +12,13 @@ add_library(ttkAlgorithms ${TTK_BUILD_TYPE} itkFiberBundleStatisticsCalculator.cxx ) -target_link_libraries(ttkAlgorithms - ITKCommon - vtkCommon - vtkIO +target_link_libraries(ttkAlgorithms + ITKCommon + ${VTK_LIBRARIES} ) if(NOT ${PROJECT_NAME}_INSTALL_NO_LIBRARIES) - install(TARGETS ttkAlgorithms + install(TARGETS ttkAlgorithms RUNTIME DESTINATION bin LIBRARY DESTINATION lib ARCHIVE DESTINATION lib diff --git a/Algorithms/itkITKTensorsToVTKTensorsFilter.txx b/Algorithms/itkITKTensorsToVTKTensorsFilter.txx index 5202b861c44c971ce1a50247858506af6822be6a..97fc08bf13d93cae7ee10800fc8e3786fef6210d 100755 --- a/Algorithms/itkITKTensorsToVTKTensorsFilter.txx +++ b/Algorithms/itkITKTensorsToVTKTensorsFilter.txx @@ -55,21 +55,21 @@ namespace itk ITKTensorsToVTKTensorsFilter< TTensorImage > ::~ITKTensorsToVTKTensorsFilter() { - m_VTKTensors->Delete(); + m_VTKTensors->Delete(); m_DirectionMatrix->Delete(); } /** - * Set a vtkStructuredPoints as input + * Set a vtkStructuredPoints as input */ template < class TTensorImage > void ITKTensorsToVTKTensorsFilter< TTensorImage > ::CopyVTKTensors( vtkStructuredPoints* p ) { - p->DeepCopy(m_VTKTensors); + p->DeepCopy(m_VTKTensors); } @@ -98,33 +98,33 @@ namespace itk ::GenerateData() { - TensorImageConstPointer input = this->GetInput(0); - if ( ! input ) - { + TensorImageConstPointer input = this->GetInput(0); + if ( ! input ) + { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: No input ITK tensor image has been set."); - } + } - // VTK only supports 3D image with 3x3 tensors - typename TensorImageType::IndexType idx = {{0,0,0}}; - TensorType t = input->GetPixel(idx); - if( input->GetImageDimension() != 3 || t.GetTensorDimension() != 3 ) - { + // VTK only supports 3D image with 3x3 tensors + typename TensorImageType::IndexType idx = {{0,0,0}}; + TensorType t = input->GetPixel(idx); + if( input->GetImageDimension() != 3 || t.GetTensorDimension() != 3 ) + { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: VTK only supports 3D images and 3x3 tensors."); - } + } - typename TensorImageType::SizeType size = input->GetLargestPossibleRegion().GetSize(); - int numVoxels = 1; - int dims[3]; - double origin[3]; - double spacing[3]; + typename TensorImageType::SizeType size = input->GetLargestPossibleRegion().GetSize(); + int numVoxels = 1; + int dims[3]; + double origin[3]; + double spacing[3]; - for( unsigned int i = 0; i < 3; i++ ) - { + for( unsigned int i = 0; i < 3; i++ ) + { numVoxels *= size[i]; dims[i] = size[i]; origin[i] = input->GetOrigin()[i]; spacing[i] = input->GetSpacing()[i]; - } + } vtkDataArray* data = 0; if (typeid(ScalarType)==typeid(double)) @@ -134,19 +134,19 @@ namespace itk else itkExceptionMacro (<<"data is not float or double, cannot convert"); - data->SetNumberOfComponents(9); - data->SetNumberOfTuples(numVoxels); + data->SetNumberOfComponents(9); + data->SetNumberOfTuples(numVoxels); - typedef ImageRegionConstIterator<TensorImageType> IteratorType; - IteratorType it(input, input->GetLargestPossibleRegion()); + typedef ImageRegionConstIterator<TensorImageType> IteratorType; + IteratorType it(input, input->GetLargestPossibleRegion()); - unsigned long step = numVoxels/100; - unsigned int ind = 0; + unsigned long step = numVoxels/100; + unsigned int ind = 0; - this->UpdateProgress ( 0.0 ); + this->UpdateProgress ( 0.0 ); - while( !it.IsAtEnd() ) - { + while( !it.IsAtEnd() ) + { TensorType tensor = it.Get(); double buffer[9]; buffer[0] = tensor.GetNthComponent(0); @@ -167,7 +167,7 @@ namespace itk { this->UpdateProgress ( double(ind)/double(numVoxels) ); } - } + } typename TTensorImage::DirectionType directions = input->GetDirection(); typename TTensorImage::PointType i_origin = input->GetOrigin(); @@ -182,16 +182,16 @@ namespace itk m_DirectionMatrix->MultiplyPoint (v_origin, v_origin2); for (int i=0; i<3; i++) m_DirectionMatrix->SetElement (i, 3, v_origin[i]-v_origin2[i]); - - this->UpdateProgress ( 1.0 ); - - m_VTKTensors->Initialize(); - m_VTKTensors->SetDimensions(dims); - m_VTKTensors->SetSpacing(spacing); - m_VTKTensors->SetOrigin(origin); - m_VTKTensors->GetPointData()->SetTensors(data); - data->Delete(); - m_VTKTensors->Update(); + + this->UpdateProgress ( 1.0 ); + + m_VTKTensors->Initialize(); + m_VTKTensors->SetDimensions(dims); + m_VTKTensors->SetSpacing(spacing); + m_VTKTensors->SetOrigin(origin); + m_VTKTensors->GetPointData()->SetTensors(data); + data->Delete(); +// m_VTKTensors->Update(); } diff --git a/Algorithms/itkRemoveNonPositiveTensorsTensorImageFilter.txx b/Algorithms/itkRemoveNonPositiveTensorsTensorImageFilter.txx index d37a160441b58136d80e6c819becd3b905da6212..e7f2c50d8dcc58a1d7eb10b1fc5c551962bb227b 100644 --- a/Algorithms/itkRemoveNonPositiveTensorsTensorImageFilter.txx +++ b/Algorithms/itkRemoveNonPositiveTensorsTensorImageFilter.txx @@ -163,13 +163,13 @@ namespace itk InputPixelType Tn = bit.GetNext (i); InputPixelType Tmn = bit.GetPrevious (i); - if( Tn.IsPositive() && Tn.IsFinite() ) + if( !Tn.IsZero() && Tn.IsPositive() && Tn.IsFinite() ) { Mean += Tn.Log() * input->GetSpacing()[i]; sum += input->GetSpacing()[i]; } - if( Tmn.IsPositive() && Tmn.IsFinite() ) + if( !Tmn.IsZero() && Tmn.IsPositive() && Tmn.IsFinite() ) { Mean += Tmn.Log() * input->GetSpacing()[i]; sum += input->GetSpacing()[i]; @@ -182,7 +182,10 @@ namespace itk else { Mean /= sum; - out = Mean.Exp(); + if ( !Mean.IsZero() && Mean.IsPositive() && Mean.IsFinite() ) + out = Mean.Exp(); + else + out = static_cast<ScalarType>( 0.0 ); } break; diff --git a/CMakeLists.txt b/CMakeLists.txt index d8a0a089e18c6a27fc260f7124d18867575dfbc3..d2a3bba6d33ae136ff1aeced735e8fe6ee3d9c40 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,13 +1,9 @@ -cmake_minimum_required(VERSION 2.4) -IF(COMMAND cmake_policy) - cmake_policy(VERSION 2.4) - cmake_policy(SET CMP0005 OLD) - cmake_policy(SET CMP0003 NEW) -ENDIF(COMMAND cmake_policy) - - -project(TTK) +cmake_minimum_required(VERSION 3.0) +set(CMAKE_MACOSX_RPATH OFF) +set(TTK_VERSION_NUMBER 1.4.0) +mark_as_advanced(TTK_VERSION_NUMBER) +project(TTK VERSION ${TTK_VERSION_NUMBER}) # On Visual Studio 8 MS deprecated C. This removes all 1.276E1265 security # warnings @@ -138,6 +134,8 @@ SET(ITKIO_LIBRARIES ITKIOBioRad ITKIOStimulate ITKIOMRC + ITKIOTIFF + ${ITKIOPhilipsREC_LIBRARIES} ) set(ITK_TRANSFORM_LIBRARIES @@ -221,7 +219,7 @@ endif(TTK_USE_VTKINRIA3D) # Header file used to pass CMake settings to the source code # ----------------------------------------------------------------------------- configure_file(${PROJECT_SOURCE_DIR}/ttkConfigure.h.in - ${PROJECT_BINARY_DIR}/ttkConfigure.h CONFIGURE IMMEDIATE) + ${PROJECT_BINARY_DIR}/ttkConfigure.h) if(NOT ${PROJECT_NAME}_INSTALL_NO_DEVELOPMENT) install(FILES ${TTK_BINARY_DIR}/ttkConfigure.h @@ -267,8 +265,8 @@ CONFIGURE_FILE(${TTK_SOURCE_DIR}/UseTTK.cmake.in # Save the compiler settings so another project can import them. SET(TTK_BUILD_SETTINGS_FILE ${TTK_BINARY_DIR}/TTKBuildSettings.cmake) -INCLUDE(${CMAKE_ROOT}/Modules/CMakeExportBuildSettings.cmake) -CMAKE_EXPORT_BUILD_SETTINGS(${TTK_BUILD_SETTINGS_FILE}) +#INCLUDE(${CMAKE_ROOT}/Modules/CMakeExportBuildSettings.cmake) +#CMAKE_EXPORT_BUILD_SETTINGS(${TTK_BUILD_SETTINGS_FILE}) # Create the TTKConfig.cmake file containing the TTK configuration. INCLUDE (${TTK_SOURCE_DIR}/ttkGenerateTTKConfig.cmake) diff --git a/Commands/CMakeLists.txt b/Commands/CMakeLists.txt index 863c6641b94170e01c834d3fd1378afed62a61b1..cca60b5b0eeabb06f4d51bf670a0319ccbca2743 100644 --- a/Commands/CMakeLists.txt +++ b/Commands/CMakeLists.txt @@ -75,10 +75,10 @@ ITKProgramFactory ITKTensor ttkAlgorithms - ${ITKIO_LIBRARIES} - ITKStatistics + ${ITKIO_LIBRARIES} + ITKStatistics ${ITK_TRANSFORM_LIBRARIES} - vtkIO + ${VTK_LIBRARIES} ${TTK_LAPACK_LIBRARIES} ) @@ -87,7 +87,7 @@ ) target_link_libraries(ttk - ttkCommands + ttkCommands ITKProgramFactory ${ITKIO_LIBRARIES} ) @@ -129,18 +129,18 @@ endif (${bigobjFile_COMPILE_FLAGS} STREQUAL NOTFOUND) endforeach(bigobjFile) endif(WIN32) - + target_link_libraries(ttkConvertCommands ${TTK_MIPS_LIBRARIES} ${TTK_VTKINRIA3D_LIBRARIES} ITKProgramFactory ITKTensor - ${ITKIO_LIBRARIES} + ${ITKIO_LIBRARIES} ITKSpatialObjects - vtkIO + ${VTK_LIBRARIES} ${TTK_LAPACK_LIBRARIES} ) - + add_executable(ttk-convert ttk-convert.cxx ) @@ -150,8 +150,8 @@ ITKProgramFactory ${ITKIO_LIBRARIES} ) - - + + @@ -193,7 +193,7 @@ ITKTensor ${ITKIO_LIBRARIES} ${ITK_TRANSFORM_LIBRARIES} - vtkIO + ${VTK_LIBRARIES} ${TTK_LAPACK_LIBRARIES} ) @@ -212,12 +212,12 @@ # A vector field converter add_executable(itkStanleyToITK itkStanleyVectorFieldToITKVectorField.cxx) - TARGET_LINK_LIBRARIES(itkStanleyToITK - ITKCommon - ${ITKIO_LIBRARIES} + TARGET_LINK_LIBRARIES(itkStanleyToITK + ITKCommon + ${ITKIO_LIBRARIES} ) - + # registration commands if (TTK_USE_GMM) @@ -236,31 +236,28 @@ target_link_libraries(AssessRegistrationTensor ${TTK_MIPS_LIBRARIES} ${TTK_VTKINRIA3D_LIBRARIES} - ${TTK_LAPACK_LIBRARIES} + ${TTK_LAPACK_LIBRARIES} ${ITKIO_LIBRARIES} ITKTensor - vtkCommon - vtkIO + ${VTK_LIBRARIES} ) target_link_libraries(DemonsRegistrationTensor ${TTK_MIPS_LIBRARIES} ${TTK_VTKINRIA3D_LIBRARIES} ${TTK_LAPACK_LIBRARIES} - ${ITKIO_LIBRARIES} + ${ITKIO_LIBRARIES} ${ITK_TRANSFORM_LIBRARIES} ITKTensor - vtkCommon - vtkIO + ${VTK_LIBRARIES} ) target_link_libraries(LogDomainDemonsRegistrationTensor ${TTK_MIPS_LIBRARIES} ${TTK_VTKINRIA3D_LIBRARIES} ${TTK_LAPACK_LIBRARIES} - ${ITKIO_LIBRARIES} + ${ITKIO_LIBRARIES} ${ITK_TRANSFORM_LIBRARIES} ITKTensor - vtkCommon - vtkIO + ${VTK_LIBRARIES} ) endif (TTK_USE_GMM) @@ -316,7 +313,7 @@ if(NOT ${PROJECT_NAME}_INSTALL_NO_RUNTIME) ) endif(NOT ${PROJECT_NAME}_INSTALL_NO_RUNTIME) - + if(NOT ${PROJECT_NAME}_INSTALL_NO_DEVELOPMENT) file(GLOB __files1 "${CMAKE_CURRENT_SOURCE_DIR}/*.h") diff --git a/Commands/itkConsolidateFiberBundleCommand.cxx b/Commands/itkConsolidateFiberBundleCommand.cxx index 9c5f6aaceb460e64e39c3c3155f01d765249c6f1..a3791c139f2a928e18d8a75cba4008f55cf43cac 100644 --- a/Commands/itkConsolidateFiberBundleCommand.cxx +++ b/Commands/itkConsolidateFiberBundleCommand.cxx @@ -30,7 +30,7 @@ namespace itk { - + ConsolidateFiberBundleCommand::ConsolidateFiberBundleCommand() { m_ShortDescription = "Consolidate a fiber bundle by retrieving point data from all fibers"; @@ -46,31 +46,31 @@ namespace itk int ConsolidateFiberBundleCommand::Execute (int narg, const char* arg[]) { - + GetPot cl(narg, const_cast<char**>(arg)); // argument parser if( cl.size() == 1 || cl.search(2, "--help", "-h") ) { std::cout << this->GetLongDescription() << std::endl; return -1; } - + const char* bundle = cl.follow("NoFile",2,"-i","-I"); const char* all_fibers = cl.follow("NoFile",2,"-a","-A"); const char* output = cl.follow("output.vtk",2,"-o","-O"); - + const bool IsInputPresent = cl.search(2,"-a","-A"); - + vtkPolyDataReader* reader2 = vtkPolyDataReader::New(); reader2->SetFileName (bundle); reader2->Update(); vtkPolyData* vtkBundle = reader2->GetOutput(); - + vtkPolyDataReader* reader1 = vtkPolyDataReader::New(); vtkPolyData* vtkAllFibers = NULL; vtkPoints* allPoints = NULL; vtkCellArray* lines = NULL; - + if (IsInputPresent) { reader1->SetFileName (all_fibers); @@ -81,65 +81,65 @@ namespace itk { vtkAllFibers = reader2->GetOutput(); } - + allPoints = vtkAllFibers->GetPoints(); - + lines = vtkBundle->GetLines(); lines->InitTraversal(); - + vtkPolyData* vtkOutput = vtkPolyData::New(); vtkOutput->Initialize(); - vtkOutput->Allocate(); - + vtkOutput->Allocate(); + vtkPoints* bundlePoints = vtkPoints::New(); bundlePoints->Initialize(); - + vtkUnsignedCharArray* colors = vtkUnsignedCharArray::New(); colors->SetNumberOfComponents(3); - + vtkIdType npt, *pto; vtkIdType cellId = lines->GetNextCell (npt, pto); - + while (cellId!=0) { - + vtkIdType* newLine = new vtkIdType[npt]; - + for( int i=0; i<npt; i++) { - double* pt = NULL; - pt = allPoints->GetPoint (pto[i]); - - newLine[i] = bundlePoints->InsertNextPoint (pt); - colors->InsertNextTuple( vtkAllFibers->GetPointData()->GetScalars()->GetTuple (pto[i])); + double* pt = NULL; + pt = allPoints->GetPoint (pto[i]); + + newLine[i] = bundlePoints->InsertNextPoint (pt); + colors->InsertNextTuple( vtkAllFibers->GetPointData()->GetScalars()->GetTuple (pto[i])); } - + vtkOutput->InsertNextCell (VTK_POLY_LINE, npt, newLine); - + // delete [] newLine; - + cellId = lines->GetNextCell (npt, pto); } - + vtkOutput->SetPoints (bundlePoints); vtkOutput->GetPointData()->SetScalars (colors); - - + + vtkPolyDataWriter* writer = vtkPolyDataWriter::New(); writer->SetFileName (output); - writer->SetInput (vtkOutput); + writer->SetInputData(vtkOutput); writer->Update(); - + reader1->Delete(); reader2->Delete(); vtkOutput->Delete(); bundlePoints->Delete(); writer->Delete(); colors->Delete(); - - + + return 0; - + } - + } diff --git a/Commands/itkRBFRotationMatrixExtrapolationCommand.cxx b/Commands/itkRBFRotationMatrixExtrapolationCommand.cxx index 1b9ac4e145a5bceafc63e9a1a036a3df471b90a4..4b09b10dd9bf26e3e8d3e73346a8070a5403b5dc 100644 --- a/Commands/itkRBFRotationMatrixExtrapolationCommand.cxx +++ b/Commands/itkRBFRotationMatrixExtrapolationCommand.cxx @@ -51,49 +51,49 @@ namespace itk RBFRotationMatrixExtrapolationCommand::~RBFRotationMatrixExtrapolationCommand() {} - + int RBFRotationMatrixExtrapolationCommand::Execute(int narg, const char* arg[]) { - + GetPot cl(narg, const_cast<char**>(arg)); // argument parser if( cl.size() == 1 || cl.search(2, "--help", "-h") ) { std::cout << this->GetLongDescription() << std::endl; return -1; } - + const bool IsInputPresent = cl.search(2,"-i","-I"); - const bool AreSourcesPresent = cl.search(1,"-rot"); + const bool AreSourcesPresent = cl.search(1,"-rot"); const bool IsOutputPresent = cl.search(2,"-o","-O"); - + if(!IsInputPresent || !AreSourcesPresent || !IsOutputPresent) { std::cerr << "Input and (or) mask and (or) output not set." << std::endl; return -1; } - + const char* filemask = cl.follow("NoFile", 2, "-i","-I"); const char* output = cl.follow("NoFile", 2, "-o","-O"); const char* matrixFile = cl.follow("NoFile",1,"-rot"); - + if( strcmp(matrixFile,"NoFile")==0 || strcmp(output,"NoFile")==0 - || strcmp(filemask,"NoFile")==0) + || strcmp(filemask,"NoFile")==0) { std::cerr << "Input and (or) mask and (or) output not set." << std::endl; return -1; } - + const double sigma = cl.follow(1.0, 2, "-s","-S"); const double gamma = cl.follow(1.0, 2, "-g","-G"); const int t = cl.follow (1, 2, "-t", "-T"); - + std::cout << "Processing RBF interpolation of file: " << matrixFile << ", with: " << std::endl; std::cout << "Sigma: " << sigma << std::endl; std::cout << "Gamma: " << gamma << std::endl; std::cout << "Mask: " << filemask << std::endl; std::cout << "Output: " << output << std::endl; std::cout << std::flush; - + // typedefs typedef double ScalarType; typedef itk::Matrix<ScalarType, 3, 3> MatrixType; @@ -106,7 +106,7 @@ namespace itk typedef FilterType::VectorOfPointsType VectorOfPointsType; typedef itk::ImageFileReader<ImageType> ImageFileReaderType; - + // read the input image ImageFileReaderType::Pointer imReader = ImageFileReaderType::New(); imReader->SetFileName(filemask); @@ -121,50 +121,50 @@ namespace itk return -1; } std::cout << " Done." << std::endl; - + ImageType::Pointer input = imReader->GetOutput(); - - + + // read the model std::cout << "Reading: " << matrixFile; vtkUnstructuredGridReader* matReader = vtkUnstructuredGridReader::New(); matReader->SetFileName(matrixFile); vtkUnstructuredGrid* matrices = matReader->GetOutput(); - matrices->Update(); + matReader->Update(); std::cout << " Done." << std::endl; - - + + if( !matrices->GetPointData()->GetArray ("Frenet") ) { std::cout << "Error: Array Frenet was not found." << std::endl; return -1; } - - + + std::cout << "Converting..."; // convert the model to a vector of matrix + points VectorOfMatrixType vecM; VectorOfPointsType vecP; int numPoints = matrices->GetNumberOfPoints(); - + vtkDoubleArray* matArray = vtkDoubleArray::SafeDownCast (matrices->GetPointData()->GetArray("Frenet")); if( !matArray ) { std::cout << "Error: Cannot cast array." << std::endl; return -1; } - + for(int i=0;i<numPoints;i++) { double pt[3]; matrices->GetPoint(i,pt); PointType p; for(int m=0;m<3;m++) - p[m]=pt[m]; - + p[m]=pt[m]; + double matCoefs[9]; matArray->GetTuple(i,matCoefs); - + MatrixType M; M (0,0) = matCoefs[0]; M (1,0) = matCoefs[1]; @@ -175,15 +175,15 @@ namespace itk M (0,2) = matCoefs[6]; M (1,2) = matCoefs[7]; M (2,2) = matCoefs[8]; - + vecM.push_back( M ); vecP.push_back( p ); - + } matReader->Delete(); std::cout << "Done." << std::endl; - - + + FilterType::Pointer myFilter = FilterType::New(); myFilter->SetInput(input); myFilter->SetRotationMatrices(vecM); @@ -191,8 +191,8 @@ namespace itk myFilter->SetSigma(sigma); myFilter->SetGamma(gamma); myFilter->SetNumberOfThreads (t); - - + + std::cout << "Filtering..." << std::flush; try { @@ -204,7 +204,7 @@ namespace itk return -1; } std::cout << "Done." << std::endl; - + /* VectorOfMatrixType test = FilterType::InterpolateValuesAt (vecM,vecP,vecP,1,1); for( unsigned int i=0;i<test.size();i++) @@ -212,36 +212,36 @@ namespace itk std::cout << "Pos: " << vecP[i] << std::endl; std::cout << "Real: " << vecM[i] << std::endl; std::cout << "Int: " << test[i] << std::endl; - + }*/ - - + + MatrixImageType::Pointer itkMatImage = myFilter->GetOutput(); //MatrixImageType::SizeType size = itkMatImage->GetLargestPossibleRegion().GetSize(); MatrixImageType::PointType origin = itkMatImage->GetOrigin(); MatrixImageType::SpacingType spacing = itkMatImage->GetSpacing(); - + //double origin2[3] = {origin[0],origin[1],origin[2]}; //double spacing2[3] = {spacing[0],spacing[1],spacing[2]}; - + //vtkStructuredPoints* matImage = vtkStructuredPoints::New(); vtkUnstructuredGrid* matImage = vtkUnstructuredGrid::New(); //matImage->SetDimensions (size[0], size[1], size[2]); //matImage->SetOrigin ( origin2 ); //matImage->SetSpacing ( spacing2 ); - + vtkDoubleArray* matImageArray = vtkDoubleArray::New(); matImageArray->SetName ("Frenet"); matImageArray->SetNumberOfComponents (9); - + vtkPoints* myPoints = vtkPoints::New(); - + typedef itk::ImageRegionConstIteratorWithIndex<MatrixImageType> IteratorType; IteratorType it (itkMatImage, itkMatImage->GetLargestPossibleRegion() ); while( !it.IsAtEnd() ) { - + MatrixType m = it.Get(); double coef[9]; coef[0] = m (0,0); @@ -253,33 +253,33 @@ namespace itk coef[6] = m (0,2); coef[7] = m (1,2); coef[8] = m (2,2); - + matImageArray->InsertNextTuple (coef); - + MatrixImageType::IndexType index = it.GetIndex(); myPoints->InsertNextPoint (origin[0]+index[0]*spacing[0], origin[1]+index[1]*spacing[1], - origin[2]+index[2]*spacing[2]); - + origin[2]+index[2]*spacing[2]); + ++it; } - + matImage->SetPoints (myPoints); matImage->GetPointData()->AddArray (matImageArray); - + //vtkStructuredPointsWriter* writer = vtkStructuredPointsWriter::New(); vtkUnstructuredGridWriter* writer = vtkUnstructuredGridWriter::New(); - writer->SetInput (matImage); + writer->SetInputData(matImage); writer->SetFileName (output); std::cout << "Writing..." << std::flush; - writer->Write(); + writer->Write(); std::cout << "Done." << std::endl; - + myPoints->Delete(); matImage->Delete(); matImageArray->Delete(); writer->Delete(); - + return 0; } - + } diff --git a/Commands/itkRBFTensorExtrapolationCommand.cxx b/Commands/itkRBFTensorExtrapolationCommand.cxx index 08798a1ef37e95ede1b3e57e968a2fb622e7e99e..f6169b2b8d67d476055fb04a74200ab28a8ce5d8 100644 --- a/Commands/itkRBFTensorExtrapolationCommand.cxx +++ b/Commands/itkRBFTensorExtrapolationCommand.cxx @@ -49,46 +49,46 @@ namespace itk RBFTensorExtrapolationCommand::~RBFTensorExtrapolationCommand() {} - + int RBFTensorExtrapolationCommand::Execute(int narg, const char* arg[]) { - + GetPot cl(narg, const_cast<char**>(arg)); // argument parser if( cl.size() == 1 || cl.search(2, "--help", "-h") ) { std::cout << this->GetLongDescription() << std::endl; return -1; } - + const bool IsInputPresent = cl.search(2,"-i","-I"); - const bool AreSourcesPresent = cl.search(1,"-tens"); + const bool AreSourcesPresent = cl.search(1,"-tens"); const bool IsOutputPresent = cl.search(2,"-o","-O"); const bool IsTangentPresent = cl.search(1,"-tan"); const bool isEucSet = cl.search(1,"-euc"); - + if(!IsInputPresent || !AreSourcesPresent || !IsOutputPresent) { std::cerr << "Input and (or) mask and (or) output not set." << std::endl; return -1; } - + const char* filemask = cl.follow("NoFile", 2, "-i","-I"); const char* output = cl.follow("NoFile", 2, "-o","-O"); const char* tensorsFile = cl.follow("NoFile",1,"-tens"); const char* tangentFile = cl.follow("NoFile",1,"-tan"); - + if( strcmp(tensorsFile,"NoFile")==0 || strcmp(output,"NoFile")==0 - || strcmp(filemask,"NoFile")==0) + || strcmp(filemask,"NoFile")==0) { std::cerr << "Input and (or) mask and (or) output not set." << std::endl; return -1; } - + const double sigma = cl.follow(1.0, 2, "-s","-S"); const double gamma = cl.follow(1.0, 2, "-g","-G"); const int t = cl.follow (1, 2, "-t", "-T"); - + std::cout << "Processing RBF interpolation of file: " << tensorsFile << ", with: " << std::endl; std::cout << "Sigma: " << sigma << std::endl; std::cout << "Gamma: " << gamma << std::endl; @@ -98,7 +98,7 @@ namespace itk std::cout << "Mask: " << filemask << std::endl; std::cout << "Output: " << output << std::endl; std::cout << std::flush; - + // typedefs typedef double ScalarType; typedef itk::TensorImageIO<ScalarType, 3, 3> IOType; @@ -112,8 +112,8 @@ namespace itk typedef FilterType::PointType PointType; typedef FilterType::VectorOfPointsType VectorOfPointsType; typedef itk::ImageFileReader<ImageType> ImageFileReaderType; - - + + // read the input image ImageFileReaderType::Pointer imReader = ImageFileReaderType::New(); @@ -129,35 +129,35 @@ namespace itk return -1; } std::cout << " Done." << std::endl; - + ImageType::Pointer input = imReader->GetOutput(); - - + + // read the model std::cout << "Reading: " << tensorsFile; vtkUnstructuredGridReader* tensReader = vtkUnstructuredGridReader::New(); tensReader->SetFileName(tensorsFile); vtkUnstructuredGrid* tensors = tensReader->GetOutput(); - tensors->Update(); + tensReader->Update(); std::cout << " Done." << std::endl; - + std::cout << "Converting..."; // convert the model to a vector of tensors + points VectorOfTensorsType vecT; VectorOfPointsType vecP; int numPoints = tensors->GetNumberOfPoints(); - + for(int i=0;i<numPoints;i++) { double pt[3]; tensors->GetPoint(i,pt); PointType p; for(int m=0;m<3;m++) - p[m]=pt[m]; - + p[m]=pt[m]; + double tensorCoefs[9]; tensors->GetPointData()->GetTensors()->GetTuple(i,tensorCoefs); - + TensorType T; T.SetNthComponent ( 0, tensorCoefs[0] ); T.SetNthComponent ( 1, tensorCoefs[3] ); @@ -165,7 +165,7 @@ namespace itk T.SetNthComponent ( 3, tensorCoefs[6] ); T.SetNthComponent ( 4, tensorCoefs[7] ); T.SetNthComponent ( 5, tensorCoefs[8] ); - + if( isEucSet ) { vecT.push_back(T); @@ -174,13 +174,13 @@ namespace itk { vecT.push_back(T.Log()); } - - + + vecP.push_back(p); } tensReader->Delete(); std::cout << "Done." << std::endl; - + VectorOfPointsType vecTangent; if(IsTangentPresent) { @@ -188,24 +188,24 @@ namespace itk vtkUnstructuredGridReader* tanReader = vtkUnstructuredGridReader::New(); tanReader->SetFileName(tangentFile); vtkUnstructuredGrid* tangent = tanReader->GetOutput(); - tangent->Update(); - + tanReader->Update(); + int nTan = tangent->GetNumberOfPoints(); - + for(int i=0;i<nTan;i++) { - double pt[3]; - tangent->GetPoint(i,pt); - PointType p; - for(int m=0;m<3;m++) + double pt[3]; + tangent->GetPoint(i,pt); + PointType p; + for(int m=0;m<3;m++) p[m]=pt[m]; - vecTangent.push_back(p); + vecTangent.push_back(p); } - + tangent->Delete(); std::cout << "Done." << std::endl; - } - + } + FilterType::Pointer myFilter = FilterType::New(); myFilter->SetInput(input); myFilter->SetTensors(vecT); @@ -214,7 +214,7 @@ namespace itk myFilter->SetSigma(sigma); myFilter->SetGamma(gamma); myFilter->SetNumberOfThreads (t); - + try { myFilter->Update(); @@ -224,29 +224,29 @@ namespace itk std::cerr << e << std::endl; return -1; } - - + + ExpFilterType::Pointer myExpFilter = 0; if( !isEucSet ) { myExpFilter = ExpFilterType::New(); myExpFilter->SetInput ( myFilter->GetOutput() ); - + std::cout << "Filtering..." << std::flush; try { - myExpFilter->Update(); + myExpFilter->Update(); } catch(itk::ExceptionObject &e) { - std::cerr << e << std::endl; - return -1; + std::cerr << e << std::endl; + return -1; } std::cout << "Done." << std::endl; } - - - + + + IOType::Pointer writer = IOType::New(); writer->SetFileName(output); if( isEucSet ) @@ -257,7 +257,7 @@ namespace itk { writer->SetInput( myExpFilter->GetOutput() ); } - + std::cout << "Writing..." << std::flush; try { @@ -269,8 +269,8 @@ namespace itk return -1; } std::cout << "Done." << std::endl; - + return 0; } - + } diff --git a/Commands/itkSparseTensorsExtrapolationCommand.cxx b/Commands/itkSparseTensorsExtrapolationCommand.cxx index b7ef200470936fba5dadcc8cbd178688de396bd7..4e535770ee25ea4ebfe51136f530d2d955b06e98 100644 --- a/Commands/itkSparseTensorsExtrapolationCommand.cxx +++ b/Commands/itkSparseTensorsExtrapolationCommand.cxx @@ -54,7 +54,7 @@ namespace itk int SparseTensorsExtrapolationCommand::Execute (int argc, const char *argv[]) { - + GetPot cl(argc, const_cast<char**>(argv)); // argument parser if( cl.size() == 1 || cl.search(2, "--help", "-h") ) { @@ -66,7 +66,7 @@ namespace itk const bool IsInputPresent = cl.search(2,"-i","-I"); const bool IsOutputPresent = cl.search(2,"-o","-O"); const bool AreSourcesPresent = cl.search(1,"-tens"); - + if(!IsInputPresent || !IsOutputPresent || !AreSourcesPresent) { std::cerr << "Error: Input and (or) output and (or) sources not set." << std::endl; @@ -77,8 +77,8 @@ namespace itk const char* output = cl.follow("NoFile", 2, "-o","-O"); const char* tensorsFile = cl.follow("NoFile",1,"-tens"); const bool tan = cl.follow(0,1,"-tan"); - - if(strcmp(tensorsFile,"NoFile")==0 || strcmp(output,"NoFile")==0 + + if(strcmp(tensorsFile,"NoFile")==0 || strcmp(output,"NoFile")==0 || strcmp(mask,"NoFile")==0) { std::cerr << "Error: Input and (or) output and (or) sources not set." << std::endl; @@ -91,33 +91,33 @@ namespace itk const double lambda = cl.follow(1.0,2,"-l","-L"); const double sigma = cl.follow(2.0,2,"-s","-S"); const double rms = cl.follow(1.0,1,"-rms"); - const int threads = cl.follow (1, 2, "-t","-T"); - - + const int threads = cl.follow (1, 2, "-t","-T"); + + std::cout << "Processing diffusion of file: " << tensorsFile << ", with: " << std::endl; std::cout << "Sigma = " << sigma << std::endl; std::cout << "Lambda = " << lambda << std::endl; - std::cout << "Delta t = " << deltat << std::endl; + std::cout << "Delta t = " << deltat << std::endl; std::cout << "Nite = " << nite << std::endl; - std::cout << "RMS = " << rms << std::endl; + std::cout << "RMS = " << rms << std::endl; std::cout << "Output: " << output << std::endl; std::cout << std::flush; - - + + // typedefs - typedef double ScalarType; + typedef double ScalarType; typedef itk::TensorImageIO<ScalarType, 3 ,3> IOType; typedef IOType::TensorImageType TensorImageType; typedef TensorImageType::PixelType TensorType; typedef itk::SparseTensorsDiffusionTensorImageFilter<TensorImageType, TensorImageType> - FilterType; + FilterType; typedef FilterType::VectorOfTensorsType VectorOfTensorsType; - typedef FilterType::PointType PointType; + typedef FilterType::PointType PointType; typedef FilterType::VectorOfPointsType VectorOfPointsType; - - + + // read the tensor field - IOType::Pointer tensReader = IOType::New(); + IOType::Pointer tensReader = IOType::New(); tensReader->SetFileName(mask); std::cout << "Reading..." << std::flush; try @@ -130,25 +130,25 @@ namespace itk return -1; } std::cout << "Done." << std::endl; - + TensorImageType::Pointer input = tensReader->GetOutput(); - - // read the model + + // read the model std::cout << "Reading: " << tensorsFile; vtkUnstructuredGridReader* reader = vtkUnstructuredGridReader::New(); reader->SetFileName(tensorsFile); vtkUnstructuredGrid* tensors = reader->GetOutput(); - tensors->Update(); + reader->Update(); std::cout << " Done." << std::endl; - - + + std::cout << "Converting..."; // convert the model to a vector of tensors + points VectorOfTensorsType vecT; VectorOfPointsType vecP; VectorOfPointsType vecTangent; int numPoints = tensors->GetNumberOfPoints(); - + for( int i=0; i<numPoints; i++) { // get the position @@ -156,12 +156,12 @@ namespace itk tensors->GetPoint(i,pt); PointType p; for(int m=0;m<3;m++) - p[m]=pt[m]; - - + p[m]=pt[m]; + + double tensorCoefs[9]; tensors->GetPointData()->GetTensors()->GetTuple(i,tensorCoefs); - + TensorType T; T.SetNthComponent ( 0, tensorCoefs[0] ); T.SetNthComponent ( 1, tensorCoefs[3] ); @@ -169,38 +169,38 @@ namespace itk T.SetNthComponent ( 3, tensorCoefs[6] ); T.SetNthComponent ( 4, tensorCoefs[7] ); T.SetNthComponent ( 5, tensorCoefs[8] ); - + if (T.GetTrace() > 0.01) { - vecT.push_back(T); - vecP.push_back(p); + vecT.push_back(T); + vecP.push_back(p); } if (tan) { - // get the tangent which is stored as the normal - double tangent[3]; - tensors->GetPointData()->GetNormals()->GetTuple (i, tangent); - PointType TAN; - for(int m=0;m<3;m++) - TAN[m]=tangent[m]; - vecTangent.push_back (TAN); - } - + // get the tangent which is stored as the normal + double tangent[3]; + tensors->GetPointData()->GetNormals()->GetTuple (i, tangent); + PointType TAN; + for(int m=0;m<3;m++) + TAN[m]=tangent[m]; + vecTangent.push_back (TAN); + } + } - reader->Delete(); + reader->Delete(); std::cout << "Done." << std::endl; - - + + // set up the pipeline typedef itk::LogTensorImageFilter<TensorImageType,TensorImageType> LogFilterType; LogFilterType::Pointer logFilter = LogFilterType::New(); logFilter->SetInput(input); logFilter->SetNumberOfThreads (threads); - + // set up the filter - FilterType::Pointer myFilter = FilterType::New(); + FilterType::Pointer myFilter = FilterType::New(); myFilter->SetTensors(vecT); myFilter->SetPoints(vecP); myFilter->SetTangents (vecTangent); @@ -210,17 +210,17 @@ namespace itk myFilter->SetTimeStep(deltat); myFilter->SetNumberOfIterations(nite); myFilter->SetRMSChange(rms); - myFilter->SetNumberOfThreads (threads); + myFilter->SetNumberOfThreads (threads); myFilter->SetInput(logFilter->GetOutput()); - - - + + + typedef itk::ExpTensorImageFilter<TensorImageType,TensorImageType> ExpFilterType; ExpFilterType::Pointer expFilter = ExpFilterType::New(); expFilter->SetInput(myFilter->GetOutput()); expFilter->SetNumberOfThreads (threads); - + std::cout << "Pipeline started."<< std::endl; try { @@ -231,10 +231,10 @@ namespace itk std::cerr << e; return -1; } - + std::cout << "Pipeline finished." << std::endl; - - + + tensReader->SetFileName(output); tensReader->SetInput(expFilter->GetOutput()); try @@ -245,10 +245,10 @@ namespace itk { std::cerr << e; return -1; - } - + } + return 0; - + } - + } diff --git a/Commands/itkTensorGaussianFilteringCommand.cxx b/Commands/itkTensorGaussianFilteringCommand.cxx index 8516d2ea754aa71851cfa95e63450e80d860aca7..8ff8bc31b23e5a399cd50b0de0da6f0777971012 100644 --- a/Commands/itkTensorGaussianFilteringCommand.cxx +++ b/Commands/itkTensorGaussianFilteringCommand.cxx @@ -44,141 +44,141 @@ namespace itk TensorGaussianFilteringCommand::~TensorGaussianFilteringCommand() {} - + int TensorGaussianFilteringCommand::Execute (int narg, const char* arg[]) { - // parse arguments - GetPot cl (narg, const_cast<char**>(arg)); - if( cl.size() == 1 || cl.search (2,"--help","-h") ) - { - std::cout << this->GetLongDescription() << std::endl; - return -1; - } - - const char* file_in = cl.follow("NoFile", 2, "-i","-I"); - const char* file_out = cl.follow("NoFile", 2, "-o","-O"); - const double sigma = cl.follow (1.0, 2, "-s", "-S"); - const bool euclidean = cl.search (1,"-euc"); - - typedef double ScalarType; - typedef itk::TensorImageIO<ScalarType, 3, 3> TensorIOType; - typedef TensorIOType::TensorImageType TensorImageType; - - - std::cout << "Reading: " << file_in; - std::cout << std::flush; - - // read in the tensor field - TensorIOType::Pointer myIO = TensorIOType::New(); - myIO->SetFileName (file_in); - try - { - myIO->Read(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - - std::cout << " Done." << std::endl; - - TensorImageType::Pointer myTensorImage = myIO->GetOutput(); - - - if( !euclidean ) - { - - std::cout << "Loging..." << std::flush; - - // log the tensor field - typedef itk::LogTensorImageFilter<TensorImageType, TensorImageType> LogFilterType; - LogFilterType::Pointer myLoger = LogFilterType::New(); - myLoger->SetInput (myTensorImage); - try - { - myLoger->Update(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - - std::cout << "Done." << std::endl; - - myTensorImage = myLoger->GetOutput(); - } - - - - std::cout << "Smoothing..." << std::flush; - typedef itk::DiscreteGaussianImageFilter<TensorImageType, TensorImageType> - GaussianFilterType; - GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New(); - gaussianFilter->SetInput ( myTensorImage ); - gaussianFilter->SetVariance ( sigma ); - gaussianFilter->SetUseImageSpacingOn (); - - try - { - gaussianFilter->Update(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - std::cout << "Done." << std::endl; - - myTensorImage = gaussianFilter->GetOutput(); - - - if( !euclidean ) - { - - std::cout << "Exping..." << std::flush; - - // log the tensor field - typedef itk::ExpTensorImageFilter<TensorImageType, TensorImageType> ExpFilterType; - ExpFilterType::Pointer myExper = ExpFilterType::New(); - myExper->SetInput ( myTensorImage ); - try - { - myExper->Update(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - - std::cout << "Done." << std::endl; - - myTensorImage = myExper->GetOutput(); - - } - - - // write the tensor field - TensorIOType::Pointer myWriter = TensorIOType::New(); - myWriter->SetFileName (file_out); - myWriter->SetInput ( myTensorImage ); - try - { - myWriter->Write(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - - std::cout << " Done." << std::endl; - +// // parse arguments +// GetPot cl (narg, const_cast<char**>(arg)); +// if( cl.size() == 1 || cl.search (2,"--help","-h") ) +// { +// std::cout << this->GetLongDescription() << std::endl; +// return -1; +// } + +// const char* file_in = cl.follow("NoFile", 2, "-i","-I"); +// const char* file_out = cl.follow("NoFile", 2, "-o","-O"); +// const double sigma = cl.follow (1.0, 2, "-s", "-S"); +// const bool euclidean = cl.search (1,"-euc"); + +// typedef double ScalarType; +// typedef itk::TensorImageIO<ScalarType, 3, 3> TensorIOType; +// typedef TensorIOType::TensorImageType TensorImageType; + + +// std::cout << "Reading: " << file_in; +// std::cout << std::flush; + +// // read in the tensor field +// TensorIOType::Pointer myIO = TensorIOType::New(); +// myIO->SetFileName (file_in); +// try +// { +// myIO->Read(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } + +// std::cout << " Done." << std::endl; + +// TensorImageType::Pointer myTensorImage = myIO->GetOutput(); + + +// if( !euclidean ) +// { + +// std::cout << "Loging..." << std::flush; + +// // log the tensor field +// typedef itk::LogTensorImageFilter<TensorImageType, TensorImageType> LogFilterType; +// LogFilterType::Pointer myLoger = LogFilterType::New(); +// myLoger->SetInput (myTensorImage); +// try +// { +// myLoger->Update(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } + +// std::cout << "Done." << std::endl; + +// myTensorImage = myLoger->GetOutput(); +// } + + + +// std::cout << "Smoothing..." << std::flush; +// typedef itk::DiscreteGaussianImageFilter<TensorImageType, TensorImageType> +// GaussianFilterType; +// GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New(); +// gaussianFilter->SetInput ( myTensorImage ); +// gaussianFilter->SetVariance ( sigma ); +// gaussianFilter->SetUseImageSpacingOn (); + +// try +// { +// gaussianFilter->Update(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } +// std::cout << "Done." << std::endl; + +// myTensorImage = gaussianFilter->GetOutput(); + + +// if( !euclidean ) +// { + +// std::cout << "Exping..." << std::flush; + +// // log the tensor field +// typedef itk::ExpTensorImageFilter<TensorImageType, TensorImageType> ExpFilterType; +// ExpFilterType::Pointer myExper = ExpFilterType::New(); +// myExper->SetInput ( myTensorImage ); +// try +// { +// myExper->Update(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } + +// std::cout << "Done." << std::endl; + +// myTensorImage = myExper->GetOutput(); + +// } + + +// // write the tensor field +// TensorIOType::Pointer myWriter = TensorIOType::New(); +// myWriter->SetFileName (file_out); +// myWriter->SetInput ( myTensorImage ); +// try +// { +// myWriter->Write(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } + +// std::cout << " Done." << std::endl; + return 0; - + } - + } diff --git a/Commands/itkTensorGradientMagnitudeCalculatorCommand.cxx b/Commands/itkTensorGradientMagnitudeCalculatorCommand.cxx index a718602c23324d15f3db7e3f3df3d3f3985e6d29..12c56f02a4e26e372dabfd0a756e8ed4040be7e5 100644 --- a/Commands/itkTensorGradientMagnitudeCalculatorCommand.cxx +++ b/Commands/itkTensorGradientMagnitudeCalculatorCommand.cxx @@ -42,144 +42,144 @@ namespace itk m_LongDescription += "-o [fileOut]\n\n"; m_LongDescription += m_ShortDescription; } - + TensorGradientMagnitudeCalculatorCommand::~TensorGradientMagnitudeCalculatorCommand() {} - + int TensorGradientMagnitudeCalculatorCommand::Execute (int narg, const char* arg[]) { - itk::Object::GlobalWarningDisplayOff(); - - // parse arguments - GetPot cl (narg, const_cast<char**>(arg)); - if( cl.size() == 1 || cl.search (2,"--help","-h") ) - { - std::cout << this->GetLongDescription() << std::endl; - return -1; - } - - const bool IsInputPresent = cl.search(2,"-i","-I"); - const bool IsOutputPresent = cl.search(2,"-o","-O"); - const bool euclidean = cl.search (1,"-euc"); - const double sigma = cl.follow (1.0, 2, "-s", "-S"); - - if(!IsInputPresent || !IsOutputPresent) - { - std::cerr << "Input file and/or output not set" << std::endl; - return -1; - } - - - const char* file_in = cl.follow("NoFile", 2, "-i","-I"); - const char* file_out = cl.follow("NoFile", 2, "-o","-O"); - if(strcmp(file_in,"NoFile")==0 || strcmp(file_out,"NoFile")==0) - { - std::cerr << "Input file and/or output not set" << std::endl; - return -1; - } - - - // Read the input: - typedef double ScalarType; - const unsigned int ImageDimension = 3; - const unsigned int TensorDimension = 3; - - typedef itk::TensorImageIO<ScalarType, TensorDimension, ImageDimension> IOType; - typedef IOType::TensorImageType TensorImageType; - - - TensorImageType::Pointer tensorImage = 0; - { - IOType::Pointer myReader = IOType::New(); - myReader->SetFileName (file_in); - - std::cout << "Reading: " << file_in; - std::cout << std::flush; - try - { - myReader->Read(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - std::cout << " Done." << std::endl; - - tensorImage = myReader->GetOutput(); - - } - - - /** Specific typedefs */ - typedef itk::LogTensorImageFilter<TensorImageType,TensorImageType> LogFilterType; - typedef itk::Image<ScalarType, ImageDimension> ScalarImageType; - - typedef itk::VectorGradientMagnitudeImageFilter<TensorImageType, ScalarType, ScalarImageType> GradientFilterType; - typedef itk::DiscreteGaussianImageFilter<TensorImageType, TensorImageType> GaussianFilterType; - - - if( !euclidean ) - { - LogFilterType::Pointer myLog = LogFilterType::New(); - myLog->SetInput ( tensorImage ); - - try - { - myLog->Update(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - - tensorImage = myLog->GetOutput(); - } - - - - GaussianFilterType::Pointer gaussian = GaussianFilterType::New(); - gaussian->SetInput ( tensorImage ); - gaussian->SetVariance ( sigma ); - - GradientFilterType::Pointer gradient = GradientFilterType::New(); - gradient->SetInput ( gaussian->GetOutput() ); - - try - { - gradient->Update(); - } - catch (itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - - - - // save the image: - typedef itk::ImageFileWriter<ScalarImageType> WriterType; - WriterType::Pointer myWriter = WriterType::New(); - myWriter->SetFileName(file_out); - myWriter->SetInput( gradient->GetOutput() ); - - std::cout << "Writing: " << file_out; - std::cout << std::flush; - try - { - myWriter->Update(); - } - catch(itk::ExceptionObject &e) - { - std::cerr << e; - return -1; - } - std::cout << " Done." << std::endl; - - - return 0; +// itk::Object::GlobalWarningDisplayOff(); + +// // parse arguments +// GetPot cl (narg, const_cast<char**>(arg)); +// if( cl.size() == 1 || cl.search (2,"--help","-h") ) +// { +// std::cout << this->GetLongDescription() << std::endl; +// return -1; +// } + +// const bool IsInputPresent = cl.search(2,"-i","-I"); +// const bool IsOutputPresent = cl.search(2,"-o","-O"); +// const bool euclidean = cl.search (1,"-euc"); +// const double sigma = cl.follow (1.0, 2, "-s", "-S"); + +// if(!IsInputPresent || !IsOutputPresent) +// { +// std::cerr << "Input file and/or output not set" << std::endl; +// return -1; +// } + + +// const char* file_in = cl.follow("NoFile", 2, "-i","-I"); +// const char* file_out = cl.follow("NoFile", 2, "-o","-O"); +// if(strcmp(file_in,"NoFile")==0 || strcmp(file_out,"NoFile")==0) +// { +// std::cerr << "Input file and/or output not set" << std::endl; +// return -1; +// } + + +// // Read the input: +// typedef double ScalarType; +// const unsigned int ImageDimension = 3; +// const unsigned int TensorDimension = 3; + +// typedef itk::TensorImageIO<ScalarType, TensorDimension, ImageDimension> IOType; +// typedef IOType::TensorImageType TensorImageType; + + +// TensorImageType::Pointer tensorImage = 0; +// { +// IOType::Pointer myReader = IOType::New(); +// myReader->SetFileName (file_in); + +// std::cout << "Reading: " << file_in; +// std::cout << std::flush; +// try +// { +// myReader->Read(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } +// std::cout << " Done." << std::endl; + +// tensorImage = myReader->GetOutput(); + +// } + + +// /** Specific typedefs */ +// typedef itk::LogTensorImageFilter<TensorImageType,TensorImageType> LogFilterType; +// typedef itk::Image<ScalarType, ImageDimension> ScalarImageType; + +// typedef itk::VectorGradientMagnitudeImageFilter<TensorImageType, ScalarType, ScalarImageType> GradientFilterType; +// typedef itk::DiscreteGaussianImageFilter<TensorImageType, TensorImageType> GaussianFilterType; + + +// if( !euclidean ) +// { +// LogFilterType::Pointer myLog = LogFilterType::New(); +// myLog->SetInput ( tensorImage ); + +// try +// { +// myLog->Update(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } + +// tensorImage = myLog->GetOutput(); +// } + + + +// GaussianFilterType::Pointer gaussian = GaussianFilterType::New(); +// gaussian->SetInput ( tensorImage ); +// gaussian->SetVariance ( sigma ); + +// GradientFilterType::Pointer gradient = GradientFilterType::New(); +// gradient->SetInput ( gaussian->GetOutput() ); + +// try +// { +// gradient->Update(); +// } +// catch (itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } + + + +// // save the image: +// typedef itk::ImageFileWriter<ScalarImageType> WriterType; +// WriterType::Pointer myWriter = WriterType::New(); +// myWriter->SetFileName(file_out); +// myWriter->SetInput( gradient->GetOutput() ); + +// std::cout << "Writing: " << file_out; +// std::cout << std::flush; +// try +// { +// myWriter->Update(); +// } +// catch(itk::ExceptionObject &e) +// { +// std::cerr << e; +// return -1; +// } +// std::cout << " Done." << std::endl; + + + return 0; } - + } diff --git a/Commands/itkTensorNormalizedGaussianInterpolationCommand.cxx b/Commands/itkTensorNormalizedGaussianInterpolationCommand.cxx index 13f76a4baa783e2207df8c95db3232fab9e733ae..b0ad57b3468dbbf3dd311f8c529d30ed77f3e2af 100644 --- a/Commands/itkTensorNormalizedGaussianInterpolationCommand.cxx +++ b/Commands/itkTensorNormalizedGaussianInterpolationCommand.cxx @@ -51,20 +51,20 @@ namespace itk int TensorNormalizedGaussianInterpolationCommand::Execute(int narg, const char* arg[]) { - + GetPot cl(narg, const_cast<char**>(arg)); // argument parser if( cl.size() == 1 || cl.search(2, "--help", "-h") ) { std::cout << this->GetLongDescription() << std::endl; return -1; } - + const bool IsInputPresent = cl.search(2,"-i","-I"); // == domain of diffusion const bool AreSourcesPresent = cl.search(1,"-tens"); const bool IsMapPresent = cl.search(1,"-map"); const bool IsOutputPresent = cl.search(1,"-o"); const bool IsTangentPresent = cl.search(1,"-tan"); - + if(!IsInputPresent || !AreSourcesPresent || !IsOutputPresent) { std::cerr << "Error: Input and (or) sources and (or) output not set." << std::endl; @@ -76,20 +76,20 @@ namespace itk const char* mapFile = cl.follow("NoFile", 1, "-map"); const char* output = cl.follow("NoFile", 2, "-o","-O"); const char* tangentFile = cl.follow("NoFile",1,"-tan"); - - + + if( strcmp(tensorsFile,"NoFile")==0 || strcmp(output,"NoFile")==0 - || strcmp(filemask,"NoFile")==0) + || strcmp(filemask,"NoFile")==0) { std::cerr << "Error: Input file and (or) mask and (or) output not set." << std::endl; return -1; } - + const double ox = cl.follow(0.0,1,"-or"); const double oy = cl.next(0.0); const double oz = cl.next(0.0); const double sigma = cl.follow(20.0,1,"-s"); - + std::cout << "Processing NG interpolation of file: " << tensorsFile << std::endl; std::cout << "Sigma: " << sigma << std::endl; if( IsTangentPresent ) @@ -98,8 +98,8 @@ namespace itk std::cout << "Mask: " << filemask << std::endl; std::cout << "Output: " << output << std::endl; std::cout << std::flush; - - + + // typedefs typedef double ScalarType; typedef itk::TensorImageIO<ScalarType, 3, 3> IOType; @@ -112,9 +112,9 @@ namespace itk typedef FilterType::PointType PointType; typedef FilterType::VectorOfPointsType VectorOfPointsType; typedef itk::ImageFileReader<ImageType> ImageFileReaderType; - - - + + + // read the input image ImageFileReaderType::Pointer imReader = ImageFileReaderType::New(); imReader->SetFileName(filemask); @@ -129,54 +129,54 @@ namespace itk return -1; } std::cout << " Done." << std::endl; - + ImageType::Pointer input = imReader->GetOutput(); - + // read the input image ImageFileReaderType::Pointer imReader2 = ImageFileReaderType::New(); if (IsMapPresent) { - + imReader2->SetFileName(mapFile); std::cout << "Reading: " << mapFile << std::flush; try { - imReader2->Update(); + imReader2->Update(); } catch(itk::ExceptionObject &e) { - std::cerr << e << std::endl; - return -1; + std::cerr << e << std::endl; + return -1; } std::cout << " Done." << std::endl; } ImageType::Pointer secondmap = imReader2->GetOutput(); - - // read the source tensors + + // read the source tensors std::cout << "Reading: " << tensorsFile; vtkUnstructuredGridReader* tensReader = vtkUnstructuredGridReader::New(); tensReader->SetFileName(tensorsFile); vtkUnstructuredGrid* tensors = tensReader->GetOutput(); - tensors->Update(); + tensReader->Update(); std::cout << " Done." << std::endl; - + std::cout << "Converting..."; // convert the model to a vector of tensors + points VectorOfTensorsType vecT; VectorOfPointsType vecP; int numPoints = tensors->GetNumberOfPoints(); - + for(int i=0; i<numPoints; i++) - { + { double pt[3]; tensors->GetPoint(i,pt); PointType p; for(int m=0;m<3;m++) - p[m]=pt[m]; - + p[m]=pt[m]; + double tensorCoefs[9]; tensors->GetPointData()->GetTensors()->GetTuple(i,tensorCoefs); - + TensorType T; T.SetNthComponent ( 0, tensorCoefs[0] ); T.SetNthComponent ( 1, tensorCoefs[3] ); @@ -184,18 +184,18 @@ namespace itk T.SetNthComponent ( 3, tensorCoefs[6] ); T.SetNthComponent ( 4, tensorCoefs[7] ); T.SetNthComponent ( 5, tensorCoefs[8] ); - + if (T.GetTrace() > 0.1) { - vecT.push_back(T); - vecP.push_back(p); + vecT.push_back(T); + vecP.push_back(p); } - + } tensReader->Delete(); std::cout << "Done." << std::endl; - - + + VectorOfPointsType vecTangent; if(IsTangentPresent) { @@ -203,25 +203,25 @@ namespace itk vtkUnstructuredGridReader* tanReader = vtkUnstructuredGridReader::New(); tanReader->SetFileName(tangentFile); vtkUnstructuredGrid* tangent = tanReader->GetOutput(); - tangent->Update(); - + tanReader->Update(); + int nTan = tangent->GetNumberOfPoints(); - + for(int i=0;i<nTan;i++) { - double pt[3]; - tangent->GetPoint(i,pt); - PointType p; - for(int m=0;m<3;m++) - p[m]=pt[m]; - vecTangent.push_back(p); + double pt[3]; + tangent->GetPoint(i,pt); + PointType p; + for(int m=0;m<3;m++) + p[m]=pt[m]; + vecTangent.push_back(p); } - + tangent->Delete(); } - - - FilterType::Pointer myFilter = FilterType::New(); + + + FilterType::Pointer myFilter = FilterType::New(); myFilter->SetTensors(vecT); myFilter->SetTangents(vecTangent); if (IsMapPresent) @@ -229,11 +229,11 @@ namespace itk myFilter->UseAuxiliaryMapOn(); myFilter->SetAuxiliaryMap (secondmap.GetPointer()); } - + myFilter->SetPoints(vecP); myFilter->SetSigma(sigma); myFilter->SetInput(input); - + std::cout << "Filtering..." << std::flush; try { @@ -245,13 +245,13 @@ namespace itk return -1; } std::cout << "Done." << std::endl; - - - + + + IOType::Pointer writer = IOType::New(); writer->SetFileName(output); writer->SetInput(myFilter->GetOutput()); - + std::cout << "Wrting..." << std::flush; try { @@ -263,8 +263,8 @@ namespace itk return -1; } std::cout << "Done." << std::endl; - + return 0; } - + } diff --git a/Commands/itkTensorsToVTKUnstructuredGridCommand.cxx b/Commands/itkTensorsToVTKUnstructuredGridCommand.cxx index a69976adaf3a32b04a99bd98734162eb8215ca67..224ec88ef8bc060a05ca049182369e17063841fd 100644 --- a/Commands/itkTensorsToVTKUnstructuredGridCommand.cxx +++ b/Commands/itkTensorsToVTKUnstructuredGridCommand.cxx @@ -43,13 +43,13 @@ namespace itk m_LongDescription += "-o [output vtk file (default : output.vtk)]\n\n"; m_LongDescription += m_ShortDescription; } - + TensorsToVTKUnstructuredGridCommand::~TensorsToVTKUnstructuredGridCommand() {} int TensorsToVTKUnstructuredGridCommand::Execute(int narg, const char* arg[]) { - + GetPot cl(narg, const_cast<char**>(arg)); // argument parser if( cl.size() == 1 || cl.search(2, "--help", "-h") ) { @@ -64,7 +64,7 @@ namespace itk std::cout << "inputfile: " << inputfile << std::endl; std::cout << "outputfile: " << outputfile << std::endl; std::cout << std::flush; - + // typedefs typedef double ScalarType; typedef itk::TensorImageIO<ScalarType, 3, 3> TensorIOType; @@ -74,20 +74,20 @@ namespace itk typedef TensorImageType::PointType PointType; typedef std::vector<TensorType> TensorArrayType; typedef std::vector<PointType> PointArrayType; - - // instantiation - + + // instantiation + // read the input tensors and put tham into a vtkUnstructuredGrid // they come from a text file listing all files to read, either vtk or itk... - + std::cout << "Reading input tensors: " << inputfile << std::endl; - + TensorArrayType vecT; PointArrayType vecP; itk::ContinuousIndex<double, 3> index; - + std::cout<<"reading list : "<<inputfile<<std::endl; - + std::ifstream inputliststream (inputfile); if(inputliststream.fail()) { @@ -97,10 +97,10 @@ namespace itk unsigned int NumberOfImageFiles = 0; inputliststream >> NumberOfImageFiles; std::cout<<"encountered : "<<NumberOfImageFiles<<std::endl; - + std::string sline = ""; itksys::SystemTools::GetLineFromStream(inputliststream, sline); - + for (unsigned int N=0; N<NumberOfImageFiles; N++) { std::string line = ""; @@ -109,90 +109,90 @@ namespace itk TensorIOType::Pointer reader = TensorIOType::New(); reader->SetFileName(line.c_str()); bool success = false; - + try { - reader->Read(); - success = true; + reader->Read(); + success = true; } catch(itk::ExceptionObject &) { - std::cerr << "cannot read with TensorIO"<< std::endl; + std::cerr << "cannot read with TensorIO"<< std::endl; } - - + + if (!success) { - vtkDataSetReader* vtkreader = vtkDataSetReader::New(); - vtkreader->SetFileName (line.c_str()); - vtkreader->Update(); - vtkPointSet* dataset = vtkPointSet::SafeDownCast (vtkreader->GetOutput()); - - if (!dataset || !dataset->GetNumberOfPoints() || !dataset->GetPointData()->GetTensors()) - { - std::cerr << "cannot read file "<<line.c_str()<< std::endl; - std::cerr << "skipping line..."<< std::endl; - vtkreader->Delete(); - continue; - } - - for (int i=0; i<dataset->GetNumberOfPoints(); i++) - { - PointType x; - for (unsigned int j=0; j<3; j++) - x[j] = dataset->GetPoint (i)[j]; - TensorType T; - double* vals = dataset->GetPointData()->GetTensors()->GetTuple (i); - - T[0] = vals[0]; - T[1] = vals[1]; - T[2] = vals[4]; - T[3] = vals[2]; - T[4] = vals[5]; - T[5] = vals[8]; - - vecT.push_back(T); - vecP.push_back(x); - } - - vtkreader->Delete(); - + vtkDataSetReader* vtkreader = vtkDataSetReader::New(); + vtkreader->SetFileName (line.c_str()); + vtkreader->Update(); + vtkPointSet* dataset = vtkPointSet::SafeDownCast (vtkreader->GetOutput()); + + if (!dataset || !dataset->GetNumberOfPoints() || !dataset->GetPointData()->GetTensors()) + { + std::cerr << "cannot read file "<<line.c_str()<< std::endl; + std::cerr << "skipping line..."<< std::endl; + vtkreader->Delete(); + continue; + } + + for (int i=0; i<dataset->GetNumberOfPoints(); i++) + { + PointType x; + for (unsigned int j=0; j<3; j++) + x[j] = dataset->GetPoint (i)[j]; + TensorType T; + double* vals = dataset->GetPointData()->GetTensors()->GetTuple (i); + + T[0] = vals[0]; + T[1] = vals[1]; + T[2] = vals[4]; + T[3] = vals[2]; + T[4] = vals[5]; + T[5] = vals[8]; + + vecT.push_back(T); + vecP.push_back(x); + } + + vtkreader->Delete(); + } - + TensorImageType::Pointer tensorimage = reader->GetOutput(); TensorIteratorType it(tensorimage, tensorimage->GetLargestPossibleRegion()); - + while( !it.IsAtEnd() ) { - PointType x; - tensorimage->TransformIndexToPhysicalPoint(it.GetIndex(), x); - - TensorType T = it.Get(); - - if ((T.GetTrace() > 0.0) && !T.HasNans()) - { - T = T.ApplyMatrix(tensorimage->GetDirection()); - vecT.push_back(T); - vecP.push_back(x); - } - ++it; + PointType x; + tensorimage->TransformIndexToPhysicalPoint(it.GetIndex(), x); + + TensorType T = it.Get(); + + if ((T.GetTrace() > 0.0) && !T.HasNans()) + { + T = T.ApplyMatrix(tensorimage->GetDirection()); + vecT.push_back(T); + vecP.push_back(x); + } + ++it; } - std::cout<<" done."<<std::endl; + std::cout<<" done."<<std::endl; } - + vtkUnstructuredGrid* crossvalidation = vtkUnstructuredGrid::New(); vtkDoubleArray* data = vtkDoubleArray::New(); vtkPoints* points = vtkPoints::New(); points->SetNumberOfPoints (vecP.size()); data->SetNumberOfComponents (9); data->SetNumberOfTuples (vecP.size()); - + for (unsigned int i=0; i<vecP.size(); i++) { double pt[3] = {vecP[i][0], - vecP[i][1], - vecP[i][2]}; - + vecP[i][1], + vecP[i][2]}; + points->SetPoint (i, pt); double vals[9]; vals[0] = vecT[i][0]; @@ -204,27 +204,27 @@ namespace itk vals[6] = vecT[i][3]; vals[7] = vecT[i][4]; vals[8] = vecT[i][5]; - + data->SetTuple (i, vals); } - + crossvalidation->SetPoints (points); crossvalidation->GetPointData()->SetTensors (data); - + vtkDataSetWriter* writer = vtkDataSetWriter::New(); - writer->SetFileName (outputfile); - writer->SetInput (crossvalidation); + writer->SetFileName(outputfile); + writer->SetInputData(crossvalidation); writer->Update(); - - + + writer->Delete(); points->Delete(); data->Delete(); crossvalidation->Delete(); - + std::cout<<" done."<<std::endl; - + return EXIT_SUCCESS; } - + } diff --git a/Commands/itkVTKFibersToITKGroupSpatialObjectCommand.cxx b/Commands/itkVTKFibersToITKGroupSpatialObjectCommand.cxx index 0f50a6ccf9132d1de3048ce7afc01e5ed4dadbfb..46d491ac9a9cfdd5a52efaabc41db660d129bb74 100644 --- a/Commands/itkVTKFibersToITKGroupSpatialObjectCommand.cxx +++ b/Commands/itkVTKFibersToITKGroupSpatialObjectCommand.cxx @@ -67,31 +67,30 @@ namespace itk const char* inrTensorFile = cl.follow ("NoFile",2,"-t","-T"); const double bvalue = cl.follow (1.0, 2, "-b","-B"); const char* itkFiberFile = cl.follow ("NoFile", 2, "-o","-O"); - + // Read the input: typedef double ScalarType; const unsigned int ImageDimension = 3; const unsigned int TensorDimension = 3; - + typedef itk::TensorImageIO<ScalarType, TensorDimension, ImageDimension> IOType; typedef IOType::TensorImageType TensorImageType; - typedef TensorImageType::PixelType TensorType; - + // read the vtk fiber file std::cout << "Reading: " << vtkFiberFile << std::endl; vtkPolyDataReader* reader = vtkPolyDataReader::New(); reader->SetFileName( vtkFiberFile ); reader->Update(); - + vtkPolyData* vtkFibers = reader->GetOutput(); - vtkFibers->Update(); - + reader->Update(); + // read the inrimage tensor file IOType::Pointer myReader = IOType::New(); myReader->SetFileName (inrTensorFile); - + std::cout << "Reading: " << inrTensorFile; std::cout << std::flush; try @@ -104,7 +103,7 @@ namespace itk return -1; } std::cout << " Done." << std::endl; - + // log the input typedef itk::LogTensorImageFilter<TensorImageType,TensorImageType> LogFilterType; @@ -131,7 +130,7 @@ namespace itk myConverter->SetInput (vtkFibers); myConverter->SetTensorImage (myTensors); myConverter->SetBVal (bvalue); - + std::cout << "Converting..."; std::cout << std::flush; try @@ -144,20 +143,20 @@ namespace itk return -1; } std::cout << "Done." << std::endl; - + GroupSpatialObjectType::Pointer output = myConverter->GetOutput(); - - + + std::cout << "Wrting..."; itk::SpatialObjectWriter<3>::Pointer fibWriter = itk::SpatialObjectWriter<3>::New(); fibWriter->SetInput(output); fibWriter->SetFileName(itkFiberFile); fibWriter->Update(); std::cout << "Done." << std::endl; - + reader->Delete(); - + return 0; } - + } diff --git a/Commands/itkWarpFibersCommand.cxx b/Commands/itkWarpFibersCommand.cxx index b4efbc6ab009f3c4261fba079f93a9d749532e64..f75c26078d5e8f522e7f9c389c83a7975161a3ba 100644 --- a/Commands/itkWarpFibersCommand.cxx +++ b/Commands/itkWarpFibersCommand.cxx @@ -51,10 +51,10 @@ namespace itk m_LongDescription += m_ShortDescription; } - + WarpFibersCommand::~WarpFibersCommand() {} - + int WarpFibersCommand::Execute (int narg, const char* arg[]) { @@ -66,32 +66,32 @@ namespace itk return -1; } - + const char* file_in = cl.follow ("",2,"-i","-I"); const char* file_out = cl.follow ("",2,"-o","-O"); const char* file_vector = cl.follow ("",2,"-v","-V"); const char* file_matrix = cl.follow ("",2,"-m","-M"); - - + + typedef itk::Vector<double, 3> VectorType; typedef itk::Image<VectorType, 3> VectorImageType; typedef itk::ImageFileReader<VectorImageType> ReaderType; typedef itk::VectorLinearInterpolateImageFunction<VectorImageType, double> InterpolateFunctionType; - - + + vtkPolyDataReader* reader = vtkPolyDataReader::New(); reader->SetFileName (file_in); - + vtkPolyData* bundle = reader->GetOutput(); - bundle->Update(); - + reader->Update(); + if ( strcmp (file_vector, "")!=0 && strcmp (file_matrix, "")!=0 ) { std::cout << "Error: please choose between vector OR matrix" << std::endl; return -1; } - - + + VectorImageType::Pointer vector = 0; if( strcmp (file_vector, "")!=0 ) { @@ -99,52 +99,52 @@ namespace itk reader_v->SetFileName (file_vector); try { - reader_v->Update(); + reader_v->Update(); } catch (itk::ExceptionObject &e) { - std::cerr <<e; - return -1; + std::cerr <<e; + return -1; } vector = reader_v->GetOutput(); - } - - - + } + + + InterpolateFunctionType::Pointer interpolator = InterpolateFunctionType::New(); if( !vector.IsNull() ) { interpolator->SetInputImage ( vector ); } - - + + typedef itk::MatrixOffsetTransformBase<double, 3, 3> LinearTransformType; typedef itk::AffineTransform<double, 3> AffineTransformType; - + LinearTransformType::Pointer transform = 0; if( strcmp (file_matrix, "")!=0 ) { itk::TransformFactory< LinearTransformType >::RegisterTransform (); itk::TransformFactory< AffineTransformType >::RegisterTransform (); - + typedef itk::TransformFileReader TransformReaderType; TransformReaderType::Pointer reader_t = TransformReaderType::New(); reader_t->SetFileName ( file_matrix ); try { - reader_t->Update(); + reader_t->Update(); } catch (itk::ExceptionObject &e) { - std::cerr << e; - return -1; + std::cerr << e; + return -1; } transform = dynamic_cast<LinearTransformType*>( reader_t->GetTransformList()->front().GetPointer() ); } - - + + vtkPoints* points = bundle->GetPoints(); - + /* vtkPolyData* new_bundle = vtkPolyData::New(); new_bundle->DeepCopy (bundle); @@ -152,11 +152,11 @@ namespace itk vtkUnsignedCharArray* new_colors = vtkUnsignedCharArray::New(); new_colors->SetNumberOfComponents( 4 ); new_colors->SetNumberOfTuples(new_points->GetNumberOfPoints()); - + vtkUnsignedCharArray* colors = vtkUnsignedCharArray::New(); colors->SetNumberOfComponents( 4 ); colors->SetNumberOfTuples(new_points->GetNumberOfPoints()); - + double n_color[4] = {0,0,255,255}; double color[4] = {255,0,0,255}; */ @@ -164,37 +164,37 @@ namespace itk { double pt[3]; points->GetPoint ( i, pt); - + VectorImageType::PointType pt_i; pt_i[0] = pt[0]; pt_i[1] = pt[1]; pt_i[2] = pt[2]; - - + + VectorImageType::PointType pt_n = 0.0; - + if( !vector.IsNull() ) { - VectorType vec = interpolator->Evaluate ( pt_i ); - pt_n = pt_i + vec; + VectorType vec = interpolator->Evaluate ( pt_i ); + pt_n = pt_i + vec; } else if (!transform.IsNull() ) { - pt_n = transform->TransformPoint (pt); + pt_n = transform->TransformPoint (pt); } pt[0] = pt_n[0]; pt[1] = pt_n[1]; pt[2] = pt_n[2]; - + points->SetPoint (i, pt); //new_colors->SetTuple (i, n_color); //colors->SetTuple (i, color); } - + //bundle->GetPointData()->SetScalars (colors); //new_bundle->GetPointData()->SetScalars ( new_colors ); - + /* vtkAppendPolyData* appender = vtkAppendPolyData::New(); appender->AddInput ( bundle ); @@ -208,14 +208,14 @@ namespace itk writer2->SetFileTypeToBinary(); writer2->Write(); */ - + //bundle->SetPoints (new_points); vtkPolyDataWriter* writer = vtkPolyDataWriter::New(); writer->SetFileName ( file_out ); - writer->SetInput ( bundle ); + writer->SetInputData( bundle ); writer->SetFileTypeToBinary(); writer->Write(); - + /* new_bundle->Delete(); new_colors->Delete(); @@ -223,11 +223,11 @@ namespace itk appender->Delete(); writer2->Delete(); */ - + reader->Delete(); writer->Delete(); - + return 0; } - + } diff --git a/Common/itkGradientFileReader.cxx b/Common/itkGradientFileReader.cxx index 42ff56289f1100f7a372074c0565b50c0e046f92..ca085934aef7a90971c69bb27084d2ca6983663a 100644 --- a/Common/itkGradientFileReader.cxx +++ b/Common/itkGradientFileReader.cxx @@ -20,7 +20,7 @@ #include "itkGradientFileReader.h" #include <itksys/SystemTools.hxx> -#include <itksys/ios/sstream> +//#include <itksys/ios/sstream> #include <sstream> #include <fstream> @@ -171,7 +171,7 @@ void GradientFileReader break; } - itksys_ios::istringstream parse ( line ); + std::istringstream parse ( line ); std::vector<ScalarType> values; @@ -230,7 +230,7 @@ void GradientFileReader continue; } - itksys_ios::istringstream parse ( line ); + std::istringstream parse ( line ); VectorType gradient; ScalarType value; @@ -285,7 +285,7 @@ void GradientFileReader // Push back itkDebugMacro ( "Name: \"" << Name << "\"" ); itkDebugMacro ( "Value: \"" << Value << "\"" ); - itksys_ios::istringstream parse ( Value ); + std::istringstream parse ( Value ); int vector_id = std::atoi (Name.c_str()); if (vector_id < 0 || vector_id > 512) diff --git a/Common/itkGradientFileWriter.h b/Common/itkGradientFileWriter.h index 4d522a72f64d936d8c599530cac9844a1e964f64..3f77541a0655bf4e2a6128d8bca229d694dfb3bf 100644 --- a/Common/itkGradientFileWriter.h +++ b/Common/itkGradientFileWriter.h @@ -22,7 +22,7 @@ #include <string> #include <vector> #include <itksys/SystemTools.hxx> -#include <itksys/ios/sstream> +//#include <itksys/ios/sstream> #include <sstream> #include <fstream> diff --git a/Common/itkLogTensorImageFilter.txx b/Common/itkLogTensorImageFilter.txx index a1d3ef9a58a747c66fe5d22fcf353f5385903adb..baf3c96f823d1d5e8c19ea58868fa6906941b6cb 100644 --- a/Common/itkLogTensorImageFilter.txx +++ b/Common/itkLogTensorImageFilter.txx @@ -49,8 +49,8 @@ namespace itk T = T.Log(); } catch( itk::ExceptionObject &) - { - throw itk::ExceptionObject(__FILE__,__LINE__,"Error in LogTensorImageFilter::ThreadedGenerateData()"); + { + throw itk::ExceptionObject(__FILE__,__LINE__,"Error in LogTensorImageFilter::ThreadedGenerateData()"); } } itOut.Set (T); diff --git a/Common/itkTensorImageIO.txx b/Common/itkTensorImageIO.txx index c7cea03a40bc99bc93a4934607731f42089c990c..69025971d92d2c2ef865153571d37be91e1c1c75 100644 --- a/Common/itkTensorImageIO.txx +++ b/Common/itkTensorImageIO.txx @@ -48,7 +48,6 @@ #include <vtkDoubleArray.h> #include <vtksys/SystemTools.hxx> -#include <vtksys/stl/string> namespace itk { @@ -76,12 +75,12 @@ namespace itk throw itk::ExceptionObject (__FILE__,__LINE__,"Error: FileName is not set."); } - vtksys_stl::string ext = vtksys::SystemTools::GetFilenameExtension (m_FileName.c_str()); + std::string ext = vtksys::SystemTools::GetFilenameExtension (m_FileName.c_str()); if( ext=="" ) { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: Extension is not set."); } - + try { if ( CheckExtension(m_FileName.c_str(), ".vtk") ) @@ -115,7 +114,7 @@ namespace itk std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Read()"); } - + } @@ -139,7 +138,7 @@ namespace itk { delete tens; throw itk::ExceptionObject(__FILE__,__LINE__,"Error while reading. "); - + } // get the important information about size, spacing, origin, etc. @@ -155,13 +154,13 @@ namespace itk RegionType region; region.SetSize (size); region.SetIndex (index); - + // Allocate the output image if(m_Output) { m_Output->Initialize(); } - + m_Output->SetRegions (region); m_Output->SetSpacing (spacing); m_Output->SetOrigin (origin); @@ -172,31 +171,31 @@ namespace itk // of the inrimage as the rotation instance is in voxel coordinates. // Note also that it SHOULD be a rigid rotation matrix, // in theory its det. should be 1... - + yav::Matrix4x4<double> matrix = tens->getTransformationMatrix(); DirectionType direction; - + for (unsigned int i=0; i<3; i++) for (unsigned int j=0; j<3; j++) - direction[i][j] = matrix[i][j]; - + direction[i][j] = matrix[i][j]; + m_Output->SetDirection(direction); try { - m_Output->Allocate(); + m_Output->Allocate(); } catch (itk::ExceptionObject &e) { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Read() during allocation."); - } + } // Copy the input to the output and free the input image tens->convertCast(yav::Inrimage::WT_DOUBLE_VECTOR); double *buffer = (double*)(tens->getData()); - + typedef ImageRegionIterator<TensorImageType> IteratorType; IteratorType it(m_Output, m_Output->GetLargestPossibleRegion()); unsigned long numPixels = m_Output->GetLargestPossibleRegion().GetNumberOfPixels(); @@ -204,21 +203,21 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while(!it.IsAtEnd()) { - unsigned int offset = ind*(TensorType::NDegreesOfFreedom); + unsigned int offset = ind*(TensorType::NDegreesOfFreedom); T* coeff = new T[TensorType::NDegreesOfFreedom]; for(unsigned int i=0;i<(TensorType::NDegreesOfFreedom);i++) coeff[i] = static_cast<T>(buffer[offset+i]); - + TensorType Tens(coeff); delete [] coeff; - + it.Set (Tens); - + ++it; ++ind; @@ -226,17 +225,17 @@ namespace itk { this->UpdateProgress ( double(ind)/double(numPixels) ); } - + } this->UpdateProgress ( 1.0 ); - - delete tens; + + delete tens; #else - std::cerr << "TTK was not compiled with inrimage support, sorry!" << std::endl; + std::cerr << "TTK was not compiled with inrimage support, sorry!" << std::endl; #endif } - + template<class T, unsigned int TensorDimension, unsigned int ImageDimension> void @@ -249,27 +248,27 @@ namespace itk { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: VTK only supports 3D images and 3x3 tensors."); } - - + + // actually read the file vtkStructuredPointsReader* reader = vtkStructuredPointsReader::New(); reader->SetFileName( filename ); - + if( !reader->IsFileStructuredPoints() ) { - reader->Delete(); + reader->Delete(); throw itk::ExceptionObject(__FILE__,__LINE__,"Error: File is not a VTK structured points."); } - + vtkStructuredPoints* tensors = reader->GetOutput(); - tensors->Update(); - + reader->Update(); + // Allocate the output image if( m_Output ) { m_Output->Initialize(); } - + typename TensorImageType::SizeType size; double spacing[3]; double origin[3]; @@ -279,17 +278,17 @@ namespace itk spacing[i] = tensors->GetSpacing()[i]; origin[i] = tensors->GetOrigin()[i]; } - + typename TensorImageType::RegionType region; typename TensorImageType::IndexType index = {{0,0,0}}; - + region.SetSize(size); region.SetIndex(index); - + m_Output->SetRegions(region); m_Output->SetSpacing(spacing); m_Output->SetOrigin(origin); - + try { m_Output->Allocate(); @@ -300,7 +299,7 @@ namespace itk reader->Delete(); throw itk::ExceptionObject(__FILE__,__LINE__,"Error during memory allocation"); } - + vtkDoubleArray* data = (vtkDoubleArray*)(tensors->GetPointData()->GetTensors()); @@ -312,13 +311,13 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while( !it.IsAtEnd() ) { - + double coefs[9]; data->GetTuple(ind,coefs); - + T newCoefs[6]; newCoefs[0] = static_cast<T>(coefs[0]); newCoefs[1] = static_cast<T>(coefs[3]); @@ -326,10 +325,10 @@ namespace itk newCoefs[3] = static_cast<T>(coefs[6]); newCoefs[4] = static_cast<T>(coefs[7]); newCoefs[5] = static_cast<T>(coefs[8]); - + TensorType tensor(newCoefs); it.Set (tensor); - + ++ind; ++it; @@ -340,186 +339,186 @@ namespace itk } this->UpdateProgress ( 1.0 ); - + reader->Delete(); } - + template<class T, unsigned int TensorDimension, unsigned int ImageDimension> void TensorImageIO<T,TensorDimension,ImageDimension> ::ReadNrrd (const char* filename) - { + { typedef double PixelType; typedef itk::Image <PixelType, ImageDimension> ImageType; - + itk::NrrdImageIO::Pointer myIO = itk::NrrdImageIO::New(); myIO->SetFileName(filename); myIO->ReadImageInformation(); - + unsigned int numComponents = myIO->GetNumberOfComponents(); - + if( numComponents!=DegreesOfFreedom ) { char message[512]; sprintf(message, "Error: Number of components should be: %d, and actually is: %d.", DegreesOfFreedom, numComponents); throw itk::ExceptionObject(__FILE__,__LINE__,message); } - + if (myIO->GetPixelType() != itk::ImageIOBase::DIFFUSIONTENSOR3D) { - // Two cases : image is stored as float or double - if (myIO->GetComponentType() == itk::ImageIOBase::DOUBLE) - { - typedef itk::Vector<double, DegreesOfFreedom> VectorType; - typedef itk::Image <VectorType, ImageDimension> VectorImageType; - - typedef itk::ImageFileReader<VectorImageType> ReaderType; - - typename ReaderType::Pointer myReader = ReaderType::New(); - + // Two cases : image is stored as float or double + if (myIO->GetComponentType() == itk::ImageIOBase::DOUBLE) + { + typedef itk::Vector<double, DegreesOfFreedom> VectorType; + typedef itk::Image <VectorType, ImageDimension> VectorImageType; + + typedef itk::ImageFileReader<VectorImageType> ReaderType; + + typename ReaderType::Pointer myReader = ReaderType::New(); + myReader->SetFileName (filename); - myReader->SetImageIO(myIO); - - try - { + myReader->SetImageIO(myIO); + + try + { myReader->Update(); - } - catch (itk::ExceptionObject &e) - { + } + catch (itk::ExceptionObject &e) + { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::ReadNrrd()"); - } - - typename VectorImageType::Pointer myImage = myReader->GetOutput(); - - if( m_Output ) - { + } + + typename VectorImageType::Pointer myImage = myReader->GetOutput(); + + if( m_Output ) + { m_Output->Initialize(); - } - - typename ImageType::RegionType region = myImage->GetLargestPossibleRegion(); - - m_Output->SetRegions (region); - m_Output->SetSpacing (myImage->GetSpacing()); - m_Output->SetOrigin (myImage->GetOrigin()); - m_Output->SetDirection (myImage->GetDirection()); - - try - { - m_Output->Allocate(); - } - catch (itk::ExceptionObject &e) - { + } + + typename ImageType::RegionType region = myImage->GetLargestPossibleRegion(); + + m_Output->SetRegions (region); + m_Output->SetSpacing (myImage->GetSpacing()); + m_Output->SetOrigin (myImage->GetOrigin()); + m_Output->SetDirection (myImage->GetDirection()); + + try + { + m_Output->Allocate(); + } + catch (itk::ExceptionObject &e) + { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Read() during allocation."); - } - - itk::ImageRegionConstIteratorWithIndex<VectorImageType> itIn (myImage, myImage->GetLargestPossibleRegion()); - itk::ImageRegionIteratorWithIndex<TensorImageType> itOut(m_Output, m_Output->GetLargestPossibleRegion()); - - unsigned long numPixels = m_Output->GetLargestPossibleRegion().GetNumberOfPixels(); - unsigned long step = numPixels/100; - unsigned int ind = 0; - - this->UpdateProgress ( 0.0 ); - - while(!itOut.IsAtEnd()) - { - + } + + itk::ImageRegionConstIteratorWithIndex<VectorImageType> itIn (myImage, myImage->GetLargestPossibleRegion()); + itk::ImageRegionIteratorWithIndex<TensorImageType> itOut(m_Output, m_Output->GetLargestPossibleRegion()); + + unsigned long numPixels = m_Output->GetLargestPossibleRegion().GetNumberOfPixels(); + unsigned long step = numPixels/100; + unsigned int ind = 0; + + this->UpdateProgress ( 0.0 ); + + while(!itOut.IsAtEnd()) + { + VectorType vec = itIn.Get(); TensorType tensor; - + for( unsigned int j=0; j<DegreesOfFreedom; j++) { tensor[j] = static_cast<typename TensorType::ValueType>(vec[j]); } - - itOut.Set (tensor); + + itOut.Set (tensor); ++itOut; ++itIn; ++ind; - + if( (ind%step)==0 ) { this->UpdateProgress ( double(ind)/double(numPixels) ); } - - } + + } this->UpdateProgress ( 1.0 ); } else { - //Case image is in float - typedef itk::Vector<float, DegreesOfFreedom> VectorType; - typedef itk::Image <VectorType, ImageDimension> VectorImageType; - - typedef itk::ImageFileReader<VectorImageType> ReaderType; - - typename ReaderType::Pointer myReader = ReaderType::New(); - - myReader->SetFileName (filename); + //Case image is in float + typedef itk::Vector<float, DegreesOfFreedom> VectorType; + typedef itk::Image <VectorType, ImageDimension> VectorImageType; + + typedef itk::ImageFileReader<VectorImageType> ReaderType; + + typename ReaderType::Pointer myReader = ReaderType::New(); + + myReader->SetFileName (filename); myReader->SetImageIO(myIO); - - try - { + + try + { myReader->Update(); - } - catch (itk::ExceptionObject &e) - { + } + catch (itk::ExceptionObject &e) + { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::ReadNrrd()"); - } - - typename VectorImageType::Pointer myImage = myReader->GetOutput(); - - if( m_Output ) - { + } + + typename VectorImageType::Pointer myImage = myReader->GetOutput(); + + if( m_Output ) + { m_Output->Initialize(); - } - - typename ImageType::RegionType region = myImage->GetLargestPossibleRegion(); - - m_Output->SetRegions (region); - m_Output->SetSpacing (myImage->GetSpacing()); - m_Output->SetOrigin (myImage->GetOrigin()); - m_Output->SetDirection (myImage->GetDirection()); - + } + + typename ImageType::RegionType region = myImage->GetLargestPossibleRegion(); + + m_Output->SetRegions (region); + m_Output->SetSpacing (myImage->GetSpacing()); + m_Output->SetOrigin (myImage->GetOrigin()); + m_Output->SetDirection (myImage->GetDirection()); + try - { - m_Output->Allocate(); - } - catch (itk::ExceptionObject &e) - { + { + m_Output->Allocate(); + } + catch (itk::ExceptionObject &e) + { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Read() during allocation."); - } - - itk::ImageRegionConstIteratorWithIndex<VectorImageType> itIn (myImage, myImage->GetLargestPossibleRegion()); - itk::ImageRegionIteratorWithIndex<TensorImageType> itOut(m_Output, m_Output->GetLargestPossibleRegion()); + } + + itk::ImageRegionConstIteratorWithIndex<VectorImageType> itIn (myImage, myImage->GetLargestPossibleRegion()); + itk::ImageRegionIteratorWithIndex<TensorImageType> itOut(m_Output, m_Output->GetLargestPossibleRegion()); unsigned long numPixels = m_Output->GetLargestPossibleRegion().GetNumberOfPixels(); - unsigned long step = numPixels/100; - unsigned int ind = 0; - - this->UpdateProgress ( 0.0 ); - - while(!itOut.IsAtEnd()) + unsigned long step = numPixels/100; + unsigned int ind = 0; + + this->UpdateProgress ( 0.0 ); + + while(!itOut.IsAtEnd()) { - VectorType vec = itIn.Get(); + VectorType vec = itIn.Get(); TensorType tensor; - + for( unsigned int j=0; j<DegreesOfFreedom; j++) { tensor[j] = static_cast<typename TensorType::ValueType>(vec[j]); } itOut.Set (tensor); - - ++itOut; + + ++itOut; ++itIn; ++ind; @@ -527,11 +526,11 @@ namespace itk { this->UpdateProgress ( double(ind)/double(numPixels) ); } - - } + + } this->UpdateProgress ( 1.0 ); - } + } } else { @@ -539,52 +538,52 @@ namespace itk { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: Nrrd format with diffusiontensor3d only supports 3x3 tensors."); } - - + + typedef itk::DiffusionTensor3D<double> DiffPixelType; typedef itk::Image <DiffPixelType, ImageDimension> DiffImageType; - + typedef itk::ImageFileReader<DiffImageType> ReaderType; - + typename ReaderType::Pointer myReader = ReaderType::New(); myReader->SetFileName (filename); myReader->SetImageIO(myIO); - + try { myReader->Update(); } catch (itk::ExceptionObject &e) { - std::cerr << e; - throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::ReadNrrd()"); + std::cerr << e; + throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::ReadNrrd()"); } - + typename DiffImageType::Pointer myImage = myReader->GetOutput(); if( m_Output ) { - m_Output->Initialize(); + m_Output->Initialize(); } typename DiffImageType::RegionType region = myImage->GetLargestPossibleRegion(); - + m_Output->SetRegions (region); m_Output->SetSpacing (myImage->GetSpacing()); m_Output->SetOrigin (myImage->GetOrigin()); m_Output->SetDirection (myImage->GetDirection()); - + try { - m_Output->Allocate(); + m_Output->Allocate(); } catch (itk::ExceptionObject &e) { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Read() during allocation."); } - + itk::ImageRegionConstIteratorWithIndex<DiffImageType> itIn (myImage, myImage->GetLargestPossibleRegion()); itk::ImageRegionIteratorWithIndex<TensorImageType> itOut(m_Output, m_Output->GetLargestPossibleRegion()); @@ -593,31 +592,31 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while(!itOut.IsAtEnd()) { - + DiffPixelType d3d = itIn.Get(); - TensorType tensor; - for( int j=0; j<3; j++) + TensorType tensor; + for( int j=0; j<3; j++) { - for(int i=0; i<=j; i++) + for(int i=0; i<=j; i++) { - tensor.SetComponent ( i, j, static_cast<T>(d3d (i,j))); + tensor.SetComponent ( i, j, static_cast<T>(d3d (i,j))); } } itOut.Set (tensor); - + ++itOut; - ++itIn; - ++ind; + ++itIn; + ++ind; if( (ind%step)==0 ) { this->UpdateProgress ( double(ind)/double(numPixels) ); } - + } this->UpdateProgress ( 1.0 ); @@ -631,30 +630,30 @@ namespace itk { typedef T PixelType; typedef itk::Image <PixelType, ImageDimension> ImageType; - + itk::NiftiImageIO::Pointer myIO = itk::NiftiImageIO::New(); myIO->SetFileName(filename); myIO->ReadImageInformation(); - - + + unsigned int numComponents = myIO->GetNumberOfComponents(); - + if( numComponents!=DegreesOfFreedom ) { char message[512]; sprintf(message, "Error: Number of components should be: %d, and actually is: %d.", DegreesOfFreedom, numComponents); throw itk::ExceptionObject(__FILE__,__LINE__,message); } - + typedef itk::Vector<PixelType, DegreesOfFreedom> VectorType; typedef itk::Image <VectorType, ImageDimension> VectorImageType; - + typedef itk::ImageFileReader<VectorImageType> ReaderType; - + typename ReaderType::Pointer myReader = ReaderType::New(); myReader->SetFileName (filename); myReader->SetImageIO(myIO); - + try { myReader->Update(); @@ -664,32 +663,32 @@ namespace itk std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::ReadNifti()"); } - - + + typename VectorImageType::Pointer myImage = myReader->GetOutput(); - + if( m_Output ) { m_Output->Initialize(); } typename ImageType::RegionType region = myImage->GetLargestPossibleRegion(); - + m_Output->SetRegions (region); m_Output->SetSpacing (myImage->GetSpacing()); m_Output->SetOrigin (myImage->GetOrigin()); m_Output->SetDirection (myImage->GetDirection()); - + try { - m_Output->Allocate(); + m_Output->Allocate(); } catch (itk::ExceptionObject &e) { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Read() during allocation."); } - + itk::ImageRegionConstIteratorWithIndex<VectorImageType> itIn (myImage, myImage->GetLargestPossibleRegion()); itk::ImageRegionIteratorWithIndex<TensorImageType> itOut(m_Output, m_Output->GetLargestPossibleRegion()); @@ -698,20 +697,20 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while(!itOut.IsAtEnd()) { - + VectorType vec = itIn.Get(); TensorType tensor; - + for( unsigned int j=0; j<DegreesOfFreedom; j++) { tensor[j] = vec[j]; } - + itOut.Set (tensor); - + ++itOut; ++itIn; ++ind; @@ -720,14 +719,14 @@ namespace itk { this->UpdateProgress ( double(ind)/double(numPixels) ); } - + } this->UpdateProgress ( 1.0 ); - + } - + template<class T, unsigned int TensorDimension, unsigned int ImageDimension> void TensorImageIO<T,TensorDimension,ImageDimension> @@ -736,30 +735,30 @@ namespace itk typedef T PixelType; typedef itk::Image <PixelType, ImageDimension> ImageType; - + itk::MetaImageIO::Pointer myIO = itk::MetaImageIO::New(); myIO->SetFileName(filename); myIO->ReadImageInformation(); - - + + unsigned int numComponents = myIO->GetNumberOfComponents(); - + if( numComponents!=DegreesOfFreedom ) { char message[512]; sprintf(message, "Error: Number of components should be: %d, and actually is: %d.", DegreesOfFreedom, numComponents); throw itk::ExceptionObject(__FILE__,__LINE__,message); } - + typedef itk::Vector<PixelType, DegreesOfFreedom> VectorType; typedef itk::Image <VectorType, ImageDimension> VectorImageType; - + typedef itk::ImageFileReader<VectorImageType> ReaderType; - + typename ReaderType::Pointer myReader = ReaderType::New(); myReader->SetFileName (filename); myReader->SetImageIO(myIO); - + try { myReader->Update(); @@ -769,32 +768,32 @@ namespace itk std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::ReadNifti()"); } - - + + typename VectorImageType::Pointer myImage = myReader->GetOutput(); - + if( m_Output ) { m_Output->Initialize(); } typename ImageType::RegionType region = myImage->GetLargestPossibleRegion(); - + m_Output->SetRegions (region); m_Output->SetSpacing (myImage->GetSpacing()); m_Output->SetOrigin (myImage->GetOrigin()); m_Output->SetDirection (myImage->GetDirection()); - + try { - m_Output->Allocate(); + m_Output->Allocate(); } catch (itk::ExceptionObject &e) { std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Read() during allocation."); } - + itk::ImageRegionConstIteratorWithIndex<VectorImageType> itIn (myImage, myImage->GetLargestPossibleRegion()); itk::ImageRegionIteratorWithIndex<TensorImageType> itOut(m_Output, m_Output->GetLargestPossibleRegion()); @@ -803,20 +802,20 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while(!itOut.IsAtEnd()) { - + VectorType vec = itIn.Get(); TensorType tensor; - + for( unsigned int j=0; j<DegreesOfFreedom; j++) { tensor[j] = static_cast<typename TensorType::ValueType>(vec[j]); } - + itOut.Set (tensor); - + ++itOut; ++itIn; ++ind; @@ -825,16 +824,16 @@ namespace itk { this->UpdateProgress ( double(ind)/double(numPixels) ); } - + } this->UpdateProgress ( 1.0 ); } - - - - + + + + template<class T, unsigned int TensorDimension, unsigned int ImageDimension> void @@ -846,20 +845,20 @@ namespace itk { throw itk::ExceptionObject(__FILE__,__LINE__,"Error: input is not set."); } - + if(m_FileName == "") { throw itk::ExceptionObject(__FILE__,__LINE__,"Error: FileName is not set."); } - - - vtksys_stl::string ext = vtksys::SystemTools::GetFilenameExtension (m_FileName.c_str()); + + + std::string ext = vtksys::SystemTools::GetFilenameExtension (m_FileName.c_str()); if( ext=="" ) { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: Extension is not set."); } - - + + try { if ( CheckExtension(m_FileName.c_str(), ".vtk") ) @@ -893,10 +892,10 @@ namespace itk std::cerr << e; throw itk::ExceptionObject (__FILE__,__LINE__,"Error in TensorImageIO::Write()"); } - + } - - + + template<class T, unsigned int TensorDimension, unsigned int ImageDimension> void @@ -909,14 +908,14 @@ namespace itk { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: VTK only supports 3D images and 3x3 tensors."); } - - + + typename TensorImageType::SizeType size = m_Input->GetLargestPossibleRegion().GetSize(); int numVoxels = 1; int dims[3]; double origin[3]; double spacing[3]; - + for(unsigned int i=0; i<3; i++) { numVoxels *= size[i]; @@ -924,8 +923,8 @@ namespace itk origin[i] = m_Input->GetOrigin()[i]; spacing[i] = m_Input->GetSpacing()[i]; } - - + + vtkDoubleArray* data = vtkDoubleArray::New(); data->SetNumberOfComponents(9); data->SetNumberOfTuples(numVoxels); @@ -937,7 +936,7 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while( !it.IsAtEnd() ) { TensorType tensor = it.Get(); @@ -951,7 +950,7 @@ namespace itk buffer[6] = tensor.GetNthComponent(3); buffer[7] = tensor.GetNthComponent(4); buffer[8] = tensor.GetNthComponent(5); - + data->SetTuple(ind,buffer); ++ind; ++it; @@ -963,19 +962,20 @@ namespace itk } this->UpdateProgress ( 1.0 ); - + vtkStructuredPoints* tensors = vtkStructuredPoints::New(); tensors->SetDimensions(dims); tensors->SetSpacing(spacing); tensors->SetOrigin(origin); tensors->GetPointData()->SetTensors(data); data->Delete(); - tensors->Update(); + // tensors dhould be updated if needed by the witer. + // tensors->Update(); vtkStructuredPointsWriter* writer = vtkStructuredPointsWriter::New(); writer->SetFileTypeToBinary(); writer->SetFileName( m_FileName.c_str() ); - writer->SetInput(tensors); + writer->SetInputData(tensors); writer->Write(); tensors->Delete(); @@ -983,13 +983,13 @@ namespace itk } - + template<class T, unsigned int TensorDimension, unsigned int ImageDimension> void TensorImageIO<T,TensorDimension,ImageDimension> ::WriteNrrd (const char* filename) - { + { typedef itk::Vector<double, DegreesOfFreedom> VectorType; typedef itk::Image <VectorType, ImageDimension> ImageType; @@ -1016,7 +1016,7 @@ namespace itk typedef ImageRegionConstIterator<TensorImageType> IteratorType; IteratorType it(m_Input,m_Input->GetLargestPossibleRegion()); - + itk::ImageRegionIteratorWithIndex<ImageType> itOut(myTensorImage, myTensorImage->GetLargestPossibleRegion()); unsigned long numPixels = myTensorImage->GetLargestPossibleRegion().GetNumberOfPixels(); @@ -1024,18 +1024,18 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while( !it.IsAtEnd() ) { TensorType tensor = it.Get(); VectorType vec; - + // storing convention is in line first for( unsigned int i=0; i<DegreesOfFreedom; i++) { vec[i] = static_cast<double>(tensor[i]); - } - + } + itOut.Set (vec); ++it; @@ -1046,23 +1046,23 @@ namespace itk { this->UpdateProgress ( double(ind)/double(numPixels) ); } - + } this->UpdateProgress ( 1.0 ); - - + + typename itk::ImageFileWriter<ImageType>::Pointer myWriter = itk::ImageFileWriter<ImageType>::New(); myWriter->SetFileName(filename); myWriter->SetInput(myTensorImage); try { - myWriter->Write(); + myWriter->Write(); } catch(itk::ExceptionObject &e) { - std::cerr << e; - throw itk::ExceptionObject(__FILE__,__LINE__,"Error in TensorImageIO::WriteNrrd()"); + std::cerr << e; + throw itk::ExceptionObject(__FILE__,__LINE__,"Error in TensorImageIO::WriteNrrd()"); } } @@ -1073,15 +1073,15 @@ namespace itk void TensorImageIO<T,TensorDimension,ImageDimension> ::WriteNifti (const char* filename) - { - + { + typedef itk::Vector<T, DegreesOfFreedom> VectorType; typedef itk::Image<VectorType, ImageDimension> VectorImageType; typename VectorImageType::Pointer myTensorImage = VectorImageType::New(); - + typename TensorImageType::RegionType region = m_Input->GetLargestPossibleRegion(); - + myTensorImage->SetRegions (region); myTensorImage->SetSpacing (m_Input->GetSpacing()); myTensorImage->SetOrigin (m_Input->GetOrigin()); @@ -1100,7 +1100,7 @@ namespace itk typedef ImageRegionConstIterator<TensorImageType> IteratorType; IteratorType it(m_Input,m_Input->GetLargestPossibleRegion()); - + itk::ImageRegionIteratorWithIndex<VectorImageType> itOut(myTensorImage, myTensorImage->GetLargestPossibleRegion()); unsigned long numPixels = myTensorImage->GetLargestPossibleRegion().GetNumberOfPixels(); @@ -1108,19 +1108,19 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while( !it.IsAtEnd() ) { TensorType tensor = it.Get(); VectorType vec; - + for( unsigned int i=0; i<DegreesOfFreedom; i++) { vec[i] = tensor[i]; } - + itOut.Set (vec); - + ++it; ++itOut; ++ind; @@ -1129,23 +1129,23 @@ namespace itk { this->UpdateProgress ( double(ind)/double(numPixels) ); } - + } this->UpdateProgress ( 1.0 ); - - + + typename itk::ImageFileWriter<VectorImageType>::Pointer myWriter = itk::ImageFileWriter<VectorImageType>::New(); myWriter->SetFileName(filename); myWriter->SetInput(myTensorImage); try { - myWriter->Write(); + myWriter->Write(); } catch(itk::ExceptionObject &e) { - std::cerr << e; - throw itk::ExceptionObject(__FILE__,__LINE__,"Error in TensorImageIO::WriteNrrd()"); + std::cerr << e; + throw itk::ExceptionObject(__FILE__,__LINE__,"Error in TensorImageIO::WriteNrrd()"); } } @@ -1156,14 +1156,14 @@ namespace itk TensorImageIO<T,TensorDimension,ImageDimension> ::WriteMha (const char* filename) { - + this->WriteNifti (filename); } - - - + + + template<class T, unsigned int TensorDimension, unsigned int ImageDimension> void TensorImageIO<T,TensorDimension,ImageDimension> @@ -1171,18 +1171,18 @@ namespace itk { #ifdef TTK_USE_MIPS SizeType size = m_Input->GetLargestPossibleRegion().GetSize(); - + int numVox=1; for( unsigned int i=0; i<ImageDimension; i++) { numVox *= size[i]; } - + if( numVox==0 ) { throw itk::ExceptionObject(__FILE__,__LINE__,"Error: Tensor field has null dimensions."); } - + SpacingType spacing = m_Input->GetSpacing(); PointType origin = m_Input->GetOrigin(); yav::Inrimage* inr = new yav::Inrimage(size[0], @@ -1204,23 +1204,23 @@ namespace itk orig[2] = origin[2]; inr->setTranslation(orig); - + // we only use the rotational part of the // inrimage transformation matrix, as the translation is in the origin. // note that we take the Transformation matrix instead of the rotation instance // of the inrimage as the rotation instance is in voxel coordinates. // Note also that it SHOULD be a rigid rotation matrix, // in theory its det. should be 1... - + DirectionType direction = m_Input->GetDirection(); - yav::Rotation3D matrix; - + yav::Rotation3D matrix; + for (unsigned int i=0; i<3; i++) for (unsigned int j=0; j<3; j++) - matrix[i][j] = direction[i][j]; - + matrix[i][j] = direction[i][j]; + inr->setRotation(matrix); - + // Get the buffer double* buffer = (double*)inr->getData(); @@ -1231,7 +1231,7 @@ namespace itk unsigned int ind = 0; this->UpdateProgress ( 0.0 ); - + while(!it.IsAtEnd()) { int offset = (int)(TensorType::NDegreesOfFreedom)*ind; @@ -1240,7 +1240,7 @@ namespace itk { buffer[offset + i] = Tens.GetNthComponent (i); } - + ++it; ++ind; @@ -1249,21 +1249,21 @@ namespace itk this->UpdateProgress ( double(ind)/double(numVox) ); } } - + try - { + { inr->write(m_FileName.c_str()); } catch (...) - { + { throw itk::ExceptionObject (__FILE__,__LINE__,"Error while writing."); } this->UpdateProgress ( 1.0 ); - + delete inr; #else - std::cerr << "TTK was not compiled with inrimage support, sorry!" << std::endl; + std::cerr << "TTK was not compiled with inrimage support, sorry!" << std::endl; #endif }