diff --git a/Applications/Diffusion/MitkDiffusion.cpp b/Applications/Diffusion/MitkDiffusion.cpp index 46f3729715..78cbd7aaed 100644 --- a/Applications/Diffusion/MitkDiffusion.cpp +++ b/Applications/Diffusion/MitkDiffusion.cpp @@ -1,134 +1,137 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include class QtSafeApplication : public QtSingleApplication { public: QtSafeApplication(int& argc, char** argv) : QtSingleApplication(argc, argv) {} /** * Reimplement notify to catch unhandled exceptions and open an error message. * * @param receiver * @param event * @return */ bool notify(QObject* receiver, QEvent* event) { QString msg; try { return QApplication::notify(receiver, event); } catch (Poco::Exception& e) { msg = QString::fromStdString(e.displayText()); } catch (std::exception& e) { msg = e.what(); } catch (...) { msg = "Unknown exception"; } QString text("An error occurred. You should save all data and quit the program to " "prevent possible data loss.\nSee the error log for details.\n\n"); text += msg; QMessageBox::critical(0, "Error", text); return false; } }; int main(int argc, char** argv) { + + QtSafeApplication::setAttribute( Qt::AA_X11InitThreads ); + // Create a QApplication instance first QtSafeApplication qSafeApp(argc, argv); qSafeApp.setApplicationName("MitkDiffusion"); qSafeApp.setOrganizationName("DKFZ"); bool showSplashScreen(true); QPixmap pixmap( ":/splash/splashscreen.png" ); QSplashScreen splash( pixmap ); splash.setMask( pixmap.mask() ); splash.setWindowFlags( Qt::SplashScreen | Qt::WindowStaysOnTopHint | Qt::FramelessWindowHint ); if (showSplashScreen) { splash.show(); qSafeApp.sendPostedEvents(); qSafeApp.processEvents(); qSafeApp.flush(); QTimer::singleShot(4000, &splash, SLOT(close()) ); } // This function checks if an instance is already running // and either sends a message to it (containing the command // line arguments) or checks if a new instance was forced by // providing the BlueBerry.newInstance command line argument. // In the latter case, a path to a temporary directory for // the new application's storage directory is returned. QString storageDir = handleNewAppInstance(&qSafeApp, argc, argv, "BlueBerry.newInstance"); // These paths replace the .ini file and are tailored for installation // packages created with CPack. If a .ini file is presented, it will // overwrite the settings in MapConfiguration Poco::Path basePath(argv[0]); basePath.setFileName(""); Poco::Path provFile(basePath); provFile.setFileName("MitkDiffusion.provisioning"); Poco::Util::MapConfiguration* diffConfig(new Poco::Util::MapConfiguration()); if (!storageDir.isEmpty()) { diffConfig->setString(berry::Platform::ARG_STORAGE_DIR, storageDir.toStdString()); } diffConfig->setString(berry::Platform::ARG_PROVISIONING, provFile.toString()); diffConfig->setString(berry::Platform::ARG_APPLICATION, "org.mitk.qt.diffusionimagingapp"); // Preload the org.mitk.gui.qt.ext plug-in (and hence also QmitkExt) to speed // up a clean-cache start. This also works around bugs in older gcc and glibc implementations, // which have difficulties with multiple dynamic opening and closing of shared libraries with // many global static initializers. It also helps if dependent libraries have weird static // initialization methods and/or missing de-initialization code. diffConfig->setString(berry::Platform::ARG_PRELOAD_LIBRARY, "liborg_mitk_gui_qt_ext"); return berry::Starter::Run(argc, argv, diffConfig); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp index c5c01f68ea..d11c505925 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp @@ -1,325 +1,327 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP #define MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP #include "mitkDWIHeadMotionCorrectionFilter.h" #include "itkSplitDWImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkPyramidImageRegistrationMethod.h" //#include "mitkRegistrationMethodITK4.h" #include "mitkImageToDiffusionImageSource.h" #include "mitkDiffusionImageCorrectionFilter.h" #include #include #include "mitkIOUtil.h" #include template< typename DiffusionPixelType> mitk::DWIHeadMotionCorrectionFilter ::DWIHeadMotionCorrectionFilter() : m_CurrentStep(0), m_Steps(100), m_IsInValidState(true), m_AbortRegistration(false), m_AverageUnweighted(true) { } template< typename DiffusionPixelType> void mitk::DWIHeadMotionCorrectionFilter ::GenerateData() { typedef itk::SplitDWImageFilter< DiffusionPixelType, DiffusionPixelType> SplitFilterType; DiffusionImageType* input = const_cast(this->GetInput(0)); typedef mitk::PyramidImageRegistrationMethod RegistrationMethod; // typedef mitk::RegistrationMethodITKV4 RegistrationMethod unsigned int numberOfSteps = input->GetVectorImage()->GetNumberOfComponentsPerPixel () ; m_Steps = numberOfSteps; // // (1) Extract the b-zero images to a 3d+t image, register them by NCorr metric and // rigid registration : they will then be used are reference image for registering // the gradient images // typedef itk::B0ImageExtractionToSeparateImageFilter< DiffusionPixelType, DiffusionPixelType> B0ExtractorType; typename B0ExtractorType::Pointer b0_extractor = B0ExtractorType::New(); b0_extractor->SetInput( input->GetVectorImage() ); b0_extractor->SetDirections( input->GetDirections() ); b0_extractor->Update(); mitk::Image::Pointer b0Image = mitk::Image::New(); b0Image->InitializeByItk( b0_extractor->GetOutput() ); b0Image->SetImportChannel( b0_extractor->GetOutput()->GetBufferPointer(), mitk::Image::CopyMemory ); // (2.1) Use the extractor to access the extracted b0 volumes mitk::ImageTimeSelector::Pointer t_selector = mitk::ImageTimeSelector::New(); t_selector->SetInput( b0Image ); t_selector->SetTimeNr(0); t_selector->Update(); // first unweighted image as reference space for the registration mitk::Image::Pointer b0referenceImage = t_selector->GetOutput(); const unsigned int numberOfb0Images = b0Image->GetTimeSteps(); // register b-zeros only if the flag to average is set, use the first one otherwise if( m_AverageUnweighted && numberOfb0Images > 1) { RegistrationMethod::Pointer registrationMethod = RegistrationMethod::New(); registrationMethod->SetFixedImage( b0referenceImage ); registrationMethod->SetTransformToRigid(); // the unweighted images are of same modality registrationMethod->SetCrossModalityOff(); // use the advanced (windowed sinc) interpolation registrationMethod->SetUseAdvancedInterpolation(true); // Initialize the temporary output image mitk::Image::Pointer registeredB0Image = b0Image->Clone(); mitk::ImageTimeSelector::Pointer t_selector2 = mitk::ImageTimeSelector::New(); t_selector2->SetInput( b0Image ); for( unsigned int i=1; iSetTimeNr(i); t_selector2->Update(); registrationMethod->SetMovingImage( t_selector2->GetOutput() ); try { MITK_INFO << " === (" << i <<"/"<< numberOfb0Images-1 << ") :: Starting registration"; registrationMethod->Update(); } catch( const itk::ExceptionObject& e) { m_IsInValidState = false; mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // import volume to the inter-results mitk::ImageWriteAccessor imac(registrationMethod->GetResampledMovingImage()); registeredB0Image->SetImportVolume( imac.GetData(), i, 0, mitk::Image::ReferenceMemory ); } // use the accumulateImageFilter as provided by the ItkAccumulateFilter method in the header file AccessFixedDimensionByItk_1(registeredB0Image, ItkAccumulateFilter, (4), b0referenceImage ); } // // (2) Split the diffusion image into a 3d+t regular image, extract only the weighted images // typename SplitFilterType::Pointer split_filter = SplitFilterType::New(); split_filter->SetInput (input->GetVectorImage() ); split_filter->SetExtractAllAboveThreshold(8, input->GetBValueMap() ); - try { split_filter->Update(); } catch( const itk::ExceptionObject &e) { m_IsInValidState = false; mitkThrow() << " Caught exception from SplitImageFilter : " << e.what(); } mitk::Image::Pointer splittedImage = mitk::Image::New(); splittedImage->InitializeByItk( split_filter->GetOutput() ); splittedImage->SetImportChannel( split_filter->GetOutput()->GetBufferPointer(), mitk::Image::CopyMemory ); // // (3) Use again the time-selector to access the components separately in order // to perform the registration of Image -> unweighted reference // RegistrationMethod::Pointer weightedRegistrationMethod = RegistrationMethod::New(); weightedRegistrationMethod->SetTransformToAffine(); weightedRegistrationMethod->SetCrossModalityOn(); // // - (3.1) Set the reference image // - a single b0 image // - average over the registered b0 images if multiple present // weightedRegistrationMethod->SetFixedImage( b0referenceImage ); // use the advanced (windowed sinc) interpolation weightedRegistrationMethod->SetUseAdvancedInterpolation(true); + weightedRegistrationMethod->SetVerboseOn(); // // - (3.2) Register all timesteps in the splitted image onto the first reference // unsigned int maxImageIdx = splittedImage->GetTimeSteps(); mitk::TimeGeometry* tsg = splittedImage->GetTimeGeometry(); mitk::ProportionalTimeGeometry* ptg = dynamic_cast(tsg); ptg->Expand(maxImageIdx+1); ptg->SetTimeStepGeometry( ptg->GetGeometryForTimeStep(0), maxImageIdx ); mitk::Image::Pointer registeredWeighted = mitk::Image::New(); registeredWeighted->Initialize( splittedImage->GetPixelType(0), *tsg ); // insert the first unweighted reference as the first volume - mitk::ImageWriteAccessor imac(b0referenceImage); - registeredWeighted->SetImportVolume( imac.GetData(), - 0,0, mitk::Image::CopyMemory ); + // in own scope to release the accessor asap after copy + { + mitk::ImageWriteAccessor imac(b0referenceImage); + registeredWeighted->SetImportVolume( imac.GetData(), + 0,0, mitk::Image::CopyMemory ); + } // mitk::Image::Pointer registeredWeighted = splittedImage->Clone(); // this time start at 0, we have only gradient images in the 3d+t file // the reference image comes form an other image mitk::ImageTimeSelector::Pointer t_selector_w = mitk::ImageTimeSelector::New(); t_selector_w->SetInput( splittedImage ); // store the rotation parts of the transformations in a vector typedef RegistrationMethod::TransformMatrixType MatrixType; std::vector< MatrixType > estimated_transforms; for( unsigned int i=0; iSetTimeNr(i); t_selector_w->Update(); weightedRegistrationMethod->SetMovingImage( t_selector_w->GetOutput() ); try { MITK_INFO << " === (" << i+1 <<"/"<< maxImageIdx << ") :: Starting registration"; weightedRegistrationMethod->Update(); } catch( const itk::ExceptionObject& e) { m_IsInValidState = false; mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // allow expansion mitk::ImageWriteAccessor imac(weightedRegistrationMethod->GetResampledMovingImage()); registeredWeighted->SetImportVolume( imac.GetData(), i+1, 0, mitk::Image::CopyMemory); estimated_transforms.push_back( weightedRegistrationMethod->GetLastRotationMatrix() ); } // // (4) Cast the resulting image back to an diffusion weighted image // typename DiffusionImageType::GradientDirectionContainerType *gradients = input->GetDirections(); typename DiffusionImageType::GradientDirectionContainerType::Pointer gradients_new = DiffusionImageType::GradientDirectionContainerType::New(); typename DiffusionImageType::GradientDirectionType bzero_vector; bzero_vector.fill(0); // compose the direction vector // - no direction for the first image // - correct ordering of the directions based on the index list gradients_new->push_back( bzero_vector ); typename SplitFilterType::IndexListType index_list = split_filter->GetIndexList(); typename SplitFilterType::IndexListType::const_iterator lIter = index_list.begin(); while( lIter != index_list.end() ) { gradients_new->push_back( gradients->at( *lIter ) ); ++lIter; } typename mitk::ImageToDiffusionImageSource< DiffusionPixelType >::Pointer caster = mitk::ImageToDiffusionImageSource< DiffusionPixelType >::New(); caster->SetImage( registeredWeighted ); caster->SetBValue( input->GetReferenceBValue() ); caster->SetGradientDirections( gradients_new.GetPointer() ); try { caster->Update(); } catch( const itk::ExceptionObject& e) { m_IsInValidState = false; MITK_ERROR << "Casting back to diffusion image failed: "; mitkThrow() << "Subprocess failed with exception: " << e.what(); } // // (5) Adapt the gradient directions according to the estimated transforms // typedef mitk::DiffusionImageCorrectionFilter< DiffusionPixelType > CorrectionFilterType; typename CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); OutputImagePointerType output = caster->GetOutput(); corrector->SetImage( output ); corrector->CorrectDirections( estimated_transforms ); // // (6) Pass the corrected image to the filters output port // m_CurrentStep += 1; this->GetOutput()->SetVectorImage(output->GetVectorImage()); this->GetOutput()->SetReferenceBValue(output->GetReferenceBValue()); this->GetOutput()->SetMeasurementFrame(output->GetMeasurementFrame()); this->GetOutput()->SetDirections(output->GetDirections()); this->GetOutput()->InitializeFromVectorImage(); - this->GetOutput()->Modified(); } #endif // MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp index 71ea41585a..ff255bb3e5 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp @@ -1,181 +1,526 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPyramidImageRegistrationMethod.h" #include "mitkException.h" #include "mitkImageAccessByItk.h" mitk::PyramidImageRegistrationMethod::PyramidImageRegistrationMethod() : m_FixedImage(NULL), m_MovingImage(NULL), m_CrossModalityRegistration(true), m_UseAffineTransform(true), m_UseWindowedSincInterpolator(false), m_UseNearestNeighborInterpolator(false), m_UseMask(false), m_EstimatedParameters(0), m_InitialParameters(0), m_Verbose(false) { } mitk::PyramidImageRegistrationMethod::~PyramidImageRegistrationMethod() { if( m_EstimatedParameters != NULL) { delete [] m_EstimatedParameters; } } void mitk::PyramidImageRegistrationMethod::SetFixedImage(mitk::Image::Pointer fixed) { if( fixed.IsNotNull() ) { m_FixedImage = fixed; } } void mitk::PyramidImageRegistrationMethod::SetMovingImage(mitk::Image::Pointer moving) { if( moving.IsNotNull() ) { m_MovingImage = moving; } } void mitk::PyramidImageRegistrationMethod::SetFixedImageMask(mitk::Image::Pointer mask) { m_FixedImageMask = mask; } void mitk::PyramidImageRegistrationMethod::Update() { if( m_MovingImage.IsNull() ) { mitkThrow() << " Moving image is null"; } if( m_FixedImage.IsNull() ) { mitkThrow() << " Moving image is null"; } unsigned int allowedDimension = 3; if( m_FixedImage->GetDimension() != allowedDimension || m_MovingImage->GetDimension() != allowedDimension ) { mitkThrow() << " Only 3D Images supported."; } // // One possibility: use the FixedDimesnionByItk, but this instantiates ALL possible // pixel type combinations! // AccessTwoImagesFixedDimensionByItk( m_FixedImage, m_MovingImage, RegisterTwoImages, 3); // in helper: TypeSubset : short, float - AccessTwoImagesFixedDimensionTypeSubsetByItk( m_FixedImage, m_MovingImage, RegisterTwoImages, 3); + AccessTwoImagesFixedDimensionTypeSubsetByItk( m_FixedImage, m_MovingImage, RegisterTwoImagesV4, 3); } mitk::PyramidImageRegistrationMethod::TransformMatrixType mitk::PyramidImageRegistrationMethod ::GetLastRotationMatrix() { TransformMatrixType output; if( m_EstimatedParameters == NULL ) { output.set_identity(); return output; } typedef itk::MatrixOffsetTransformBase< double, 3, 3> BaseTransformType; BaseTransformType::Pointer base_transform = BaseTransformType::New(); if( this->m_UseAffineTransform ) { typedef itk::AffineTransform< double > TransformType; TransformType::Pointer transform = TransformType::New(); TransformType::ParametersType affine_params( TransformType::ParametersDimension ); this->GetParameters( &affine_params[0] ); transform->SetParameters( affine_params ); base_transform = transform; } else { typedef itk::Euler3DTransform< double > RigidTransformType; RigidTransformType::Pointer rtransform = RigidTransformType::New(); RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); this->GetParameters( &rigid_params[0] ); rtransform->SetParameters( rigid_params ); base_transform = rtransform; } return base_transform->GetMatrix().GetVnlMatrix(); } mitk::Image::Pointer mitk::PyramidImageRegistrationMethod ::GetResampledMovingImage() { mitk::Image::Pointer output = mitk::Image::New(); //output->Initialize( this->m_FixedImage ); AccessFixedDimensionByItk_1( this->m_MovingImage, ResampleMitkImage, 3, output ); return output; } mitk::Image::Pointer mitk::PyramidImageRegistrationMethod::GetResampledMovingImage(mitk::Image::Pointer movingImage, double* transform) { mitk::Image::Pointer output = mitk::Image::New(); unsigned int dim = 12; if( !m_UseAffineTransform ) dim = 6; if (m_EstimatedParameters == NULL) m_EstimatedParameters = new double[dim]; double tmpParams[12]; // save and set temporal transformation values for( unsigned int i=0; i +#include +#include +#include +#include +#include +#include + +template +void mitk::PyramidImageRegistrationMethod:: +RegisterTwoImagesV4(itk::Image* itkImage1, itk::Image* itkImage2) +{ + // Basic typedefs needed + typedef typename itk::Image ItkImageTypeFixed; + typedef typename itk::Image ItkImageTypeMoving; + + typedef itk::MatrixOffsetTransformBase< double, VImageDimension1, VImageDimension2 > BaseTransformType; + + typedef itk::Image< double, 3> ItkRegistrationImageType; + typedef itk::CastImageFilter< ItkImageTypeFixed, ItkRegistrationImageType> FixedCastFilterType; + typedef itk::CastImageFilter< ItkImageTypeMoving, ItkRegistrationImageType> MovingCastFilterType; + + typedef typename itk::ImageRegistrationMethodv4< ItkRegistrationImageType, ItkRegistrationImageType, AffineTransformType> RegistrationType; + typedef typename itk::MattesMutualInformationImageToImageMetricv4< ItkRegistrationImageType, ItkRegistrationImageType> MIMetricType; + typedef typename itk::CorrelationImageToImageMetricv4< ItkRegistrationImageType, ItkRegistrationImageType > NCMetricType; + typedef typename itk::ImageToImageMetricv4< ItkRegistrationImageType, ItkRegistrationImageType > BaseMetricType; + + typedef itk::GradientDescentLineSearchOptimizerv4 OptimizerType; + + typename ItkImageTypeFixed::Pointer referenceImage = itkImage1; + typename ItkImageTypeMoving::Pointer movingImage = itkImage2; + + // initial parameter dimension ( affine = default ) + unsigned int paramDim = 12; + typename BaseTransformType::Pointer transform; + if( m_UseAffineTransform ) + { + transform = AffineTransformType::New(); + } + else + { + transform = RigidTransformType::New(); + paramDim = 6; + } + + typename BaseTransformType::ParametersType initialParams(paramDim); + + initialParams.Fill(0); + if( m_UseAffineTransform ) + { + initialParams[0] = initialParams[4] = initialParams[8] = 1; + } + + + // [Prepare registration] + // The ITKv4 Methods ( the MI Metric ) require double-type images so we need to perform cast first + typename FixedCastFilterType::Pointer caster_f = FixedCastFilterType::New(); + typename MovingCastFilterType::Pointer caster_m = MovingCastFilterType::New(); + + try + { + caster_f->SetInput(0, referenceImage ); + caster_m->SetInput(0, movingImage ); + + caster_f->Update(); + caster_m->Update(); + + } + catch( const itk::ExceptionObject &e ) + { + return ; + } + + // [Prepare Registration] + // Estimate the size of the pyramid based on image dimension + // here the image size on the topmost level must remain higher than the threshold + unsigned int max_pyramid_lvl = 4; + unsigned int max_schedule_val = 8; + const unsigned int min_required_image_size = 12; + typename ItkImageTypeFixed::RegionType::SizeType image_size = referenceImage->GetLargestPossibleRegion().GetSize(); + unsigned int min_value = std::min( image_size[0], std::min( image_size[1], image_size[2])); + + // Adapt also the optimizer settings to the number of levels + float optmaxstep = 10; + float optminstep = 0.1f; + unsigned int iterations = 40; + unsigned int convergence_win_size = 8; + double convergence_tolerance = 1e-5; + while( min_value / max_schedule_val < min_required_image_size ) + { + max_schedule_val /= 2; + max_pyramid_lvl--; + optmaxstep -= 2; + optminstep *= 0.5f; + convergence_win_size += 4; + convergence_tolerance *= 0.5; + } + + typename RegistrationType::ShrinkFactorsArrayType shrink_factors(max_pyramid_lvl); + shrink_factors.Fill(1); + + shrink_factors[0] = max_schedule_val; + for( unsigned int i=1; i::Pointer mi_estimator = + itk::RegistrationParameterScalesFromPhysicalShift< MIMetricType >::New(); + mi_estimator->SetMetric( metric ); + + optimizer->SetScalesEstimator( mi_estimator ); + + base_metric = metric; + } + else + { + typename NCMetricType::Pointer metric = NCMetricType::New(); + + typename itk::RegistrationParameterScalesFromPhysicalShift< NCMetricType >::Pointer nc_estimator = + itk::RegistrationParameterScalesFromPhysicalShift< NCMetricType >::New(); + + nc_estimator->SetMetric( metric ); + + optimizer->SetScalesEstimator( nc_estimator ); + + base_metric = metric; + } + + optimizer->SetDoEstimateLearningRateOnce( false ); + optimizer->SetDoEstimateLearningRateAtEachIteration( true ); + optimizer->SetNumberOfIterations( iterations ); + optimizer->SetConvergenceWindowSize( convergence_win_size ); + optimizer->SetMaximumStepSizeInPhysicalUnits( optmaxstep ); + optimizer->SetMinimumConvergenceValue( convergence_tolerance ); + + // line search options + optimizer->SetEpsilon( 1e-3 ); + optimizer->SetLowerLimit( 0.3 ); + optimizer->SetUpperLimit( 1.7 ); + optimizer->SetMaximumLineSearchIterations( 20 ); + + // add observer tag if verbose + unsigned long vopt_tag = 0; + if(m_Verbose) + { + MITK_INFO << " :: Starting at " << initialParams; + + OptimizerIterationCommandv4< OptimizerType >::Pointer iterationObserver = + OptimizerIterationCommandv4::New(); + + vopt_tag = optimizer->AddObserver( itk::IterationEvent(), iterationObserver ); + } + + + // Initializing by Geometry + if( !m_UseAffineTransform && m_InitializeByGeometry ) + { + typedef itk::CenteredVersorTransformInitializer< ItkImageTypeFixed, ItkImageTypeMoving> TransformInitializerType; + typename TransformInitializerType::TransformType::Pointer rigidTransform = TransformInitializerType::TransformType::New() ; + MITK_INFO << "Initializer starting at : " << rigidTransform->GetParameters(); + + typename TransformInitializerType::Pointer initializer = TransformInitializerType::New(); + initializer->SetTransform( rigidTransform); + initializer->SetFixedImage( referenceImage.GetPointer() ); + initializer->SetMovingImage( movingImage.GetPointer() ); + initializer->MomentsOn(); + initializer->InitializeTransform(); + + MITK_INFO << "Initialized Rigid position : " << rigidTransform->GetParameters(); + initialParams[3] = rigidTransform->GetParameters()[3]; + initialParams[4] = rigidTransform->GetParameters()[4]; + initialParams[5] = rigidTransform->GetParameters()[5]; + } + + // [Prepare Registration] + // Masking (Optional) + if( m_UseMask ) + { + typedef itk::Image BinaryImageType; + BinaryImageType::Pointer itkFixedImageMask = BinaryImageType::New(); + itk::ImageMaskSpatialObject< 3 >::Pointer fixedMaskSpatialObject = itk::ImageMaskSpatialObject< 3 >::New(); + + CastToItkImage( m_FixedImageMask, itkFixedImageMask); + itk::NotImageFilter::Pointer notFilter = itk::NotImageFilter::New(); + notFilter->SetInput(itkFixedImageMask); + notFilter->Update(); + fixedMaskSpatialObject->SetImage( notFilter->GetOutput() ); + base_metric->SetFixedImageMask( fixedMaskSpatialObject ); + } + + // [Prepare Registration] + // combine all components to set-up registration + typename RegistrationType::Pointer registration = RegistrationType::New(); + + registration->SetFixedImage( 0, caster_f->GetOutput() ); + registration->SetMovingImage( 0, caster_m->GetOutput() ); + registration->SetMetric( base_metric ); + registration->SetOptimizer( optimizer ); + registration->SetMovingInitialTransform( transform.GetPointer() ); + registration->SetNumberOfLevels(max_pyramid_lvl); + registration->SetShrinkFactorsPerLevel( shrink_factors ); + + // observe the pyramid level change in order to adapt parameters + typename PyramidOptControlCommandv4::Pointer pyramid_observer = + PyramidOptControlCommandv4::New(); + unsigned int pyramid_tag = registration->AddObserver( itk::InitializeEvent(), pyramid_observer ); + + + try + { + registration->Update(); + } + catch( const itk::ExceptionObject &e) + { + registration->Print( std::cout ); + + MITK_ERROR << "[Registration Update] Caught ITK exception: "; + MITK_ERROR("itk.exception") << e.what(); + + mitkThrow() << "Registration failed with exception: " << e.what(); + } + + // [Post Registration] + // Retrieve the last transformation parameters from the performed registration task + if( m_EstimatedParameters != NULL) + { + delete [] m_EstimatedParameters; + } + + m_EstimatedParameters = new double[paramDim]; + + typename BaseTransformType::ParametersType finalParameters = registration->GetOptimizer()->GetCurrentPosition(); + for( unsigned int i=0; iGetValue() << " :: " << finalParameters; + + if( m_Verbose ) + { + MITK_INFO("Termination") << optimizer->GetStopConditionDescription(); + } + + // remove optimizer tag if used + if( vopt_tag ) + { + optimizer->RemoveObserver( vopt_tag ); + } + + if( pyramid_tag ) + { + registration->RemoveObserver( pyramid_tag ); + } + + +} + +template< typename TPixel, unsigned int VDimension> +void mitk::PyramidImageRegistrationMethod:: +ResampleMitkImage( itk::Image* itkImage, + mitk::Image::Pointer& outputImage ) +{ + typedef typename itk::Image< TPixel, VDimension> ImageType; + typedef typename itk::ResampleImageFilter< ImageType, ImageType, double> ResampleImageFilterType; + + typedef itk::LinearInterpolateImageFunction< ImageType, double > InterpolatorType; + typename InterpolatorType::Pointer linear_interpolator = InterpolatorType::New(); + + + typedef itk::NearestNeighborInterpolateImageFunction< ImageType, double > NearestNeighborInterpolatorType; + typename NearestNeighborInterpolatorType::Pointer nn_interpolator = NearestNeighborInterpolatorType::New(); + + + typedef itk::WindowedSincInterpolateImageFunction< ImageType, 7> WindowedSincInterpolatorType; + typename WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New(); + + typename ImageType::Pointer reference_image = ImageType::New(); + CastToItkImage(m_FixedImage,reference_image); + + typedef itk::MatrixOffsetTransformBase< double, 3, 3> BaseTransformType; + BaseTransformType::Pointer base_transform = BaseTransformType::New(); + + if( this->m_UseAffineTransform ) + { + typedef itk::AffineTransform< double > TransformType; + TransformType::Pointer transform = TransformType::New(); + + TransformType::ParametersType affine_params( TransformType::ParametersDimension ); + this->GetParameters( &affine_params[0] ); + + transform->SetParameters( affine_params ); + + base_transform = transform; + } + // Rigid + else + { + typedef itk::Euler3DTransform< double > RigidTransformType; + RigidTransformType::Pointer rtransform = RigidTransformType::New(); + + RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); + this->GetParameters( &rigid_params[0] ); + + rtransform->SetParameters( rigid_params ); + + base_transform = rtransform; + } + + typename ResampleImageFilterType::Pointer resampler = ResampleImageFilterType::New(); + + resampler->SetInterpolator( linear_interpolator ); + if( m_UseWindowedSincInterpolator ) + resampler->SetInterpolator( sinc_interpolator ); + if ( m_UseNearestNeighborInterpolator) + resampler->SetInterpolator ( nn_interpolator ); + + //resampler->SetNumberOfThreads(1); + + resampler->SetInput( itkImage ); + resampler->SetTransform( base_transform ); + resampler->SetReferenceImage( reference_image ); + resampler->UseReferenceImageOn(); + + resampler->Update(); + + mitk::GrabItkImageMemory( resampler->GetOutput(), outputImage); +} diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h index 9e71d844b9..084b88a3db 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h @@ -1,602 +1,341 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKPYRAMIDIMAGEREGISTRATION_H #define MITKPYRAMIDIMAGEREGISTRATION_H #include #include #include "itkImageMaskSpatialObject.h" #include "itkNotImageFilter.h" #include #include #include "mitkImage.h" #include "mitkBaseProcess.h" #include "mitkPyramidRegistrationMethodHelper.h" #include #include "mitkImageToItk.h" #include "mitkImageCast.h" #include "mitkImageCaster.h" #include "mitkITKImageImport.h" #include "mitkIOUtil.h" namespace mitk { /** * @brief The PyramidImageRegistration class implements a multi-scale registration method * * The PyramidImageRegistration class is suitable for aligning (f.e.) brain MR images. The method offers two * transform types * - Rigid: optimizing translation and rotation only and * - Affine ( default ): with scaling in addition ( 12 DOF ) * - * The error metric is internally chosen based on the selected task type ( @sa SetCrossModalityOn ) + * The error metric is internally chosen based on the selected task type : \sa SetCrossModalityOn + * * It uses * - MattesMutualInformation for CrossModality=on ( default ) and * - NormalizedCorrelation for CrossModality=off. */ class MitkDiffusionCore_EXPORT PyramidImageRegistrationMethod : public itk::Object { public: /** Typedefs */ mitkClassMacro(PyramidImageRegistrationMethod, itk::Object) /** Smart pointer support */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) typedef itk::OptimizerParameters ParametersType; /** Typedef for the transformation matrix, corresponds to the InternalMatrixType from ITK transforms */ typedef vnl_matrix_fixed< double, 3, 3> TransformMatrixType; typedef itk::AffineTransform< double > AffineTransformType; typedef itk::Euler3DTransform< double > RigidTransformType; /** Registration is between modalities - will take MattesMutualInformation as error metric */ void SetCrossModalityOn() { m_CrossModalityRegistration = true; } /** Registration is between modalities - will take NormalizedCorrelation as error metric */ void SetCrossModalityOff() { m_CrossModalityRegistration = false; } /** Turn the cross-modality on/off */ void SetCrossModality(bool flag) { if( flag ) this->SetCrossModalityOn(); else this->SetCrossModalityOff(); } /** Controll the verbosity of the registration output * * for true, each iteration step of the optimizer is watched and passed to the std::cout */ void SetVerbose(bool flag) { if( flag ) this->SetVerboseOn(); else this->SetVerboseOff(); } /** Turn verbosity on, \sa SetVerbose */ void SetVerboseOn() { m_Verbose = true; } /** Turn verbosity off, \sa SetVerbose */ void SetVerboseOff() { m_Verbose = false; } /** A rigid ( 6dof : translation + rotation ) transform is optimized */ void SetTransformToRigid() { m_UseAffineTransform = false; } /** An affine ( 12dof : Rigid + scale + skew ) transform is optimized */ void SetTransformToAffine() { m_UseAffineTransform = true; } /** Input image, the reference one */ void SetFixedImage( mitk::Image::Pointer ); /** Input image, the one to be transformed */ void SetMovingImage( mitk::Image::Pointer ); /** Fixed image mask, excludes the masked voxels from the registration metric*/ void SetFixedImageMask( mitk::Image::Pointer mask); - + /** \brief Use the itk::CenteredVersorTransformInitializer to perform initialization step for the registration + * + * The initializer takes the geometries and their features to perform an initialization for the successive registration + */ void SetInitializeByGeometry(bool flag) { m_InitializeByGeometry = flag; } void Update(); /** * @brief Get the number of parameters optimized ( 12 or 6 ) * * @return number of paramters */ unsigned int GetNumberOfParameters() { unsigned int retValue = 12; if(!m_UseAffineTransform) retValue = 6; return retValue; } /** * @brief Copies the estimated parameters to the given array * @param paramArray, target array for copy, make sure to allocate enough space */ void GetParameters( double* paramArray) { if( m_EstimatedParameters == NULL ) { mitkThrow() << "No parameters were estimated yet, call Update() first."; } unsigned int dim = 12; if( !m_UseAffineTransform ) dim = 6; for( unsigned int i=0; i * * It returns identity if the internal parameters are not-yet allocated ( i.e. the filter did not run yet ) * * @return \sa TransformMatrixType */ TransformMatrixType GetLastRotationMatrix(); protected: PyramidImageRegistrationMethod(); ~PyramidImageRegistrationMethod(); /** Fixed image, used as reference for registration */ mitk::Image::Pointer m_FixedImage; /** Moving image, will be transformed */ mitk::Image::Pointer m_MovingImage; mitk::Image::Pointer m_FixedImageMask; bool m_CrossModalityRegistration; bool m_UseAffineTransform; bool m_UseWindowedSincInterpolator; bool m_UseNearestNeighborInterpolator; bool m_UseMask; double* m_EstimatedParameters; ParametersType m_InitialParameters; /** Control the verbosity of the regsitistration output */ bool m_Verbose; bool m_InitializeByGeometry; - template - void RegisterTwoImages(itk::Image* itkImage1, itk::Image* itkImage2) - { - typedef typename itk::Image< TPixel1, VImageDimension1> FixedImageType; - typedef typename itk::Image< TPixel2, VImageDimension2> MovingImageType; - typedef typename itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType; - typedef itk::RegularStepGradientDescentOptimizer OptimizerType; - - typedef itk::MatrixOffsetTransformBase< double, VImageDimension1, VImageDimension2 > BaseTransformType; - - typedef typename itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType > MMIMetricType; - typedef typename itk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType > NCMetricType; - typedef typename itk::ImageToImageMetric< FixedImageType, MovingImageType> BaseMetricType; - - typename itk::LinearInterpolateImageFunction::Pointer interpolator = - itk::LinearInterpolateImageFunction::New(); - - - typename BaseMetricType::Pointer metric; - - unsigned int glob_max_threads = itk::MultiThreader::GetGlobalMaximumNumberOfThreads(); - itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1); - - if( m_CrossModalityRegistration ) - { - metric = MMIMetricType::New(); - //metric->SetNumberOfThreads( 6 ); - } - else - { - metric = NCMetricType::New(); - } - - - - // initial parameter dimension ( affine = default ) - unsigned int paramDim = 12; - typename BaseTransformType::Pointer transform; - if( m_UseAffineTransform ) - { - transform = AffineTransformType::New(); - } - else - { - transform = RigidTransformType::New(); - paramDim = 6; - } - - typename BaseTransformType::ParametersType initialParams(paramDim); - - initialParams.Fill(0); - if( m_UseAffineTransform ) - { - initialParams[0] = initialParams[4] = initialParams[8] = 1; - } - -/* FIXME : Review removal, occured during rebasing on master <30 Jan> - if( m_InitialParameters.Size() == paramDim ) - { - initialParams = m_InitialParameters; - } -*/ - typename FixedImageType::Pointer referenceImage = itkImage1; - typename MovingImageType::Pointer movingImage = itkImage2; - typename FixedImageType::RegionType maskedRegion = referenceImage->GetLargestPossibleRegion(); - - typename RegistrationType::Pointer registration = RegistrationType::New(); - unsigned int max_pyramid_lvl = 3; - unsigned int max_schedule_val = 4; - typename FixedImageType::RegionType::SizeType image_size = referenceImage->GetLargestPossibleRegion().GetSize(); - unsigned int min_value = std::min( image_size[0], std::min( image_size[1], image_size[2])); - - // condition for the top level pyramid image - float optmaxstep = 6; - float optminstep = 0.05f; - unsigned int iterations = 100; - if( min_value / max_schedule_val < 12 ) - { - max_schedule_val /= 2; - max_pyramid_lvl--; - optmaxstep *= 0.25f; - optminstep *= 0.1f; - } - - typename RegistrationType::ScheduleType fixedSchedule(max_pyramid_lvl,3); - fixedSchedule.Fill(1); - - fixedSchedule[0][0] = max_schedule_val; - fixedSchedule[0][1] = max_schedule_val; - fixedSchedule[0][2] = max_schedule_val; - - for( unsigned int i=1; iSetScales( optScales ); - optimizer->SetInitialPosition( initialParams ); - optimizer->SetNumberOfIterations( iterations ); - optimizer->SetGradientMagnitudeTolerance( 1e-4 ); - optimizer->SetRelaxationFactor( 0.7 ); - optimizer->SetMaximumStepLength( optmaxstep ); - optimizer->SetMinimumStepLength( optminstep ); - - // Initializing by Geometry - if (!m_UseAffineTransform && m_InitializeByGeometry) - { - typedef itk::CenteredVersorTransformInitializer< FixedImageType, MovingImageType> TransformInitializerType; - typename TransformInitializerType::TransformType::Pointer rigidTransform = TransformInitializerType::TransformType::New() ; - MITK_INFO << "Initializer starting at : " << rigidTransform->GetParameters(); - - typename TransformInitializerType::Pointer initializer = TransformInitializerType::New(); - initializer->SetTransform( rigidTransform); - initializer->SetFixedImage( referenceImage.GetPointer() ); - initializer->SetMovingImage( movingImage.GetPointer() ); - initializer->MomentsOn(); - initializer->InitializeTransform(); - MITK_INFO << "Initialized Rigid position : " << rigidTransform->GetParameters(); - initialParams[3]=rigidTransform->GetParameters()[3]; - initialParams[4]=rigidTransform->GetParameters()[4]; - initialParams[5]=rigidTransform->GetParameters()[5]; - } - - // add observer tag if verbose */ - unsigned long vopt_tag = 0; - if(m_Verbose) - { - MITK_INFO << " :: Starting at " << initialParams; - MITK_INFO << " :: Scales = " << optScales; - - OptimizerIterationCommand< OptimizerType >::Pointer iterationObserver = - OptimizerIterationCommand::New(); - - vopt_tag = optimizer->AddObserver( itk::IterationEvent(), iterationObserver ); - } - - - // INPUT IMAGES - registration->SetFixedImage( referenceImage ); - registration->SetFixedImageRegion( maskedRegion ); - registration->SetMovingImage( movingImage ); - registration->SetSchedules( fixedSchedule, fixedSchedule); - // SET MASKED AREA - typedef itk::Image BinaryImageType; - BinaryImageType::Pointer itkFixedImageMask = BinaryImageType::New(); - itk::ImageMaskSpatialObject< 3 >::Pointer fixedMaskSpatialObject = itk::ImageMaskSpatialObject< 3 >::New(); - if (m_UseMask) - { - CastToItkImage(m_FixedImageMask,itkFixedImageMask); - itk::NotImageFilter::Pointer notFilter = itk::NotImageFilter::New(); - notFilter->SetInput(itkFixedImageMask); - notFilter->Update(); - fixedMaskSpatialObject->SetImage( notFilter->GetOutput() ); - metric->SetFixedImageMask( fixedMaskSpatialObject ); - } - // OTHER INPUTS - registration->SetMetric( metric ); - registration->SetOptimizer( optimizer ); - registration->SetTransform( transform.GetPointer() ); - registration->SetInterpolator( interpolator ); - registration->SetInitialTransformParameters( initialParams ); - - typename PyramidOptControlCommand::Pointer pyramid_observer = - PyramidOptControlCommand::New(); - - - unsigned int pyramid_tag = 0; - if(m_Verbose) - { - pyramid_tag = registration->AddObserver( itk::IterationEvent(), pyramid_observer ); - } - - try - { - registration->Update(); - } - catch (itk::ExceptionObject &e) - { - registration->Print( std::cout ); - /*registration->Print( std::cout ); - MITK_INFO << "========== reference "; - itkImage1->Print( std::cout ); - - MITK_INFO << "========== moving "; - itkImage2->Print( std::cout );*/ - - MITK_ERROR << "[Registration Update] Caught ITK exception: "; - mitkThrow() << "Registration failed with exception: " << e.what(); - } - - if( m_EstimatedParameters != NULL) - { - delete [] m_EstimatedParameters; - } - - m_EstimatedParameters = new double[paramDim]; - - typename BaseTransformType::ParametersType finalParameters = registration->GetLastTransformParameters(); - for( unsigned int i=0; iGetValue() << " :: " << finalParameters; - - // remove optimizer tag if used - if( vopt_tag ) - { - optimizer->RemoveObserver( vopt_tag ); - } - - if( pyramid_tag ) - { - registration->RemoveObserver( pyramid_tag ); - } - - itk::MultiThreader::SetGlobalMaximumNumberOfThreads( glob_max_threads ); - - } + /** + * @brief The method takes two itk::Images and performs a multi-scale registration on them + * + * Can be accessed by the AccessTwoImagesFixedDimensionByItk macro. + * + * @note Currently the access is currently reduced only for short and float ( the only types occuring with mitk::DiffusionImage ) and + * hence the method is accessed by means of the \sa AccessTwoImagesFixedDimensionTypeSubsetByItk macro. + */ +template + void RegisterTwoImagesV4(itk::Image* itkImage1, itk::Image* itkImage2); + /** + * @brief ResampleMitkImage applies the functionality of an the itk::ResampleImageFilter to the given + * mitk::Image. + * + * The API of the function is conform to the interface of \sa AccessByItk_1 macros. + */ template< typename TPixel, unsigned int VDimension> void ResampleMitkImage( itk::Image* itkImage, - mitk::Image::Pointer& outputImage ) - { - typedef typename itk::Image< TPixel, VDimension> ImageType; - typedef typename itk::ResampleImageFilter< ImageType, ImageType, double> ResampleImageFilterType; - - typedef itk::LinearInterpolateImageFunction< ImageType, double > InterpolatorType; - typename InterpolatorType::Pointer linear_interpolator = InterpolatorType::New(); - - - typedef itk::NearestNeighborInterpolateImageFunction< ImageType, double > NearestNeighborInterpolatorType; - typename NearestNeighborInterpolatorType::Pointer nn_interpolator = NearestNeighborInterpolatorType::New(); - - - typedef itk::WindowedSincInterpolateImageFunction< ImageType, 7> WindowedSincInterpolatorType; - typename WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New(); - - typename ImageType::Pointer reference_image = ImageType::New(); - CastToItkImage(m_FixedImage,reference_image); - - typedef itk::MatrixOffsetTransformBase< double, 3, 3> BaseTransformType; - BaseTransformType::Pointer base_transform = BaseTransformType::New(); - - if( this->m_UseAffineTransform ) - { - typedef itk::AffineTransform< double > TransformType; - TransformType::Pointer transform = TransformType::New(); - - TransformType::ParametersType affine_params( TransformType::ParametersDimension ); - this->GetParameters( &affine_params[0] ); - - transform->SetParameters( affine_params ); - - base_transform = transform; - } - // Rigid - else - { - typedef itk::Euler3DTransform< double > RigidTransformType; - RigidTransformType::Pointer rtransform = RigidTransformType::New(); - - RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); - this->GetParameters( &rigid_params[0] ); - - rtransform->SetParameters( rigid_params ); - - base_transform = rtransform; - } - - typename ResampleImageFilterType::Pointer resampler = ResampleImageFilterType::New(); - - resampler->SetInterpolator( linear_interpolator ); - if( m_UseWindowedSincInterpolator ) - resampler->SetInterpolator( sinc_interpolator ); - if ( m_UseNearestNeighborInterpolator) - resampler->SetInterpolator ( nn_interpolator ); - - //resampler->SetNumberOfThreads(1); - - resampler->SetInput( itkImage ); - resampler->SetTransform( base_transform ); - resampler->SetReferenceImage( reference_image ); - resampler->UseReferenceImageOn(); - - resampler->Update(); - - mitk::GrabItkImageMemory( resampler->GetOutput(), outputImage); - } + mitk::Image::Pointer& outputImage ); }; } // end namespace #endif // MITKPYRAMIDIMAGEREGISTRATION_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h index 8a772722fe..2bf6d6097c 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h @@ -1,160 +1,212 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKPYRAMIDREGISTRATIONMETHODHELPER_H #define MITKPYRAMIDREGISTRATIONMETHODHELPER_H #include #include #include #include #include #include #include #include #include #include #include #include #include #include "mitkImageAccessByItk.h" /** * @brief Provides same functionality as \a AccessTwoImagesFixedDimensionByItk for a subset of types * * For now, the subset defined is only short and float. */ #define AccessTwoImagesFixedDimensionTypeSubsetByItk(mitkImage1, mitkImage2, itkImageTypeFunction, dimension) \ { \ const mitk::PixelType& pixelType1 = mitkImage1->GetPixelType(); \ const mitk::PixelType& pixelType2 = mitkImage2->GetPixelType(); \ const mitk::Image* constImage1 = mitkImage1; \ const mitk::Image* constImage2 = mitkImage2; \ mitk::Image* nonConstImage1 = const_cast(constImage1); \ mitk::Image* nonConstImage2 = const_cast(constImage2); \ nonConstImage1->Update(); \ nonConstImage2->Update(); \ _checkSpecificDimension(mitkImage1, (dimension)); \ _checkSpecificDimension(mitkImage2, (dimension)); \ _accessTwoImagesByItkForEach(itkImageTypeFunction, ((short, dimension))((unsigned short, dimension))((float, dimension))((double, dimension)), ((short, dimension))((unsigned short, dimension))((float, dimension))((double, dimension)) ) \ { \ std::string msg("Pixel type "); \ msg.append(pixelType1.GetComponentTypeAsString() ); \ msg.append(" or pixel type "); \ msg.append(pixelType2.GetComponentTypeAsString() ); \ msg.append(" is not in " MITK_PP_STRINGIZE(MITK_ACCESSBYITK_TYPES_DIMN_SEQ(dimension))); \ throw mitk::AccessByItkException(msg); \ } \ } /** * @brief The PyramidOptControlCommand class stears the step lenght of the * gradient descent optimizer in multi-scale registration methods */ template class PyramidOptControlCommand : public itk::Command { public: typedef itk::RegularStepGradientDescentOptimizer OptimizerType; mitkClassMacro(PyramidOptControlCommand, itk::Command) itkFactorylessNewMacro(Self) itkCloneMacro(Self) void Execute(itk::Object *caller, const itk::EventObject & /*event*/) { RegistrationType* registration = dynamic_cast< RegistrationType* >( caller ); + if( registration == NULL) + return; + MITK_DEBUG << "\t - Pyramid level " << registration->GetCurrentLevel(); if( registration->GetCurrentLevel() == 0 ) - return; + { MITK_WARN("OptCommand") << "Cast to registration failed"; + return; + } OptimizerType* optimizer = dynamic_cast< OptimizerType* >(registration->GetOptimizer()); + if( optimizer == NULL) + { MITK_WARN("OptCommand") << "Cast to optimizer failed"; + return; + } + MITK_INFO /*<< optimizer->GetStopConditionDescription() << "\n"*/ << optimizer->GetValue() << " : " << optimizer->GetCurrentPosition(); optimizer->SetMaximumStepLength( optimizer->GetMaximumStepLength() * 0.25f ); optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() * 0.1f ); // optimizer->SetNumberOfIterations( optimizer->GetNumberOfIterations() * 1.5f ); } void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/){} }; +#include + +template +class PyramidOptControlCommandv4 : public itk::Command +{ +public: + + typedef itk::GradientDescentLineSearchOptimizerv4 OptimizerType; + + mitkClassMacro(PyramidOptControlCommandv4, itk::Command) + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + void Execute(itk::Object *caller, const itk::EventObject & /*event*/) + { + RegistrationType* registration = dynamic_cast< RegistrationType* >( caller ); + + if( registration == NULL) + return; + + MITK_DEBUG << "\t - Pyramid level " << registration->GetCurrentLevel(); + if( registration->GetCurrentLevel() == 0 ) + return; + + OptimizerType* optimizer = dynamic_cast< OptimizerType* >( registration->GetOptimizer() ); + + if( optimizer == NULL) + { MITK_WARN("OptCommand4") << "Cast to optimizer failed"; + return; + } + + optimizer->SetNumberOfIterations( optimizer->GetNumberOfIterations() * 2.5 ); + optimizer->SetMaximumStepSizeInPhysicalUnits( optimizer->GetMaximumStepSizeInPhysicalUnits() * 0.4); + + MITK_INFO("Pyramid.Command.Iter") << optimizer->GetNumberOfIterations(); + MITK_INFO("Pyramid.Command.MaxStep") << optimizer->GetMaximumStepSizeInPhysicalUnits(); + + } + + void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/){} +}; + template class OptimizerIterationCommand : public itk::Command { public: mitkClassMacro(OptimizerIterationCommand, itk::Command) itkFactorylessNewMacro(Self) itkCloneMacro(Self) void Execute(itk::Object *caller, const itk::EventObject & /*event*/) { OptimizerType* optimizer = dynamic_cast< OptimizerType* >( caller ); unsigned int currentIter = optimizer->GetCurrentIteration(); MITK_INFO << "[" << currentIter << "] : " << optimizer->GetValue() << " : " << optimizer->GetCurrentPosition(); } void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/) { } }; template class OptimizerIterationCommandv4 : public itk::Command { public: itkNewMacro( OptimizerIterationCommandv4 ) void Execute(itk::Object *object, const itk::EventObject & event) { OptimizerType* optimizer = dynamic_cast< OptimizerType* >( object ); if( typeid( event ) != typeid( itk::IterationEvent ) ) { return; } unsigned int currentIter = optimizer->GetCurrentIteration(); MITK_INFO << "[" << currentIter << "] : " << optimizer->GetCurrentMetricValue() << " : " << optimizer->GetMetric()->GetParameters() ; //<< " : " << optimizer->GetScales(); } void Execute(const itk::Object * object, const itk::EventObject & event) { } }; #endif // MITKPYRAMIDREGISTRATIONMETHODHELPER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx index d4cf021597..75baf432c6 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.txx @@ -1,631 +1,631 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkImageRegionIterator.h" #include "itkImageRegionConstIterator.h" #include "mitkImageCast.h" #include "mitkDiffusionImage.h" #include "itkImageDuplicator.h" #include "itkTestingComparisonImageFilter.h" template mitk::DiffusionImage::DiffusionImage() : m_VectorImage(0), m_B_Value(-1.0), m_OriginalDirections(0), m_Directions(0) { MeasurementFrameType mf; for(int i=0; i<3; i++) for(int j=0; j<3; j++) mf[i][j] = 0; for(int i=0; i<3; i++) mf[i][i] = 1; m_MeasurementFrame = mf; } template mitk::DiffusionImage::DiffusionImage(const DiffusionImage & orig) : mitk::Image(orig), m_VectorImage( 0 ), m_B_Value(0), m_OriginalDirections(0), m_Directions(0) { // Deep copy VectorImage typename itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage( orig.m_VectorImage ); duplicator->Update(); GradientDirectionContainerType::ConstIterator origIt; GradientDirectionContainerType::ConstIterator origItEnd; // Deep Copy OrignalDirectioncontainer GradientDirectionContainerType::Pointer OriginalDirectionscontainer = GradientDirectionContainerType::New(); origIt = orig.GetDirectionsWithoutMeasurementFrame()->Begin(); origItEnd = orig.GetDirectionsWithoutMeasurementFrame()->End(); for(;origIt != origItEnd; ++origIt) OriginalDirectionscontainer->push_back(origIt.Value()); // Deep Copy Directioncontainer GradientDirectionContainerType::Pointer DirectionsContainer = GradientDirectionContainerType::New(); origIt = orig.GetDirections()->Begin(); origItEnd = orig.GetDirections()->End(); for(;origIt != origItEnd; ++origIt) DirectionsContainer->push_back(origIt.Value()); // Set member of the new clone m_VectorImage = duplicator->GetOutput(); m_B_Value = orig.GetReferenceBValue(); m_OriginalDirections = OriginalDirectionscontainer; m_Directions = DirectionsContainer; m_MeasurementFrame = orig.GetMeasurementFrame(); ApplyMeasurementFrame(); UpdateBValueMap(); InitializeFromVectorImage(); // m_VectorImage = duplicator->GetOutput(); // m_B_Value = orig.GetB_Value(); // m_MeasurementFrame = orig.GetMeasurementFrame(); // m_Directions = copyInContainer; // m_OriginalDirections = copyInContainer; } template mitk::DiffusionImage::~DiffusionImage() { // Remove Observer for m_Directions //RemoveObserver(); } template void mitk::DiffusionImage ::InitializeFromVectorImage() { if(!m_VectorImage || !m_Directions || m_B_Value==-1.0) { MITK_INFO << "DiffusionImage could not be initialized. Set all members first!" << std::endl; // Fehjlermedlung updaten return; } // find bzero index int firstZeroIndex = -1; for(GradientDirectionContainerType::ConstIterator it = m_Directions->Begin(); it != m_Directions->End(); ++it) { firstZeroIndex++; GradientDirectionType g = it.Value(); if(g[0] == 0 && g[1] == 0 && g[2] == 0 ) break; } typedef itk::Image ImgType; typename ImgType::Pointer img = ImgType::New(); img->SetSpacing( m_VectorImage->GetSpacing() ); // Set the image spacing img->SetOrigin( m_VectorImage->GetOrigin() ); // Set the image origin img->SetDirection( m_VectorImage->GetDirection() ); // Set the image direction img->SetLargestPossibleRegion( m_VectorImage->GetLargestPossibleRegion()); img->SetBufferedRegion( m_VectorImage->GetLargestPossibleRegion() ); img->Allocate(); int vecLength = m_VectorImage->GetVectorLength(); InitializeByItk( img.GetPointer(), 1, vecLength ); itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); itw.GoToBegin(); itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); itr.GoToBegin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(firstZeroIndex)); ++itr; ++itw; } // init SetImportVolume(img->GetBufferPointer()); m_DisplayIndex = firstZeroIndex; MITK_INFO << "Diffusion-Image successfully initialized."; } template void mitk::DiffusionImage ::SetDisplayIndexForRendering(int displayIndex) { int index = displayIndex; int vecLength = m_VectorImage->GetVectorLength(); index = index > vecLength-1 ? vecLength-1 : index; if( m_DisplayIndex != index ) { typedef itk::Image ImgType; typename ImgType::Pointer img = ImgType::New(); CastToItkImage(this, img); itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); itw.GoToBegin(); itk::ImageRegionConstIterator itr (m_VectorImage, m_VectorImage->GetLargestPossibleRegion() ); itr.GoToBegin(); while(!itr.IsAtEnd()) { itw.Set(itr.Get().GetElement(index)); ++itr; ++itw; } } m_DisplayIndex = index; } template bool mitk::DiffusionImage::AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision) { GradientDirectionType diff = g1 - g2; GradientDirectionType diff2 = g1 + g2; return diff.two_norm() < precision || diff2.two_norm() < precision; } template mitk::DiffusionImage::GradientDirectionContainerType::Pointer mitk::DiffusionImage::CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions) { // save old and construct new direction container GradientDirectionContainerType::Pointer newDirections = GradientDirectionContainerType::New(); // fill new direction container for(GradientDirectionContainerType::ConstIterator gdcitOld = directions->Begin(); gdcitOld != directions->End(); ++gdcitOld) { // already exists? bool found = false; for(GradientDirectionContainerType::ConstIterator gdcitNew = newDirections->Begin(); gdcitNew != newDirections->End(); ++gdcitNew) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { found = true; break; } } // if not found, add it to new container if(!found) { newDirections->push_back(gdcitOld.Value()); } } return newDirections; } template void mitk::DiffusionImage::AverageRedundantGradients(double precision) { GradientDirectionContainerType::Pointer newDirs = CalcAveragedDirectionSet(precision, m_Directions); GradientDirectionContainerType::Pointer newOriginalDirs = CalcAveragedDirectionSet(precision, m_OriginalDirections); // if sizes equal, we do not need to do anything in this function if(m_Directions->size() == newDirs->size() || m_OriginalDirections->size() == newOriginalDirs->size()) return; GradientDirectionContainerType::Pointer oldDirections = m_OriginalDirections; m_Directions = newDirs; m_OriginalDirections = newOriginalDirs; // new image typename ImageType::Pointer oldImage = m_VectorImage; m_VectorImage = ImageType::New(); m_VectorImage->SetSpacing( oldImage->GetSpacing() ); // Set the image spacing m_VectorImage->SetOrigin( oldImage->GetOrigin() ); // Set the image origin m_VectorImage->SetDirection( oldImage->GetDirection() ); // Set the image direction m_VectorImage->SetLargestPossibleRegion( oldImage->GetLargestPossibleRegion() ); m_VectorImage->SetVectorLength( m_Directions->size() ); m_VectorImage->SetBufferedRegion( oldImage->GetLargestPossibleRegion() ); m_VectorImage->Allocate(); // average image data that corresponds to identical directions itk::ImageRegionIterator< ImageType > newIt(m_VectorImage, m_VectorImage->GetLargestPossibleRegion()); newIt.GoToBegin(); itk::ImageRegionIterator< ImageType > oldIt(oldImage, oldImage->GetLargestPossibleRegion()); oldIt.GoToBegin(); // initial new value of voxel typename ImageType::PixelType newVec; newVec.SetSize(m_Directions->size()); newVec.AllocateElements(m_Directions->size()); std::vector > dirIndices; for(GradientDirectionContainerType::ConstIterator gdcitNew = m_Directions->Begin(); gdcitNew != m_Directions->End(); ++gdcitNew) { dirIndices.push_back(std::vector(0)); for(GradientDirectionContainerType::ConstIterator gdcitOld = oldDirections->Begin(); gdcitOld != oldDirections->End(); ++gdcitOld) { if(AreAlike(gdcitNew.Value(), gdcitOld.Value(), precision)) { //MITK_INFO << gdcitNew.Value() << " " << gdcitOld.Value(); dirIndices[gdcitNew.Index()].push_back(gdcitOld.Index()); } } } //int ind1 = -1; while(!newIt.IsAtEnd()) { // progress //typename ImageType::IndexType ind = newIt.GetIndex(); //ind1 = ind.m_Index[2]; // init new vector with zeros newVec.Fill(0.0); // the old voxel value with duplicates typename ImageType::PixelType oldVec = oldIt.Get(); for(unsigned int i=0; i void mitk::DiffusionImage::ApplyMeasurementFrame() { if(m_OriginalDirections.IsNull()) { // original direction container was not set return; } m_Directions = GradientDirectionContainerType::New(); int c = 0; for(GradientDirectionContainerType::ConstIterator gdcit = m_OriginalDirections->Begin(); gdcit != m_OriginalDirections->End(); ++gdcit) { vnl_vector vec = gdcit.Value(); vec = vec.pre_multiply(m_MeasurementFrame); m_Directions->InsertElement(c, vec); c++; } } template void mitk::DiffusionImage::UpdateBValueMap() { if(!m_B_ValueMap.empty()) m_B_ValueMap.clear(); GradientDirectionContainerType::ConstIterator gdcit; for( gdcit = this->m_Directions->Begin(); gdcit != this->m_Directions->End(); ++gdcit) - m_B_ValueMap[GetB_Value(gdcit.Index())].push_back(gdcit.Index()); + m_B_ValueMap[this->GetB_Value(gdcit.Index())].push_back(gdcit.Index()); } template float mitk::DiffusionImage::GetB_Value(unsigned int i) { if(i > m_Directions->Size()-1) return -1; if(m_Directions->ElementAt(i).one_norm() <= 0.0) { return 0; } else { double twonorm = m_Directions->ElementAt(i).two_norm(); double bval = m_B_Value*twonorm*twonorm; if (bval<0) bval = ceil(bval - 0.5); else bval = floor(bval + 0.5); return bval; } } template void mitk::DiffusionImage::SetDirections( GradientDirectionContainerType::Pointer directions ) { this->m_OriginalDirections = directions; ApplyMeasurementFrame(); UpdateBValueMap(); } template void mitk::DiffusionImage::SetDirections(const std::vector > &directions) { GradientDirectionContainerType::Pointer tempContainer = GradientDirectionContainerType::New(); for(unsigned int i=0; iInsertElement( i, directions[i].GetVnlVector() ); SetDirections(tempContainer); } template typename mitk::DiffusionImage::MeasurementFrameType mitk::DiffusionImage::GetMeasurementFrame() const { return m_MeasurementFrame; } template void mitk::DiffusionImage::SetMeasurementFrame( const MeasurementFrameType & mFrame ) { m_MeasurementFrame = mFrame; ApplyMeasurementFrame(); } template bool mitk::DiffusionImage::GetNumberOfBValues() const { return m_B_ValueMap.size(); } template const typename mitk::DiffusionImage::BValueMap & mitk::DiffusionImage::GetBValueMap() const { return m_B_ValueMap; } template float mitk::DiffusionImage::GetReferenceBValue() const { return m_B_Value; } template void mitk::DiffusionImage::SetReferenceBValue( float val ) { m_B_Value = val; } template typename mitk::DiffusionImage::GradientDirectionContainerTypePointer mitk::DiffusionImage::GetDirections() const { return m_Directions; } template typename mitk::DiffusionImage::GradientDirectionContainerTypePointer mitk::DiffusionImage::GetDirectionsWithoutMeasurementFrame() const { return m_OriginalDirections; } template typename itk::VectorImage::Pointer mitk::DiffusionImage::GetVectorImage() { return m_VectorImage; } template const typename itk::VectorImage::Pointer mitk::DiffusionImage::GetVectorImage() const { return m_VectorImage; } template void mitk::DiffusionImage::SetVectorImage(typename ImageType::Pointer image ) { m_VectorImage = image; } template inline bool mitk::Equal(const mitk::DiffusionImage& leftHandSide, const mitk::DiffusionImage &rightHandSide, ScalarType eps, bool verbose ) { bool returnValue = true; if(leftHandSide.GetReferenceBValue() != rightHandSide.GetReferenceBValue()) { if(verbose) MITK_INFO << "[( DiffusionImage )] Reference BValue is not Equal."; returnValue = false; } if(leftHandSide.GetMeasurementFrame() != rightHandSide.GetMeasurementFrame()) { if(verbose) MITK_INFO << "[( DiffusionImage )] MeasurementFrame is not Equal."; returnValue = false; } if(leftHandSide.GetBValueMap().size() != rightHandSide.GetBValueMap().size()) { if(verbose) MITK_INFO << "[( DiffusionImage )] BValue Map: Size is not Equal."; returnValue = false; } if(leftHandSide.GetNumberOfBValues() != rightHandSide.GetNumberOfBValues()) { if(verbose) MITK_INFO << "[( DiffusionImage )] BValue Map (GetNumberOfBValues): Size is not Equal."; returnValue = false; } // Slow testing area const mitk::Image* img1 = dynamic_cast(&leftHandSide); const mitk::Image* img2 = dynamic_cast(&rightHandSide); if(mitk::Equal(*img1,*img2,eps,verbose) == false) { if(verbose) MITK_INFO << "[( DiffusionImage )] Base-Class (mitk::Image) is not Equal."; returnValue = false; } { typename mitk::DiffusionImage::GradientDirectionContainerType::Iterator lhsIt = leftHandSide.GetDirectionsWithoutMeasurementFrame()->Begin(); typename mitk::DiffusionImage::GradientDirectionContainerType::Iterator lhsItEnd = leftHandSide.GetDirectionsWithoutMeasurementFrame()->End(); typename mitk::DiffusionImage::GradientDirectionContainerType::Iterator rhsIt = rightHandSide.GetDirectionsWithoutMeasurementFrame()->Begin(); for(;lhsIt != lhsItEnd;) { bool vectorNotEqual = false; for(unsigned int i = 0 ; i < lhsIt.Value().size(); i++) vectorNotEqual |= !mitk::Equal(lhsIt.Value().get(i),rhsIt.Value().get(i),0.0001); if(vectorNotEqual) { if(verbose) MITK_INFO << "[( DiffusionImage )] Original GradientDirections are not Equal."; returnValue = false; break; } ++rhsIt; ++lhsIt; } } { typename mitk::DiffusionImage::GradientDirectionContainerType::Iterator lhsIt = leftHandSide.GetDirections()->Begin(); typename mitk::DiffusionImage::GradientDirectionContainerType::Iterator lhsItEnd = leftHandSide.GetDirections()->End(); typename mitk::DiffusionImage::GradientDirectionContainerType::Iterator rhsIt = rightHandSide.GetDirections()->Begin(); for(;lhsIt != lhsItEnd;) { bool vectorNotEqual = false; for(unsigned int i = 0 ; i < lhsIt.Value().size(); i++) vectorNotEqual |= !mitk::Equal(lhsIt.Value().get(i),rhsIt.Value().get(i),0.0001); if(vectorNotEqual) { if(verbose) MITK_INFO << "[( DiffusionImage )] GradientDirections are not Equal."; returnValue = false; break; } ++rhsIt; ++lhsIt; } } { typename mitk::DiffusionImage::BValueMap rhsMap = rightHandSide.GetBValueMap(); typename mitk::DiffusionImage::BValueMap lhsMap = leftHandSide.GetBValueMap(); typename mitk::DiffusionImage::BValueMap::const_iterator lhsIt = lhsMap.begin(); typename mitk::DiffusionImage::BValueMap::const_iterator lhsItEnd = lhsMap.end(); typename mitk::DiffusionImage::BValueMap::const_iterator rhsIt = rhsMap.begin(); for(;lhsIt != lhsItEnd;) { if(lhsIt->first != rhsIt->first) { if(verbose) MITK_INFO << "[( DiffusionImage )] BValue Map: lhsKey " << lhsIt->first << " != rhsKey " << rhsIt->first << " is not Equal."; returnValue = false; break; } if(lhsIt->second.size() != rhsIt->second.size()) { if(verbose) MITK_INFO << "[( DiffusionImage )] BValue Map: Indices vector size is not Equal. (Key: " << lhsIt->first <<")"; returnValue = false; break; } typename mitk::DiffusionImage::IndicesVector::const_iterator lhsIndVecIt = lhsIt->second.begin(); typename mitk::DiffusionImage::IndicesVector::const_iterator lhsIndVecItEnd = lhsIt->second.end(); typename mitk::DiffusionImage::IndicesVector::const_iterator rhsIndVecIt = rhsIt->second.begin(); for(;lhsIndVecIt != lhsIndVecItEnd;) { if(*lhsIndVecIt != *rhsIndVecIt) { if(verbose) MITK_INFO << "[( DiffusionImage )] BValue Map: Indices are not Equal. (Key: " << lhsIt->first <<")"; returnValue = false; break; } ++rhsIndVecIt; ++lhsIndVecIt; } lhsIt++; rhsIt++; } } try{ itk::ImageRegionIterator< itk::VectorImage< TPixelType, 3 > > it1(leftHandSide.GetVectorImage(), leftHandSide.GetVectorImage()->GetLargestPossibleRegion()); itk::ImageRegionIterator< itk::VectorImage< TPixelType, 3 > > it2(rightHandSide.GetVectorImage(), rightHandSide.GetVectorImage()->GetLargestPossibleRegion()); while(!it1.IsAtEnd()) { if (it1.Get()!=it2.Get()){ if(verbose) MITK_INFO << "[( DiffusionImage )] Capsulated itk::VectorImage is not Equal."; returnValue = false; break; } ++it1; ++it2; } } catch(...) { if(verbose) MITK_INFO << "[( DiffusionImage )] Comparision of itk Vector image abort by exception."; returnValue = false; } return returnValue; } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp index 3771fe793b..7805c72b3c 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp @@ -1,60 +1,60 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "mitkIOUtil.h" #include "mitkDWIHeadMotionCorrectionFilter.h" typedef short DiffusionPixelType; typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; /** * @brief Custom test to provide CMD-line access to the mitk::DWIHeadMotionCorrectionFilter * * @param argv : Input and Output image full path */ int mitkDWHeadMotionCorrectionTest( int argc, char* argv[] ) { MITK_TEST_BEGIN("mitkDWHeadMotionCorrectionTest"); MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output."); - itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1); +// itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1); mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); DiffusionImageType* dwimage = static_cast( inputImage.GetPointer() ); mitk::DWIHeadMotionCorrectionFilter::Pointer corrfilter = mitk::DWIHeadMotionCorrectionFilter::New(); corrfilter->SetInput( dwimage ); corrfilter->SetAverageUnweighted(false); corrfilter->Update(); try { mitk::IOUtil::SaveBaseData(corrfilter->GetOutput(), argv[2]); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Catched exception: " << e.what(); mitkThrow() << "Failed with exception from subprocess!"; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp index 93aa88088e..96f1695eb4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp @@ -1,103 +1,103 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "mitkIOUtil.h" #include #include "mitkDWIHeadMotionCorrectionFilter.h" #include "mitkNrrdDiffusionImageWriter.h" typedef short DiffusionPixelType; typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; int mitkExtractSingleShellTest( int argc, char* argv[] ) { MITK_TEST_BEGIN("mitkExtractSingleShellTest"); MITK_TEST_CONDITION_REQUIRED( argc > 3, "Specify input and output and the shell to be extracted"); /* 1. Get input data */ mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); DiffusionImageType* dwimage = static_cast( inputImage.GetPointer() ); MITK_TEST_CONDITION_REQUIRED( dwimage != NULL, "Input is a dw-image"); unsigned int extract_value = 0; std::istringstream input(argv[3]); input >> extract_value; typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter FilterType; typedef DiffusionImageType::BValueMap BValueMap; // GetShellSelection from GUI BValueMap shellSelectionMap; - BValueMap originalShellMap = dwimage->GetB_ValueMap(); + BValueMap originalShellMap = dwimage->GetBValueMap(); std::vector newNumGradientDirections; shellSelectionMap[extract_value] = originalShellMap[extract_value]; newNumGradientDirections.push_back( originalShellMap[extract_value].size() ) ; DiffusionImageType::GradientDirectionContainerType::Pointer gradientContainer = dwimage->GetDirections(); FilterType::Pointer filter = FilterType::New(); filter->SetInput(dwimage->GetVectorImage()); filter->SetOriginalGradientDirections(gradientContainer); filter->SetNumGradientDirections(newNumGradientDirections); filter->SetOriginalBValueMap(originalShellMap); filter->SetShellSelectionBValueMap(shellSelectionMap); try { filter->Update(); } catch( const itk::ExceptionObject& e) { MITK_TEST_FAILED_MSG( << "Failed due to ITK exception: " << e.what() ); } DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( filter->GetOutput() ); - image->SetB_Value(dwimage->GetB_Value()); + image->SetReferenceBValue(dwimage->GetReferenceBValue()); image->SetDirections(filter->GetGradientDirections()); image->SetMeasurementFrame(dwimage->GetMeasurementFrame()); image->InitializeFromVectorImage(); /* * 3. Write output data **/ mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::Pointer dwiwriter = mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::New(); dwiwriter->SetInput( image ); dwiwriter->SetFileName( argv[2] ); try { dwiwriter->Update(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Catched exception: " << e.what(); mitkThrow() << "Failed with exception from subprocess!"; } MITK_TEST_END(); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingUserManual.dox b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingUserManual.dox index 4c37b51b46..769b664537 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingUserManual.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkDiffusionImagingUserManual.dox @@ -1,133 +1,145 @@ /** \page org_mitk_gui_qt_diffusionimaging MITK Diffusion Imaging (MITK-DI) \tableofcontents This module provides means to diffusion weighted image reconstruction, visualization and quantification. Diffusion tensors as well as different q-ball reconstruction schemes are supported. Q-ball imaging aims at recovering more detailed information about the orientations of fibers from diffusion MRI measurements and, in particular, to resolve the orientations of crossing fibers. \section QmitkDiffusionImagingUserManualIssues Known Issues \li Dicom Import: The dicom import has so far only been implemented for Siemens dicom images. MITK-DI is capable of reading the nrrd format, which is documented elsewhere [1, 2]. These files can be created by combining the raw image data with a corresponding textual header file. The file extension should be changed from *.nrrd to *.dwi or from *.nhdr to *.hdwi respectively in order to let MITK-DI recognize the diffusion related header information provided in the files. \section QmitkDiffusionImagingUserManualPreprocessing Preprocessing The preprocessing view gives an overview over the important features of a diffusion weighted image like the number of gradient directions, b-value and the measurement frame. Additionally it allows the extraction of the B0 image, reduction of gradient directions and the generation of a binary brain mask. The image volume can be modified by applying a new mesurement frame, which is useful if the measurement frame is not set correctly in the image header, or by averaging redundant gradient directions. \imageMacro{prepro1.png,"Preprocessing",9.97} +\section QmitkDiffusionImagingUserManualDWIRegistration DWI Registration + +The DWI registration view allows head-motion and eddy-current correction. First the b=0 images (if more than 1) are registered and averaged. The averaged image then serves as fixed image for the alignment of the remaining (weighted, b>0) images. + +\imageMacro{registration_basic.png,"Basic DWI Registration Interface",8.00} + +In the basic settings, a selected DW image in the Data Manager is marked as Input Data and the Start Head Motion Correction button is enabled. If more than one DW images are selected, they will be processed in a consecutive manner. For larger sets of dw images to be processed, the Advanced settings can be used. + +\imageMacro{registration_batch.png,"Batch Processing DWI Registration Interface",8.00} + +Here all valid .dwi data located in the given input folder will be processed through the head motion correction and the output will be stored in the current working directory or in a custom folder if specified in the user interface. The output image will have the _MC postfix attached to the input's name. + \section QmitkDiffusionImagingUserManualTensorReconstruction Tensor Reconstruction The tensor reconstruction view allows ITK based tensor reconstruction [3]. The advanced settings for ITK reconstruction let you configure a manual threshold on the non-diffusion weighted image. All voxels below this threshold will not be reconstructed and left blank. It is also possible to check for negative eigenvalues. The according voxels are also left blank. -\imageMacro{tensor1.png,"ITK tensor reconstruction",9.97} +\imageMacro{tensor1.png,"ITK tensor reconstruction", 9.97} A few seconds (depending on the image size) after the reconstruction button is hit, a colored image should appear in the main window. \imageMacro{tensor4.png,"Tensor image after reconstruction",16.00} To assess the quality of the tensor fit it has been proposed to calculate the model residual [9]. This calculates the residual between the measured signal and the signal predicted by the model. Large residuals indicate an inadequacy of the model or the presence of artefacts in the signal intensity (noise, head motion, etc.). To use this option: Select a DWI dataset, estimate a tensor, select both the DWI node and the tensor node in the datamanager and press Residual Image Calculation. MITK-Diffusion can show the residual for every voxel averaged over all volumes or (in the plot widget) summarized per volume or for every slice in every volume. Clicking in the widget where the residual is shown per slice will automatically let the cross-hair jump to that position in the DWI dataset. If Percentage of outliers is checked, the per volume plot will show the percentage of outliers per volume. Otherwise it will show the mean together with the first and third quantile of residuals. See [9] for more information. \imageMacro{residuals.png,"The residual widget",16.00} The view also allows the generation of artificial diffusion weighted or Q-Ball images from the selected tensor image. The ODFs of the Q-Ball image are directly initialized from the tensor values and afterwards normalized. The diffusion weighted image is estimated using the l2-norm image of the tensor image as B0. The gradient images are afterwards generated using the standard tensor equation. \section QmitkDiffusionImagingUserManualQBallReconstruction Q-Ball Reconstruction The q-ball reonstruction view implements a variety of reconstruction methods. The different reconstruction methods are described in the following: \li Numerical: The original, numerical q-ball reconstruction presented by Tuch et al. [5] \li Standard (SH): Descoteaux's reconstruction based on spherical harmonic basis functions [6] \li Solid Angle (SH): Aganj's reconstruction with solid angle consideration [7] \li ADC-profile only: The ADC-profile reconstructed with spherical harmonic basis functions \li Raw signal only: The raw signal reconstructed with spherical harmonic basis functions \imageMacro{qballs1.png,"The q-ball resonstruction view",9.99} B0 threshold works the same as in tensor reconstruction. The maximum l-level configures the size of the spherical harmonics basis. Larger l-values (e.g. l=8) allow higher levels of detail, lower levels are more stable against noise (e.g. l=4). Lambda is a regularisation parameter. Set it to 0 for no regularisation. lambda = 0.006 has proven to be a stable choice under various settings. \imageMacro{qballs2.png,"Advanced q-ball reconstruction settings",8.00} This is how a q-ball image should initially look after reconstruction. Standard q-balls feature a relatively low GFA and thus appear rather dark. Adjust the level-window to solve this. \imageMacro{qballs3.png,"q-ball image after reconstruction",16.00} \section QmitkDiffusionImagingUserManualDicomImport Dicom Import The dicom import does not cover all hardware manufacturers but only Siemens dicom images. MITK-DI is also capable of reading the nrrd format, which is documented elsewhere [1, 2]. These files can be created by combining the raw image data with a corresponding textual header file. The file extension should be changed from *.nrrd to *.dwi or from *.nhdr to *.hdwi respectively in order to let MITK-DI recognize the diffusion related header information provided in the files. In case your dicom images are readable by MITK-DI, select one or more input dicom folders and click import. Each input folder must only contain DICOM-images that can be combined into one vector-valued 3D output volume. Different patients must be loaded from different input-folders. The folders must not contain other acquisitions (e.g. T1,T2,localizer). In case many imports are performed at once, it is recommended to set the the optional output folder argument. This prevents the images from being kept in memory. \imageMacro{dicom1.png,"Dicom import",9.59} The option "Average duplicate gradients" accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "blur radius" > 0 is configured. \section QmitkDiffusionImagingUserManualFslImport FSL Import FSL diffusion data can be imported with MITK Diffusion. FSL diffusion datasets consist of 3 files: a nifty file (filename.nii.gz or filename.nii), a bvecs file (filename.bvecs), which is a text file containing the gradient vectors, and a bvals file (filename.bvecs), containing the b-values. Due to the system that selects suitable file readers, MITK will not recognize these files as diffusion datasets. In order to make MITK recognize it as diffusion, the extension must be changed from .nii.gz to .fslgz (so the new name is filename.fslgz) or from filename.nii to filename.fsl. The bvecs and bvals files have to be renamed as well(to filename.fsl.bvecs/filenames.fsl.bvecs or to filename.fslgz.bvecs/filename.fslgz.bvals). MITK can also save diffusion weighted images in FSL format. To do this the extension of the new file should be changed to .fsl or .fslgz upon saving the file. \imageMacro{fslsave.png,"Save a dwi dataset as fsl",16.00} \section QmitkDiffusionImagingUserManualQuantification Quantification The quantification view allows the derivation of different scalar anisotropy measures for the reconstructed tensors (Fractional Anisotropy, Relative Anisotropy, Axial Diffusivity, Radial Diffusivity) or q-balls (Generalized Fractional Anisotropy). \imageMacro{quantification.png,"Anisotropy quantification",2} \section QmitkDiffusionImagingUserManualVisualizationSettings ODF Visualization Setting In this small view, the visualization of ODFs and diffusion images can be configured. Depending on the selected image in the data storage, different options are shown here. For tensor or q-ball images, the visibility of glyphs in the different render windows (T)ransversal, (S)agittal, and (C)oronal can be configured here. The maximal number of glyphs to display can also be configured here for. This is usefull to keep the system response time during rendering feasible. The other options configure normalization and scaling of the glyphs. In diffusion images, a slider lets you choose the desired image channel from the vector of images (each gradient direction one image) for rendering. Furthermore reinit can be performed and texture interpolation toggled. This is how a visualization with activated glyphs should look like: \imageMacro{visualization3.png,"Q-ball image with ODF glyph visibility toggled ON",16.00} \section QmitkDiffusionImagingUserManualReferences References 1. http://teem.sourceforge.net/nrrd/format.html 2. http://www.cmake.org/Wiki/Getting_Started_with_the_NRRD_Format 3. C.F.Westin, S.E.Maier, H.Mamata, A.Nabavi, F.A.Jolesz, R.Kikinis, "Processing and visualization for Diffusion tensor MRI", Medical image Analysis, 2002, pp 93-108 5. Tuch, D.S., 2004. Q-ball imaging. Magn Reson Med 52, 1358-1372. 6. Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2007. Regularized, fast, and robust analytical Q-ball imaging. Magn Reson Med 58, 497-510. 7. Aganj, I., Lenglet, C., Sapiro, G., 2009. ODF reconstruction in q-ball imaging with solid angle consideration. Proceedings of the Sixth IEEE International Symposium on Biomedical Imaging Boston, MA. 8. Goh, A., Lenglet, C., Thompson, P.M., Vidal, R., 2009. Estimating Orientation Distribution Functions with Probability Density Constraints and Spatial Regularity. Med Image Comput Comput Assist Interv Int Conf Med Image Comput Comput Assist Interv LNCS 5761, 877 ff. 9. J.-D. Tournier, S. Mori, A. Leemans., 2011. Diffusion Tensor Imaging and Beyond. Magn Reson Med 65, 1532-1556. \section QmitkDiffusionImagingUserManualTechnicalDetail Technical Information for Developers The diffusion imaging module uses additional properties beside the ones in use in other modules, for further information see \ref DiffusionImagingPropertiesPage . \section QmitkDiffusionImagingUserManualSubManuals Manuals of componentes The MITK Diffusion tools consist of further components, which have their own documentation, see: \li \subpage org_mitk_views_fiberprocessing \li \subpage org_mitk_views_gibbstracking \li \subpage org_mitk_views_odfdetails \li \subpage org_mitk_views_partialvolumeanalysisview \li \subpage org_mitk_views_screenshotmaker \li \subpage org_mitk_views_stochasticfibertracking \li \subpage org_mitk_views_ivim \li \subpage org_mitk_diffusionimagingapp_perspectives_connectomics \li \subpage org_mitk_views_tractbasedspatialstatistics \li \subpage org_mitk_views_fiberextraction \li \subpage org_mitk_views_fiberprocessing \li \subpage org_mitk_views_odfmaximaextraction \li \subpage org_mitk_views_streamlinetracking \li \subpage org_mitk_views_fiberfoxview \li \subpage org_mitk_views_fieldmapgenerator \li \subpage org_mitk_views_denoisingview */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkFiberfoxViewUserManual.dox b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkFiberfoxViewUserManual.dox index db22724bc4..994e4b943c 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkFiberfoxViewUserManual.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkFiberfoxViewUserManual.dox @@ -1,115 +1,115 @@ /** \page org_mitk_views_fiberfoxview Fiberfox This view provides the user interface for Fiberfox [1,2,3], an interactive simulation tool for defining artificial white matter fibers and generating corresponding diffusion weighted images. Arbitrary fiber configurations like bent, crossing, kissing, twisting, and fanning bundles can be intuitively defined by positioning only a few 3D waypoints to trigger the automated generation of synthetic fibers. From these fibers a diffusion weighted signal is simulated using a flexible combination of various diffusion models. It can be modified using specified acquisition settings such as gradient direction, b-value, signal-to-noise ratio, image size, and resolution. Additionally it enables the simulation of magnetic resonance artifacts including thermal noise, Gibbs ringing, N/2 ghosting, susceptibility distortions and motion artifacts. The employed parameters can be saved and loaded as xml file with the ending ".ffp" (Fiberfox parameters). Available sections: - \ref QmitkFiberfoxViewUserManualFiberDefinition - \ref QmitkFiberfoxViewUserManualSignalGeneration - \ref QmitkFiberfoxViewUserManualKnownIssues - \ref QmitkFiberfoxViewUserManualReferences \imageMacro{Fiberfox.png, "Fig. 1: Screenshot of the Fiberfox framework. The four render windows display an axial (top left)\, sagittal (top right) and coronal (bottom left) 2D cut as well as a 3D view of a synthetic fiber helix and the fiducials used to define its shape. In the 2D views the helix is superimposing the baseline volume of the corresponding diffusion weighted image. The sagittal render window shows a close-up view on one of the circular fiducials.",16} \section QmitkFiberfoxViewUserManualFiberDefinition Fiber Definition Fiber strands are defined simply by placing markers in a 3D image volume. The fibers are then interpolated between these fiducials. Example: \li Chose an image volume to place the markers used to define the fiber pathway. If you don't have such an image available switch to the "Signal Generation" tab, define the size and spacing of the desired image and click "Generate Image". If no fiber bundle is selected, this will generate a dummy image that can be used to place the fiducials. \li Start placing fiducials at the desired positions to define the fiber pathway. To do that, click on the button with the circle pictogram, then click at the desired position and plane in the image volume and drag your mouse while keeping the button pressed to generate a circular shape. Adjust the shape using the control points (Fig. 2). The position of control point D introduces a twist of the fibers between two successive fiducials. The actual fiber generation is triggered automatically as soon as you place the second control point. \li In some cases the fibers are entangled in a way that can't be resolved by introducing an additional fiber twist. Fiberfox tries to avoid these situations, which arise from different normal orientations of succeeding fiducials, automatically. In rare cases this is not successful. Use the double-arrow button to flip the fiber positions of the selected fiducial in one dimension. Either the problem is resolved now or you can resolve it manually by adjusting the twist-control point. \li To create non elliptical fiber profile shapes switch to the Fiber Extraction View. This view provides tools to extract subesets of fibers from fiber bundles and enables to cut out arbitrary polygonal fiber shapes from existing bundles. \imageMacro{Fiberfox-Fiducial.png, "Fig. 2: Control points defining the actual shape of the fiducial. A specifies the fiducials position in space\, B and C the two ellipse radii and D the twisting angle between two successive fiducials.",10} Fiber Options: \li Real Time Fibers: If checked, each parameter adjustment (fiducial position, number of fibers, ...) will be directly applied to the selected fiber bundle. If unchecked, the fibers will only be generated if the corresponding button "Generate Fibers" is clicked. \li Advanced Options: Show/hide advanced options -\li \#Fibers: Specifies the number of fibers that will be generated for the selected bundle. +\li \# Fibers: Specifies the number of fibers that will be generated for the selected bundle. \li Fiber Sampling: Adjusts the distenace of the fiber sampling points (in mm). A higher sampling rate is needed if high curvatures are modeled. \li Tension, Continuity, Bias: Parameters controlling the shape of the splines interpolation the fiducials. See Wikipedia for details. Fiducial Options: \li Use Constant Fiducial Radius: If checked, all fiducials are treated as circles with the same radius. The first fiducial of the bundle defines the radius of all other fiducials. \li Align with grid: Click to shift all fiducial center points to the next voxel center. Operations: \li Rotation: Define the rotation of the selected fiber bundle around each axis (in degree). \li Translation: Define the translation of the selected fiber bundle along each axis (in mm). \li Scaling: Define a scaling factor for the selected fiber bundle in each dimension. \li Transform Selection: Apply specified rotation, translation and scaling to the selected Bundle/Fiducial \li Copy Bundles: Add copies of the selected fiber bundles to the datamanager. \li Join Bundles: Add new bundle to the datamanager that contains all fibers from the selected bundles. \li Include Fiducials: If checked, the specified transformation is also applied to the fiducials belonging to the selected fiber bundle and the fiducials are also copied. \imageMacro{FiberfoxExamples.png, "Fig. 3: Examples of artificial crossing (a\,b)\, fanning (c\,d)\, highly curved (e\,f)\, kissing (g\,h) and twisting (i\,j) fibers as well as of the corresponding tensor images generated with Fiberfox.",6} \section QmitkFiberfoxViewUserManualSignalGeneration Signal Generation To generate an artificial signal from the input fibers we follow the concepts recently presented by Panagiotaki et al. in a review and taxonomy of different compartment models: a flexible model combining multiple compartments is used to simulate the anisotropic diffusion inside (intra-axonal compartment) and between axons (inter-axonal compartment), isotropic diffusion outside of the axons (extra-axonal compartment 1) and the restricted diffusion in other cell types (extra-axonal compartment 2) weighted according to their respective volume fraction. A diffusion weighted image is generated from the fibers by selecting the according fiber bundle in the datamanager and clicking "Generate Image". If some other diffusion weighted image is selected together with the fiber bundle, Fiberfox directly uses the parameters of the selected image (size, spacing, gradient directions, b-values) for the signal generation process. Additionally a binary image can be selected that defines the tissue area. Voxels outside of this mask will contain no signal, only noise. Basic Image Settings: \li Image Dimensions: Specifies actual image size (number of voxels in each dimension). \li Image Spacing: Specifies voxel size in mm. Beware that changing the voxel size also changes the signal strength, e.g. increasing the resolution from 2x2x2 mm to 1x1x1 mm decreases the signal obtained for each voxel by a factor 8. \li Gradient Directions: Number of gradients directions distributed equally over the half sphere. 10% baseline images are automatically added. \li b-Value: Diffusion weighting in s/mm². If an existing diffusion weighted image is used to set the basic parameters, the b-value is defined by the gradient direction magnitudes of this image, which also enables the use of multiple b-values. Advanced Image Settings (activate checkbox "Advanced Options"): \li Repetitions: Specifies the number of averages used for the acquisition to reduce noise. \li Signal Scale: Additional scaling factor for the signal in each voxel. The default value of 125 results in a maximum signal amplitude of 1000 for 2x2x2 mm voxels. Beware that changing this value without changing the noise variance results in a changed SNR. Adjustment of this value might be needed if the overall signal values are much too high or much too low (depends on a variety of factors like voxel size and relaxation times). \li Echo Time TE: Time between the 90° excitation pulse and the first spin echo. Increasing this time results in a stronger T2-relaxation effect (Wikipedia). \li Line Readout Time: Time to read one line in k-space. Increasing this time results in a stronger T2* effect which causes an attenuation of the higher frequencies in phase direction (here along y-axis) which again results in a blurring effect of sharp edges perpendicular to the phase direction. \li Tinhom Relaxation: Time constant specifying the signal decay due to magnetic field inhomogeneities (also called T2'). Together with the tissue specific relaxation time constant T2 this defines the T2* decay constant: T2*=(T2 T2')/(T2+T2') \li Fiber Radius (in µm): Used to calculate the volume fractions of the used compartments (fiber, water, etc.). If set to 0 (default) the fiber radius is set automatically so that the voxel containing the most fibers is filled completely. A realistic axon radius ranges from about 5 to 20 microns. Using the automatic estimation the resulting value might very well be much larger or smaller than this range. \li Simulate Signal Relaxation: If checked, the relaxation induced signal decay is simulated, other wise the parameters TE, Line Readout Time, Tinhom, and T2 are ignored. \li Disable Partial Volume Effects: If checked, the actual volume fractions of the single compartments are ignored. A voxel will either be filled by the intra axonal compartment completely or will contain no fiber at all. \li Output Volume Fractions: Output a double image for each compartment. The voxel values correspond to the volume fraction of the respective compartment. Compartment Settings: The group-boxes "Intra-axonal Compartment", "Inter-axonal Compartment" and "Extra-axonal Compartments" allow the specification which model to use and the corresponding model parameters. Currently the following models are implemented: \li Stick: The “stick” model describes diffusion in an idealized cylinder with zero radius. Parameter: Diffusivity d \li Zeppelin: Cylindrically symmetric diffusion tensor. Parameters: Parallel diffusivity d|| and perpendicular diffusivity d⊥ \li Tensor: Full diffusion tensor. Parameters: Parallel diffusivity d|| and perpendicular diffusivity constants d⊥1 and d⊥2 \li Ball: Isotropic compartment. Parameter: Diffusivity d \li Astrosticks: Consists of multiple stick models pointing in different directions. The single stick orientations can either be distributed equally over the sphere or are sampled randomly. The model represents signal coming from a type of glial cell called astrocytes, or populations of axons with arbitrary orientation. Parameters: randomization of the stick orientations and diffusivity of the sticks d. \li Dot: Isotropically restricted compartment. No parameter. For a detailed description of the single models, please refer to Panagiotaki et al. "Compartment models of the diffusion MR signal in brain white matter: A taxonomy and comparison". Additionally to the model parameters, each compartment has its own T2 signal relaxation constant (in ms). Noise and Artifacts: \li Noise: Add Rician or Chi-Square distributed noise with the specified variance to the signal. \li Spikes: Add signal spikes to the k-space signal resulting in stripe artifacts across the corresponding image slice. \li Aliasing: Aliasing artifacts occur if the FOV in phase direction is smaller than the imaged object. The parameter defines the percentage by which the FOV is shrunk. \li N/2 Ghosts: Specify the offset between successive lines in k-space. This offset causes ghost images in distance N/2 in phase direction due to the alternating EPI readout directions. \li Distortions: Simulate distortions due to magnetic field inhomogeneities. This is achieved by adding an additional phase during the readout process. The input is a frequency map specifying the inhomogeneities. The "Fieldmap Generator" view provides an interface to generate simple artificial frequency maps. \li Motion Artifacts: To simulate motion artifacts, the fiber configuration is moved between the signal simulation of the individual gradient volumes. The motion can be performed randomly, where the parameters are used to define the +/- maximum of the corresponding motion, or linearly, where the parameters define the maximum rotation/translation around/along the corresponding axis at the and of the simulated acquisition. \li Eddy Currents: EXPERIMENTAL! This feature is currently being tested and might not yet behave as expected! \li Gibbs Ringing: Ringing artifacts occurring on edges in the image due to the frequency low-pass filtering caused by the limited size of the k-space. \section QmitkFiberfoxViewUserManualKnownIssues Known Issues \li If fiducials are created in one of the marginal slices of the underlying image, a position change of the fiducial can be observed upon selection/deselection. If the fiducial is created in any other slice this bug does not occur. \li If a scaling factor is applied to the selcted fiber bundle, the corresponding fiducials are not scaled accordingly. \li In some cases the automatic update of the selected fiber bundle is not triggered even if "Real Time Fibers" is checked, e.g. if a fiducial is deleted. If this happens on can always force an update by pressing the "Generate Fibers" button. If any other issues or feature requests arises during the use of Fiberfox, please don't hesitate to send us an e-mail or directly report the issue in our bugtracker: http://bugs.mitk.org/ \section QmitkFiberfoxViewUserManualReferences References [1] Peter F. Neher, Frederik B. Laun, Bram Stieltjes, and Klaus H. Fritzsche: Fiberfox: Facilitating the creation of realistic white matter software phantoms, Magn Reson Med, DOI: 10.1002/mrm.25045. [2] Peter F. Neher, Frederik B. Laun, Bram Stieltjes, and Klaus H. Fritzsche: Fiberfox: An extensible system for generating realistic white matter software phantoms, MICCAI CDMRI Workshop, Nagoya; 09/2013 [3] Peter F. Neher, Bram Stieltjes, Frederik B. Laun, Hans-Peter Meinzer, and Klaus H. Fritzsche: Fiberfox: Fiberfox: A novel tool to generate software phantoms of complex fiber geometries, ISMRM, Salt Lake City; 04/2013 */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox index fcf0ff18c5..4bd7ea5194 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/QmitkOdfMaximaExtractionViewUserManual.dox @@ -1,49 +1,49 @@ /** \page org_mitk_views_odfmaximaextraction Peak Extraction View This view provides the user interface to extract the peaks of tensors and the spherical harmonic representation of Q-Balls. Available sections: - \ref OdfMaxUserManualInputData - \ref OdfMaxUserManualOutputData - \ref OdfMaxUserManualMethods - \ref OdfMaxUserManualParameters - \ref OdfMaxUserManualReferences \section OdfMaxUserManualInputData Input Data Mandatory Input: \li DTI image or image containing the spherical harmonic coefficients. The SH coefficient images can be obtain from the Q-Ball reconstruction view by enabling the checkbox in the advanced options. Optional Input: \li Binary mask to define the extraction area. \section OdfMaxUserManualOutputData Output Data \li Vector field: 3D representation of the resulting peaks. Only for visualization purposes (the peaks are scaled additionally to the specified normalization to improve the visualization)! -\li \#Directions per Voxel: Image containing the number of extracted peaks per voxel as image value. +\li \# Directions per Voxel: Image containing the number of extracted peaks per voxel as image value. \li Direction Images: One image for each of the extracted peaks per voxel. Each voxel contains one direction vector as image value. Use this output for evaluation purposes of the extracted peaks. \section OdfMaxUserManualMethods Peak Extraction Methods \li If a tensor image is used as input, the output is simply the largest eigenvector of each voxel. \li The finite differences extraction uses a higly sampled version of the image ODFs, extracts all local maxima and clusters the resulting directions that point in a similar direction. \li For details about the analytical method (experimental) please refer to [1]. \imageMacro{crossingmaxima.png,"Peaks of a fiber crossing extracted using finite differences method.",10.12} \section OdfMaxUserManualParameters Input Parameters \li Vector normalization method (no normalization, maximum normalization of the vecors of one voxel and independent normalization of each vecor). \li SH Order: Specify the order of the spherical harmonic coefficients. \li Maximum number of peaks to extract. If more peaks are found only the largest are kept. \li Threshold to discard small peaks. Value relative to the largest peak of the respective voxel. \li Absolute threshold on the peak size. To evaluate this threshold the peaks are additionally weighted by their GFA (low GFA voxels are more likely to be discarded). This threshold is only used for the finite differences extraction method. \section OdfMaxUserManualReferences References [1] Aganj et al. Proceedings of the Thirteenth International Conference on Medical Image Computing and Computer Assisted Intervention 2010 */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/registration_basic.png b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/registration_basic.png new file mode 100644 index 0000000000..457dabe392 Binary files /dev/null and b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/registration_basic.png differ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/registration_batch.png b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/registration_batch.png new file mode 100644 index 0000000000..0489a1626e Binary files /dev/null and b/Plugins/org.mitk.gui.qt.diffusionimaging/documentation/UserManual/registration_batch.png differ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.cpp index ca4c053c18..dd721c7515 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.cpp @@ -1,365 +1,489 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //misc #define _USE_MATH_DEFINES #include // Blueberry #include #include // Qmitk #include "QmitkDiffusionRegistrationView.h" #include // MITK #include #include #include #include #include // Qt #include #include #include +#include #include - - - -/* -#include "QmitkDataStorageComboBox.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -*/ - -#include +#include #define _USE_MATH_DEFINES #include - - QmitkRegistrationWorker::QmitkRegistrationWorker(QmitkDiffusionRegistrationView* view) : m_View(view) { } void QmitkRegistrationWorker::run() { typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::BValueMap BValueMap; - for(unsigned int i=0; i< m_View->m_SelectedDiffusionNodes.size(); i++) + unsigned int totalImagesCount; + if( !m_View->m_IsBatch ) { + totalImagesCount = m_View->m_SelectedDiffusionNodes.size(); + } + else + { + totalImagesCount = m_View->m_BatchList.size(); + } + m_View->m_TotalFiles = totalImagesCount; + + QString inputPath = m_View->m_Controls->m_InputFolderTextbox->text(); + QString outputPath = m_View->m_Controls->m_OutputFolderTextbox->text(); + + + for(unsigned int i=0; i< totalImagesCount; i++) + { + if(m_View->m_IsAborted){ + m_View->m_RegistrationThread.quit(); + return; + } + m_View->m_CurrentFile = i+1; m_View->m_GlobalRegisterer = QmitkDiffusionRegistrationView::DWIHeadMotionCorrectionFilterType::New(); - mitk::DataNode::Pointer node = m_View->m_SelectedDiffusionNodes.at(i); - DiffusionImageType::Pointer inImage = dynamic_cast*>(node->GetData()); + //mitk::DataNode::Pointer node = m_View->m_SelectedDiffusionNodes.at(i); + DiffusionImageType::Pointer inImage; + mitk::DataNode::Pointer node; + + if( !m_View->m_IsBatch ) + { + node = m_View->m_SelectedDiffusionNodes.at(i); + inImage = dynamic_cast*>(node->GetData()); + } + else + { + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( m_View->m_BatchList.at(i).toStdString() ); + inImage = static_cast( inputImage.GetPointer() ); + } + + if(inImage.IsNull()) { MITK_ERROR << "Error occured: can't get input image. \nAborting"; return; } m_View->m_GlobalRegisterer->SetInput(inImage); - try{ + try + { m_View->m_GlobalRegisterer->Update(); } catch( mitk::Exception e ) { MITK_ERROR << "Internal error occured: " << e.what() << "\nAborting"; } if( m_View->m_GlobalRegisterer->GetIsInValidState() ) { - DiffusionImageType::Pointer image = m_View->m_GlobalRegisterer->GetOutput(); - mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); - imageNode->SetData( image ); - QString name = node->GetName().c_str(); - imageNode->SetName((name+"_MC").toStdString().c_str()); - m_View->GetDataStorage()->Add(imageNode); - } + if(! m_View->m_IsBatch) + { + DiffusionImageType::Pointer image = m_View->m_GlobalRegisterer->GetOutput(); + mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); + imageNode->SetData( image ); + QString name = node->GetName().c_str(); + imageNode->SetName((name+"_MC").toStdString().c_str()); + m_View->GetDataStorage()->Add(imageNode); + } + else + { + mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::Pointer dwiwriter = + mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::New(); + + dwiwriter->SetInput( m_View->m_GlobalRegisterer->GetOutput() ); + + QString name = m_View->m_BatchList.at(i); + name = name.replace(".dwi", "_MC.dwi", Qt::CaseInsensitive); + name = name.replace(inputPath, outputPath, Qt::CaseInsensitive); + dwiwriter->SetFileName( name.toStdString().c_str() ); + + try + { + dwiwriter->Update(); + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Catched exception: " << e.what(); + mitkThrow() << "Failed with exception from subprocess!"; + } + } + } } - - m_View->m_RegistrationThread.quit(); } const std::string QmitkDiffusionRegistrationView::VIEW_ID = "org.mitk.views.diffusionregistrationview"; QmitkDiffusionRegistrationView::QmitkDiffusionRegistrationView() : QmitkAbstractView() , m_Controls( 0 ) , m_DiffusionImage( NULL ) , m_ThreadIsRunning(false) , m_Steps(100) , m_LastStep(0) , m_GlobalRegisterer(NULL) , m_RegistrationWorker(this) { m_RegistrationWorker.moveToThread(&m_RegistrationThread); connect(&m_RegistrationThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_RegistrationThread, SIGNAL(started()), &m_RegistrationWorker, SLOT(run())); connect(&m_RegistrationThread, SIGNAL(finished()), this, SLOT(AfterThread())); connect(&m_RegistrationThread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_RegistrationTimer = new QTimer(this); } // Destructor QmitkDiffusionRegistrationView::~QmitkDiffusionRegistrationView() { delete m_RegistrationTimer; } // update Registration status and generate fiber bundle void QmitkDiffusionRegistrationView::TimerUpdate() { int currentStep = m_GlobalRegisterer->GetCurrentStep(); mitk::ProgressBar::GetInstance()->Progress(currentStep-m_LastStep); UpdateRegistrationStatus(); m_LastStep = currentStep; } // update gui elements after registration is finished void QmitkDiffusionRegistrationView::AfterThread() { m_ThreadIsRunning = false; m_RegistrationTimer->stop(); mitk::ProgressBar::GetInstance()->Progress(m_GlobalRegisterer->GetSteps()-m_LastStep+1); UpdateGUI(); if( !m_GlobalRegisterer->GetIsInValidState() ) { QMessageBox::critical( NULL, "Registration", "An internal error occured, or user canceled the Registration.\n Please check the log for details." ); return; } UpdateRegistrationStatus(); m_GlobalRegisterer = 0; } // start Registration timer and update gui elements before Registration is started void QmitkDiffusionRegistrationView::BeforeThread() { m_ThreadIsRunning = true; m_RegistrationTime = QTime::currentTime(); m_ElapsedTime = 0; m_RegistrationTimer->start(1000); m_LastStep = 0; UpdateGUI(); } void QmitkDiffusionRegistrationView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkDiffusionRegistrationViewControls; m_Controls->setupUi( parent ); AdvancedSettings(); connect( m_RegistrationTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->m_RegistrationStopButton, SIGNAL(clicked()), this, SLOT(StopRegistration()) ); connect( m_Controls->m_RegistrationStartButton, SIGNAL(clicked()), this, SLOT(StartRegistration()) ); connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); + connect( m_Controls->m_SelectInputButton, SIGNAL(clicked()), this, SLOT(AddInputFolderName()) ); + connect( m_Controls->m_SelectOutputButton, SIGNAL(clicked()), this, SLOT(AddOutputFolderName()) ); + connect( m_Controls->m_StartBatchButton, SIGNAL(clicked()), this, SLOT(StartBatch()) ); + this->m_Parent = parent; } } // show/hide advanced settings frame void QmitkDiffusionRegistrationView::AdvancedSettings() { - //m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); - m_Controls->m_AdvancedFrame->setVisible(false); - m_Controls->m_AdvancedSettingsCheckbox->setVisible(false); + m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); } void QmitkDiffusionRegistrationView::OnSelectionChanged( berry::IWorkbenchPart::Pointer, const QList& nodes ) { if (m_ThreadIsRunning) return; bool foundDwiVolume = false; QString tempSelectedNames = ""; m_DiffusionImage = NULL; m_SelectedDiffusionNodes.clear(); // iterate selection for( int i=0; i*>(node->GetData()) ) { foundDwiVolume = true; m_SelectedDiffusionNodes.push_back(node); if(m_SelectedDiffusionNodes.size() > 0){tempSelectedNames += "\n";} tempSelectedNames += (node->GetName().c_str()); } } m_Controls->m_RegistrationStartButton->setEnabled(foundDwiVolume); if (foundDwiVolume) { m_Controls->m_DiffusionImageLabel->setText(tempSelectedNames); m_Controls->m_InputData->setTitle("Input Data"); } else { m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_Controls->m_InputData->setTitle("Please Select Input Data"); } UpdateGUI(); } // update gui elements displaying Registrations status void QmitkDiffusionRegistrationView::UpdateRegistrationStatus() { if (m_GlobalRegisterer.IsNull()) return; m_ElapsedTime += m_RegistrationTime.elapsed()/1000; m_RegistrationTime.restart(); unsigned long hours = m_ElapsedTime/3600; unsigned long minutes = (m_ElapsedTime%3600)/60; unsigned long seconds = m_ElapsedTime%60; m_Controls->m_RegistrationTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); m_Controls->m_CurrentStepLabel->setText( QString::number((int)(100*(float)(m_GlobalRegisterer->GetCurrentStep()-1)/m_GlobalRegisterer->GetSteps()))+"%" ); - + m_Controls->m_CurrentFileLabel->setText( QString::number(m_CurrentFile)+" / "+QString::number(m_TotalFiles) ); } void QmitkDiffusionRegistrationView::UpdateGUI() { if (!m_ThreadIsRunning && (m_SelectedDiffusionNodes.size() > 0) ) { m_Controls->m_RegistrationStopButton->setEnabled(false); m_Controls->m_RegistrationStartButton->setEnabled(true); + m_Controls->m_StartBatchButton->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_RegistrationStopButton->setText("Stop"); m_Controls->m_RegistrationStartButton->setToolTip("Start Registration"); m_Controls->m_RegistrationStopButton->setToolTip(""); } else if (!m_ThreadIsRunning) { m_Controls->m_RegistrationStopButton->setEnabled(false); m_Controls->m_RegistrationStartButton->setEnabled(false); + m_Controls->m_StartBatchButton->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_RegistrationStopButton->setText("Stop"); m_Controls->m_RegistrationStartButton->setToolTip("No Diffusion image selected."); m_Controls->m_RegistrationStopButton->setToolTip(""); } else { m_Controls->m_RegistrationStopButton->setEnabled(true); m_Controls->m_RegistrationStartButton->setEnabled(false); + m_Controls->m_StartBatchButton->setEnabled(false); m_Controls->m_AdvancedFrame->setEnabled(false); - m_Controls->m_AdvancedFrame->setVisible(false); - m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); m_Controls->m_RegistrationStartButton->setToolTip("Registration in progress."); m_Controls->m_RegistrationStopButton->setToolTip("Cancel Registration"); } + } void QmitkDiffusionRegistrationView::SetFocus() { m_Controls->m_RegistrationStartButton->setFocus(); } void QmitkDiffusionRegistrationView::StartRegistration() { if(m_ThreadIsRunning) { MITK_WARN("QmitkDiffusionRegistrationView")<<"Thread already running!"; return; } m_GlobalRegisterer = NULL; if (m_SelectedDiffusionNodes.size()<1) { QMessageBox::information( NULL, "Warning", "Please load and select a diffusion image before starting image processing."); return; } + m_IsBatch = false; + m_IsAborted = false; m_Controls->m_RegistrationStartButton->setEnabled(false); + m_Controls->m_StartBatchButton->setEnabled(false); // start worker thread m_RegistrationThread.start(QThread::NormalPriority); return; } void QmitkDiffusionRegistrationView::StopRegistration() { if (m_GlobalRegisterer.IsNull()) return; + m_IsAborted = true; m_GlobalRegisterer->SetAbortRegistration(true); m_Controls->m_RegistrationStopButton->setEnabled(false); m_Controls->m_RegistrationStopButton->setText("Stopping ..."); return; } + +void QmitkDiffusionRegistrationView::AddInputFolderName() +{ + + // SELECT FOLDER DIALOG + QFileDialog* w = new QFileDialog( m_Parent, QString("Select the input folder with DWI files within") ); + w->setFileMode( QFileDialog::Directory ); + + // RETRIEVE SELECTION + if ( w->exec() != QDialog::Accepted ) + return; + + m_Controls->m_InputFolderTextbox->setText(w->selectedFiles()[0]); +} + +void QmitkDiffusionRegistrationView::AddOutputFolderName() +{ + // SELECT FOLDER DIALOG + QFileDialog* w = new QFileDialog( m_Parent, QString("Select the output folder") ); + w->setFileMode( QFileDialog::Directory ); + + // RETRIEVE SELECTION + if ( w->exec() != QDialog::Accepted ) + return; + + m_Controls->m_OutputFolderTextbox->setText(w->selectedFiles()[0]); +} + +void QmitkDiffusionRegistrationView::StartBatch() +{ + QString inputPath = m_Controls->m_InputFolderTextbox->text(); + QString outputPath = m_Controls->m_OutputFolderTextbox->text(); + + if(inputPath == outputPath){ + QMessageBox::information( NULL, "Error", "Input and Output folders can't be the same"); + return; + } + + QStringList list, filters; + filters<<"*.dwi"; + QDirIterator dirIterator(inputPath, + filters, + QDir::Files|QDir::NoSymLinks); + + while (dirIterator.hasNext()) + { + dirIterator.next(); + + list.append(dirIterator.fileInfo().absoluteFilePath()); + std::cout << dirIterator.fileInfo().absoluteFilePath().toStdString() << endl; + m_BatchList = list; + + m_IsBatch = true; + m_IsAborted = false; + + if(m_ThreadIsRunning) + { + MITK_WARN("QmitkDiffusionRegistrationView")<<"Thread already running!"; + return; + } + m_GlobalRegisterer = NULL; + + if (m_BatchList.size()<1) + { + QMessageBox::information( NULL, "Error", "No diffusion images were found in the selected input folder."); + return; + } + + m_Controls->m_RegistrationStartButton->setEnabled(false); + m_Controls->m_StartBatchButton->setEnabled(false); + + // start worker thread + m_RegistrationThread.start(QThread::NormalPriority); + } +} diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.h index 48ed41c49e..e21ca6bb75 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationView.h @@ -1,124 +1,135 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include "ui_QmitkDiffusionRegistrationViewControls.h" #include "mitkDiffusionImage.h" #include #include #include typedef short DiffusionPixelType; /*! \brief View for diffusion image registration / head motion correction \sa QmitkFunctionality \ingroup Functionalities */ // Forward Qt class declarations using namespace std; class QmitkDiffusionRegistrationView; class QmitkRegistrationWorker : public QObject { Q_OBJECT public: QmitkRegistrationWorker(QmitkDiffusionRegistrationView* view); public slots: void run(); private: QmitkDiffusionRegistrationView* m_View; }; class QmitkDiffusionRegistrationView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const string VIEW_ID; QmitkDiffusionRegistrationView(); virtual ~QmitkDiffusionRegistrationView(); virtual void CreateQtPartControl(QWidget *parent); void SetFocus(); typedef mitk::DWIHeadMotionCorrectionFilter< DiffusionPixelType > DWIHeadMotionCorrectionFilterType; protected slots: void StartRegistration(); void StopRegistration(); void AfterThread(); ///< update gui etc. after registrations has finished void BeforeThread(); ///< start timer etc. void TimerUpdate(); + void AddInputFolderName(); + void AddOutputFolderName(); + void StartBatch(); + void AdvancedSettings(); protected: /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer, const QList&); Ui::QmitkDiffusionRegistrationViewControls* m_Controls; mitk::DiffusionImage::Pointer m_DiffusionImage; std::vector< mitk::DataNode::Pointer > m_SelectedDiffusionNodes; private: void UpdateRegistrationStatus(); ///< update textual status display of the Registration process void UpdateGUI(); ///< update button activity etc. dpending on current datamanager selection /** flags etc. */ + bool m_IsBatch, m_IsAborted; + QStringList m_BatchList; bool m_ThreadIsRunning; QTimer* m_RegistrationTimer; QTime m_RegistrationTime; unsigned long m_ElapsedTime; unsigned long m_Steps; int m_LastStep; + unsigned int m_CurrentFile; + unsigned int m_TotalFiles; + + // the Qt parent of our GUI (NOT of this object) + QWidget* m_Parent; /** global Registerer and friends */ itk::SmartPointer m_GlobalRegisterer; QmitkRegistrationWorker m_RegistrationWorker; QThread m_RegistrationThread; friend class QmitkRegistrationWorker; }; diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationViewControls.ui index 91659dd4aa..e68125a3a4 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDiffusionRegistrationViewControls.ui @@ -1,226 +1,298 @@ QmitkDiffusionRegistrationViewControls 0 0 435 744 Form 0 0 DWI Head Motion Correction Qt::Vertical 20 40 false Start DWI registration/Head Motion Correction Start Head Motion Correction Please Select Input Data Diffusion Image <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Monitor - + + + + 0 + 0 + + - Process Time: + Process Time: - + Progress: - + - - + - + + + + File: + + + + + + + - + + + false Qt::LeftToRight Stop :/QtWidgetsExt/stop.xpm:/QtWidgetsExt/stop.xpm + + + + Advanced + + + true QFrame::StyledPanel QFrame::Raised 9 0 9 0 4 - + - Advanced settings + Batch mode + + + + + + + This will process all the .dwi files inside the input folder + + + + + + + + + + + + + + 0 + 0 + + + + START BATCH PROCESSING + + + + + + + + 0 + 0 + + + + Select input folder + + + + + + + + 0 + 0 + + + + Select output folder - - - - Advanced Settings - - - diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp index fb4fccac05..e07728c7ca 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp @@ -1,971 +1,970 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkTensorReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include #include #include #include #include // itk includes #include "itkTimeProbe.h" //#include "itkTensor.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "mitkTeemDiffusionTensor3DReconstructionImageFilter.h" #include "itkDiffusionTensor3DReconstructionImageFilter.h" #include "itkTensorImageToDiffusionImageFilter.h" #include "itkPointShell.h" #include "itkVector.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkTensorReconstructionWithEigenvalueCorrectionFilter.h" //#include "itkFreeWaterEliminationFilter.h" #include "mitkProperties.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "mitkDiffusionImageMapper.h" #include "mitkLookupTableProperty.h" #include "mitkLookupTable.h" #include "mitkImageStatisticsHolder.h" #include #include #include #include const std::string QmitkTensorReconstructionView::VIEW_ID = "org.mitk.views.tensorreconstruction"; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; using namespace berry; QmitkTensorReconstructionView::QmitkTensorReconstructionView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL) { m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); } QmitkTensorReconstructionView::~QmitkTensorReconstructionView() { } void QmitkTensorReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkTensorReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); Advanced1CheckboxClicked(); } } void QmitkTensorReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkTensorReconstructionView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkTensorReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_StartReconstruction), SIGNAL(clicked()), this, SLOT(Reconstruct()) ); connect( (QObject*)(m_Controls->m_Advanced1), SIGNAL(clicked()), this, SLOT(Advanced1CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_TensorsToDWIButton), SIGNAL(clicked()), this, SLOT(TensorsToDWI()) ); connect( (QObject*)(m_Controls->m_TensorsToQbiButton), SIGNAL(clicked()), this, SLOT(TensorsToQbi()) ); connect( (QObject*)(m_Controls->m_ResidualButton), SIGNAL(clicked()), this, SLOT(ResidualCalculation()) ); connect( (QObject*)(m_Controls->m_PerSliceView), SIGNAL(pointSelected(int, int)), this, SLOT(ResidualClicked(int, int)) ); } } void QmitkTensorReconstructionView::ResidualClicked(int slice, int volume) { // Use image coord to reset crosshair // Find currently selected diffusion image // Update Label // to do: This position should be modified in order to skip B0 volumes that are not taken into account // when calculating residuals // Find the diffusion image mitk::DiffusionImage* diffImage; mitk::DataNode::Pointer correctNode; mitk::BaseGeometry* geometry; if (m_DiffusionImage.IsNotNull()) { diffImage = static_cast*>(m_DiffusionImage->GetData()); geometry = diffImage->GetGeometry(); // Remember the node whose display index must be updated correctNode = mitk::DataNode::New(); correctNode = m_DiffusionImage; } if(diffImage != NULL) { typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; GradientDirectionContainerType::Pointer dirs = diffImage->GetDirections(); for(unsigned int i=0; iSize() && i<=volume; i++) { GradientDirectionType grad = dirs->ElementAt(i); // check if image is b0 weighted if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001) { volume++; } } QString pos = "Volume: "; pos.append(QString::number(volume)); pos.append(", Slice: "); pos.append(QString::number(slice)); m_Controls->m_PositionLabel->setText(pos); if(correctNode) { int oldDisplayVal; correctNode->GetIntProperty("DisplayChannel", oldDisplayVal); std::string oldVal = QString::number(oldDisplayVal).toStdString(); std::string newVal = QString::number(volume).toStdString(); correctNode->SetIntProperty("DisplayChannel",volume); correctNode->SetSelected(true); this->FirePropertyChanged("DisplayChannel", oldVal, newVal); correctNode->UpdateOutputInformation(); mitk::Point3D p3 = m_MultiWidget->GetCrossPosition(); itk::Index<3> ix; geometry->WorldToIndex(p3, ix); // ix[2] = slice; mitk::Vector3D vec; vec[0] = ix[0]; vec[1] = ix[1]; vec[2] = slice; mitk::Vector3D v3New; geometry->IndexToWorld(vec, v3New); mitk::Point3D origin = geometry->GetOrigin(); mitk::Point3D p3New; p3New[0] = v3New[0] + origin[0]; p3New[1] = v3New[1] + origin[1]; p3New[2] = v3New[2] + origin[2]; m_MultiWidget->MoveCrossToPosition(p3New); m_MultiWidget->RequestUpdate(); } } } void QmitkTensorReconstructionView::Advanced1CheckboxClicked() { bool check = m_Controls-> m_Advanced1->isChecked(); m_Controls->frame->setVisible(check); } void QmitkTensorReconstructionView::Activated() { QmitkFunctionality::Activated(); } void QmitkTensorReconstructionView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkTensorReconstructionView::ResidualCalculation() { // Extract dwi and dti from current selection // In case of multiple selections, take the first one, since taking all combinations is not meaningful mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); mitk::DiffusionImage::Pointer diffImage = mitk::DiffusionImage::New(); TensorImageType::Pointer tensorImage; std::string nodename; if(m_DiffusionImage.IsNotNull()) { diffImage = static_cast*>(m_DiffusionImage->GetData()); } else return; if(m_TensorImage.IsNotNull()) { mitk::TensorImage* mitkVol; mitkVol = static_cast(m_TensorImage->GetData()); mitk::CastToItkImage(mitkVol, tensorImage); m_TensorImage->GetStringProperty("name", nodename); } else return; typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; mitk::DiffusionImage::GradientDirectionContainerType* gradients = diffImage->GetDirections(); // Find the min and the max values from a baseline image mitk::ImageStatisticsHolder *stats = diffImage->GetStatistics(); //Initialize filter that calculates the modeled diffusion weighted signals FilterType::Pointer filter = FilterType::New(); filter->SetInput( tensorImage ); filter->SetBValue(diffImage->GetReferenceBValue()); filter->SetGradientList(gradients); filter->SetMin(stats->GetScalarValueMin()); filter->SetMax(stats->GetScalarValueMax()); filter->Update(); // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); - image->SetReferenceBValue(diffImage->GetReferenceBValue()); image->SetDirections(gradients); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_DWI"); node->SetName(newname.toAscii()); GetDefaultDataStorage()->Add(node); mitk::DiffusionImage::BValueMap map =image->GetBValueMap(); mitk::DiffusionImage::IndicesVector b0Indices = map[0]; typedef itk::ResidualImageFilter ResidualImageFilterType; ResidualImageFilterType::Pointer residualFilter = ResidualImageFilterType::New(); residualFilter->SetInput(diffImage->GetVectorImage()); residualFilter->SetSecondDiffusionImage(image->GetVectorImage()); residualFilter->SetGradients(gradients); residualFilter->SetB0Index(b0Indices[0]); residualFilter->SetB0Threshold(30); residualFilter->Update(); itk::Image::Pointer residualImage = itk::Image::New(); residualImage = residualFilter->GetOutput(); mitk::Image::Pointer mitkResImg = mitk::Image::New(); mitk::CastToMitkImage(residualImage, mitkResImg); stats = mitkResImg->GetStatistics(); float min = stats->GetScalarValueMin(); float max = stats->GetScalarValueMax(); mitk::LookupTableProperty::Pointer lutProp = mitk::LookupTableProperty::New(); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(min, max); // If you don't want to use the whole color range, you can use // SetValueRange, SetHueRange, and SetSaturationRange lookupTable->Build(); vtkSmartPointer reversedlookupTable = vtkSmartPointer::New(); reversedlookupTable->SetTableRange(min+1, max); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookupTable->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } lut->SetVtkLookupTable(lookupTable); lutProp->SetLookupTable(lut); // Create lookuptable mitk::DataNode::Pointer resNode=mitk::DataNode::New(); resNode->SetData( mitkResImg ); resNode->SetName("Residual Image"); resNode->SetProperty("LookupTable", lutProp); bool b; resNode->GetBoolProperty("use color", b); resNode->SetBoolProperty("use color", false); GetDefaultDataStorage()->Add(resNode); m_MultiWidget->RequestUpdate(); // Draw Graph std::vector means = residualFilter->GetMeans(); std::vector q1s = residualFilter->GetQ1(); std::vector q3s = residualFilter->GetQ3(); std::vector percentagesOfOUtliers = residualFilter->GetPercentagesOfOutliers(); m_Controls->m_ResidualAnalysis->SetMeans(means); m_Controls->m_ResidualAnalysis->SetQ1(q1s); m_Controls->m_ResidualAnalysis->SetQ3(q3s); m_Controls->m_ResidualAnalysis->SetPercentagesOfOutliers(percentagesOfOUtliers); if(m_Controls->m_PercentagesOfOutliers->isChecked()) { m_Controls->m_ResidualAnalysis->DrawPercentagesOfOutliers(); } else { m_Controls->m_ResidualAnalysis->DrawMeans(); } // Draw Graph for volumes per slice in the QGraphicsView std::vector< std::vector > outliersPerSlice = residualFilter->GetOutliersPerSlice(); int xSize = outliersPerSlice.size(); if(xSize == 0) { return; } int ySize = outliersPerSlice[0].size(); // Find maximum in outliersPerSlice double maxOutlier= 0.0; for(int i=0; imaxOutlier) { maxOutlier = outliersPerSlice[i][j]; } } } // Create some QImage QImage qImage(xSize, ySize, QImage::Format_RGB32); QImage legend(1, 256, QImage::Format_RGB32); QRgb value; vtkSmartPointer lookup = vtkSmartPointer::New(); lookup->SetTableRange(0.0, maxOutlier); lookup->Build(); reversedlookupTable->SetTableRange(0, maxOutlier); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookup->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } // Fill qImage for(int i=0; iMapValue(out); int r, g, b; r = _rgba[0]; g = _rgba[1]; b = _rgba[2]; value = qRgb(r, g, b); qImage.setPixel(i,j,value); } } for(int i=0; i<256; i++) { double* rgba = lookup->GetTableValue(i); int r, g, b; r = rgba[0]*255; g = rgba[1]*255; b = rgba[2]*255; value = qRgb(r, g, b); legend.setPixel(0,255-i,value); } QString upper = QString::number(maxOutlier, 'g', 3); upper.append(" %"); QString lower = QString::number(0.0); lower.append(" %"); m_Controls->m_UpperLabel->setText(upper); m_Controls->m_LowerLabel->setText(lower); QGraphicsScene* scene = new QGraphicsScene; QGraphicsScene* scene2 = new QGraphicsScene; QPixmap pixmap(QPixmap::fromImage(qImage)); QGraphicsPixmapItem *item = new QGraphicsPixmapItem( pixmap, 0, scene); item->scale(10.0, 3.0); QPixmap pixmap2(QPixmap::fromImage(legend)); QGraphicsPixmapItem *item2 = new QGraphicsPixmapItem( pixmap2, 0, scene2); item2->scale(20.0, 1.0); m_Controls->m_PerSliceView->SetResidualPixmapItem(item); m_Controls->m_PerSliceView->setScene(scene); m_Controls->m_LegendView->setScene(scene2); m_Controls->m_PerSliceView->show(); m_Controls->m_PerSliceView->repaint(); m_Controls->m_LegendView->setHorizontalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->setVerticalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->show(); m_Controls->m_LegendView->repaint(); } void QmitkTensorReconstructionView::Reconstruct() { int method = m_Controls->m_ReconctructionMethodBox->currentIndex(); switch (method) { case 0: ItkTensorReconstruction(m_DiffusionImages); break; case 1: TensorReconstructionWithCorr(m_DiffusionImages); break; default: ItkTensorReconstruction(m_DiffusionImages); } } void QmitkTensorReconstructionView::TensorReconstructionWithCorr (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); while ( itemiter != itemiterend ) // for all items { typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientDirectionContainerType; DiffusionImageType* vols = static_cast((*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // TENSOR RECONSTRUCTION MITK_INFO << "Tensor reconstruction with correction for negative eigenvalues"; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::TensorReconstructionWithEigenvalueCorrectionFilter< DiffusionPixelType, TTensorPixelType > ReconstructionFilter; float b0Threshold = m_Controls->m_TensorReconstructionThreshold->value(); GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = vols->GetDirections()->Begin(); it != vols->GetDirections()->End(); it++) { gradientContainerCopy->push_back(it.Value()); } ReconstructionFilter::Pointer reconFilter = ReconstructionFilter::New(); reconFilter->SetGradientImage( gradientContainerCopy, vols->GetVectorImage() ); reconFilter->SetBValue(vols->GetReferenceBValue()); reconFilter->SetB0Threshold(b0Threshold); reconFilter->Update(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer outputTensorImg = reconFilter->GetOutput(); typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(outputTensorImg, outputTensorImg->GetRequestedRegion()); tensorIt.GoToBegin(); int negatives = 0; while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iInitializeByItk( outputTensorImg.GetPointer() ); image->SetVolume( outputTensorImg->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); SetDefaultNodeProperties(node, nodename+"_EigenvalueCorrected_DT"); GetDefaultDataStorage()->Add(node, *itemiter); mitk::ProgressBar::GetInstance()->Progress(); ++itemiter; } mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); } } void QmitkTensorReconstructionView::ItkTensorReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // TENSOR RECONSTRUCTION clock.Start(); MITK_DEBUG << "Tensor reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::DiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientDirectionContainerType; GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = vols->GetDirections()->Begin(); it != vols->GetDirections()->End(); it++) { gradientContainerCopy->push_back(it.Value()); } tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, vols->GetVectorImage() ); tensorReconstructionFilter->SetBValue(vols->GetReferenceBValue()); tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreshold->value() ); tensorReconstructionFilter->Update(); clock.Stop(); MITK_DEBUG << "took " << clock.GetMean() << "s."; // TENSORS TO DATATREE mitk::TensorImage::Pointer image = mitk::TensorImage::New(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer tensorImage; tensorImage = tensorReconstructionFilter->GetOutput(); // Check the tensor for negative eigenvalues if(m_Controls->m_CheckNegativeEigenvalues->isChecked()) { typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetRequestedRegion()); tensorIt.GoToBegin(); while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; //typedef itk::Tensor TensorType2; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iSetDirection( vols->GetVectorImage()->GetDirection() ); image->InitializeByItk( tensorImage.GetPointer() ); image->SetVolume( tensorReconstructionFilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); SetDefaultNodeProperties(node, nodename+"_WeightedLinearLeastSquares_DT"); GetDefaultDataStorage()->Add(node, *itemiter); mitk::ProgressBar::GetInstance()->Progress(); ++itemiter; } mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return; } } void QmitkTensorReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name) { node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) ); node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) ); node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New()); node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New()); node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2)); node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1)); node->SetProperty( "visible", mitk::BoolProperty::New( true ) ); node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) ); node->SetProperty ("layer", mitk::IntProperty::New(100)); node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) ); node->SetProperty( "name", mitk::StringProperty::New(name) ); } void QmitkTensorReconstructionView::TensorsToDWI() { DoTensorsToDWI(m_TensorImages); } void QmitkTensorReconstructionView::TensorsToQbi() { for (unsigned int i=0; isize(); i++) { mitk::DataNode::Pointer tensorImageNode = m_TensorImages->at(i); MITK_INFO << "starting Q-Ball estimation"; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(dynamic_cast(tensorImageNode->GetData()), itkvol); typedef itk::TensorImageToQBallImageFilter< TTensorPixelType, TTensorPixelType > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->Update(); typedef itk::Vector OutputPixelType; typedef itk::Image OutputImageType; mitk::QBallImage::Pointer image = mitk::QBallImage::New(); OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName(tensorImageNode->GetName()+"_Qball"); GetDefaultDataStorage()->Add(node, tensorImageNode); } } void QmitkTensorReconstructionView::OnSelectionChanged( std::vector nodes ) { m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); bool foundDwiVolume = false; bool foundTensorVolume = false; m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_DiffusionImage = NULL; m_TensorImage = NULL; m_Controls->m_InputData->setTitle("Please Select Input Data"); // iterate selection for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if (node.IsNull()) continue; // only look at interesting types if(dynamic_cast*>(node->GetData())) { foundDwiVolume = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_DiffusionImages->push_back(node); m_DiffusionImage = node; } else if(dynamic_cast(node->GetData())) { foundTensorVolume = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_TensorImages->push_back(node); m_TensorImage = node; } } m_Controls->m_StartReconstruction->setEnabled(foundDwiVolume); m_Controls->m_TensorsToDWIButton->setEnabled(foundTensorVolume); m_Controls->m_TensorsToQbiButton->setEnabled(foundTensorVolume); if (foundDwiVolume || foundTensorVolume) m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_ResidualButton->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PercentagesOfOutliers->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PerSliceView->setEnabled(foundDwiVolume && foundTensorVolume); } template itk::VectorContainer >::Pointer QmitkTensorReconstructionView::MakeGradientList() { itk::VectorContainer >::Pointer retval = itk::VectorContainer >::New(); vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); for(int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval->push_back(v); } // Add 0 vector for B0 vnl_vector_fixed v; v.fill(0.0); retval->push_back(v); return retval; } void QmitkTensorReconstructionView::DoTensorsToDWI(mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); while ( itemiter != itemiterend ) // for all items { std::string nodename; (*itemiter)->GetStringProperty("name", nodename); mitk::TensorImage* vol = static_cast((*itemiter)->GetData()); typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(vol, itkvol); typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; FilterType::GradientListPointerType gradientList = FilterType::GradientListType::New(); switch(m_Controls->m_TensorsToDWINumDirsSelect->currentIndex()) { case 0: gradientList = MakeGradientList<12>(); break; case 1: gradientList = MakeGradientList<42>(); break; case 2: gradientList = MakeGradientList<92>(); break; case 3: gradientList = MakeGradientList<162>(); break; case 4: gradientList = MakeGradientList<252>(); break; case 5: gradientList = MakeGradientList<362>(); break; case 6: gradientList = MakeGradientList<492>(); break; case 7: gradientList = MakeGradientList<642>(); break; case 8: gradientList = MakeGradientList<812>(); break; case 9: gradientList = MakeGradientList<1002>(); break; default: gradientList = MakeGradientList<92>(); } double bVal = m_Controls->m_TensorsToDWIBValueEdit->text().toDouble(); // DWI ESTIMATION clock.Start(); MBI_INFO << "DWI Estimation "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "DWI Estimation for %s", nodename.c_str()).toAscii()); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->SetBValue(bVal); filter->SetGradientList(gradientList); //filter->SetNumberOfThreads(1); filter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMean() << "s."; // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetReferenceBValue(bVal); image->SetDirections(gradientList); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); node->SetName(nodename+"_DWI"); GetDefaultDataStorage()->Add(node, *itemiter); mitk::ProgressBar::GetInstance()->Progress(); ++itemiter; } mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "DWI estimation failed:", ex.GetDescription()); return ; } }