diff --git a/CMakeLists.txt b/CMakeLists.txt index 7dcf12f..d48386d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,49 +1,69 @@ -cmake_minimum_required(VERSION 2.8.9) +cmake_minimum_required(VERSION 3.5) #----------------------------------------------------------------------------- -if(NOT Slicer_SOURCE_DIR) - set(EXTENSION_NAME MABMIS) - set(EXTENSION_HOMEPAGE "http://wiki.slicer.org/slicerWiki/index.php/Documentation/Nightly/Extensions/MABMIS") - set(EXTENSION_CATEGORY "Segmentation") - set(EXTENSION_CONTRIBUTORS "Xiaofeng Liu (GE GLobal Research), Minjeong Kim (UNC), Dinggang Shen (UNC), Jim Miller (GE GLobal Research)") - set(EXTENSION_ICONURL "http://wiki.slicer.org/slicerWiki/images/e/e2/MABMIS_Icon.png") - set(EXTENSION_DESCRIPTION "Multi-Atlas Based Group Segmentation") - set(EXTENSION_SCREENSHOTURLS "http://wiki.slicer.org/slicerWiki/images/2/2c/MABMIS_trainning_GUI.png http://wiki.slicer.org/slicerWiki/images/6/64/MABMIS_testing_GUI.png http://wiki.slicer.org/slicerWiki/images/c/c3/MABMIS_algorithm.png") -endif() +project(MABMIS) +set(EXTENSION_HOMEPAGE "http://wiki.slicer.org/slicerWiki/index.php/Documentation/Nightly/Extensions/MABMIS") +set(EXTENSION_CATEGORY "Segmentation") +set(EXTENSION_CONTRIBUTORS "Xiaofeng Liu (GE Global Research), Minjeong Kim (UNC), Dinggang Shen (UNC), Jim Miller (GE Global Research)") +set(EXTENSION_ICONURL "http://wiki.slicer.org/slicerWiki/images/e/e2/MABMIS_Icon.png") +set(EXTENSION_DESCRIPTION "Multi-Atlas Based Group Segmentation") +set(EXTENSION_SCREENSHOTURLS "http://wiki.slicer.org/slicerWiki/images/2/2c/MABMIS_trainning_GUI.png http://wiki.slicer.org/slicerWiki/images/6/64/MABMIS_testing_GUI.png http://wiki.slicer.org/slicerWiki/images/c/c3/MABMIS_algorithm.png") + +find_package(Slicer REQUIRED) +include(${Slicer_USE_FILE}) + +set(MABMIS_ITK_COMPONENTS + ITKIOImageBase + ITKTransform + ITKIOTransformBase + ITKIOXML + ITKImageGrid + ITKImageFunction + ITKPDEDeformableRegistration + ITKSmoothing + ) +find_package(ITK 4.6 COMPONENTS ${MABMIS_ITK_COMPONENTS} REQUIRED) +set(ITK_NO_IO_FACTORY_REGISTER_MANAGER 1) # See Libs/ITKFactoryRegistration/CMakeLists.txt +include(${ITK_USE_FILE}) -#----------------------------------------------------------------------------- -if(NOT Slicer_SOURCE_DIR) - find_package(Slicer REQUIRED) - include(${Slicer_USE_FILE}) -endif() #----------------------------------------------------------------------------- set(MODULE_TARGET_LIBRARIES ${ITK_LIBRARIES} ) - -#----------------------------------------------------------------------------- -set(MODULE_NAME IGR3D_MABMIS_Training) +message(STATUS ${MODULE_TARGET_LIBRARIES}) #----------------------------------------------------------------------------- SEMMacroBuildCLI( - NAME ${MODULE_NAME} + NAME IGR3D_MABMIS_Training TARGET_LIBRARIES ${MODULE_TARGET_LIBRARIES} ADDITIONAL_SRCS itkMABMISAtlasXMLFile.cxx + Training.cxx + commonMABMIS.cxx #EXECUTABLE_ONLY ) - message(STATUS ${ITK_LIBRARIES}) #----------------------------------------------------------------------------- -set(MODULE_NAME IGR3D_MABMIS_Testing) +SEMMacroBuildCLI( + NAME IGR3D_MABMIS_Testing + TARGET_LIBRARIES ${MODULE_TARGET_LIBRARIES} + ADDITIONAL_SRCS + itkMABMISAtlasXMLFile.cxx + Testing.cxx + commonMABMIS.cxx + #EXECUTABLE_ONLY + ) #----------------------------------------------------------------------------- SEMMacroBuildCLI( - NAME ${MODULE_NAME} + NAME IGR3D_MABMIS_Direct_Invoke TARGET_LIBRARIES ${MODULE_TARGET_LIBRARIES} ADDITIONAL_SRCS itkMABMISAtlasXMLFile.cxx + Testing.cxx + Training.cxx + commonMABMIS.cxx #EXECUTABLE_ONLY ) @@ -54,11 +74,5 @@ if(BUILD_TESTING) endif() #----------------------------------------------------------------------------- -if(NOT Slicer_SOURCE_DIR) - include(${Slicer_EXTENSION_CPACK}) -endif() - - - - -################################################ +include(${Slicer_EXTENSION_GENERATE_CONFIG}) +include(${Slicer_EXTENSION_CPACK}) diff --git a/Data/Input/TestData/TestData.xml b/Data/Input/TestData/TestData.xml index 2567297..c933de8 100644 --- a/Data/Input/TestData/TestData.xml +++ b/Data/Input/TestData/TestData.xml @@ -4,10 +4,10 @@ 4 - - - - + + + + diff --git a/IGR3D_MABMIS_Direct_Invoke.cxx b/IGR3D_MABMIS_Direct_Invoke.cxx new file mode 100644 index 0000000..fd30f24 --- /dev/null +++ b/IGR3D_MABMIS_Direct_Invoke.cxx @@ -0,0 +1,468 @@ +#include "IGR3D_MABMIS_Direct_InvokeCLP.h" + +#include "Testing.h" +#include "Training.h" + +#include "itkAffineTransform.h" +#include "itkTransformFileReader.h" +#include "itkResampleImageFilter.h" +#include "itkWindowedSincInterpolateImageFunction.h" +#include "itkBSplineInterpolateImageFunction.h" +#include "itkNearestNeighborInterpolateImageFunction.h" +#include "itkLabelImageGaussianInterpolateImageFunction.h" +#include "itkMinimumMaximumImageCalculator.h" +#include "itkPluginUtilities.h" + +const unsigned int Dimension = 3; +typedef itk::Image ShortImageType; +typedef itk::AffineTransform AffineType; +typedef itk::Transform TransformType; +AffineType::Pointer identity = AffineType::New(); + +// interpolationQuality: +// -1 = linear (default and fall-back) +// 0 = nearest neighbor +// 1 = LabelImageGaussianInterpolate (smoother alternative to nearest neighbor) +// 3 = cubic BSpline +// 5 = WindowedSinc with radius of 5 voxels +template +void resampleAndWrite(const std::string & inFile, const std::string & outFile, + ShortImageType::RegionType region, + ShortImageType::PointType origin, + ShortImageType::SpacingType spacing, + ShortImageType::DirectionType direction, + TransformType::Pointer transform, + int interpolationQuality) +{ + typedef itk::Image ImageType; + + typedef itk::ImageFileReader ReaderType; + typename ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(inFile); + reader->Update(); + + typedef itk::MinimumMaximumImageCalculator MinMaxType; + typename MinMaxType::Pointer minMax = MinMaxType::New(); + minMax->SetImage(reader->GetOutput()); + minMax->ComputeMinimum(); + + typedef itk::ResampleImageFilter ResampleFilterType; + typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New(); + switch (interpolationQuality) + { + case 0: + resampleFilter->SetInterpolator(itk::NearestNeighborInterpolateImageFunction::New()); + break; + case 1: + resampleFilter->SetInterpolator(itk::LabelImageGaussianInterpolateImageFunction::New()); + break; + case 3: + resampleFilter->SetInterpolator(itk::BSplineInterpolateImageFunction::New()); //order 3 (=cubic) by default + break; + case 5: + resampleFilter->SetInterpolator(itk::WindowedSincInterpolateImageFunction::New()); + break; + case -1: + default: + // ResampleImageFilter uses linear by default + break; + } + resampleFilter->SetInput(reader->GetOutput()); + resampleFilter->SetSize(region.GetSize()); + resampleFilter->SetOutputOrigin(origin); + resampleFilter->SetOutputSpacing(spacing); + resampleFilter->SetOutputDirection(direction); + resampleFilter->SetTransform(transform); + resampleFilter->SetDefaultPixelValue(minMax->GetMinimum()); + resampleFilter->Update(); + + typedef itk::ImageFileWriter WriterType; + typename WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(outFile); + writer->SetInput(resampleFilter->GetOutput()); + writer->SetUseCompression(true); + writer->Update(); +} + +//does the template-based dispatch +void resampleAndWrite(const std::string & inFile, const std::string & outFile, + ShortImageType::RegionType region, + ShortImageType::PointType origin, + ShortImageType::SpacingType spacing, + ShortImageType::DirectionType direction, + TransformType::Pointer transform, + int interpolationQuality, + itk::ImageIOBase::IOComponentType componentType) +{ +switch (componentType) + { + case itk::ImageIOBase::UCHAR: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::CHAR: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::USHORT: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::SHORT: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::UINT: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::INT: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::ULONG: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::LONG: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + //case itk::ImageIOBase::ULONGLONG: + // resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + // break; + //case itk::ImageIOBase::LONGLONG: + // resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + // break; + case itk::ImageIOBase::FLOAT: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::DOUBLE: + resampleAndWrite(inFile, outFile, region, origin, spacing, direction, transform, interpolationQuality); + break; + case itk::ImageIOBase::UNKNOWNCOMPONENTTYPE: + default: + itkGenericExceptionMacro("Error: unknown component type in image " << inFile); + break; + } +} + +// resample inFile into outFile so it matches the image grid of refFile +// uses the pixel type of the inFile +void resampleOrCopy(const std::string & inFile, + const std::string & outFile, + const std::string & refFile, + TransformType::Pointer transform, + int interpolationQuality) +{ + typedef itk::ImageFileReader ShortReaderType; + ShortReaderType::Pointer shortReader = ShortReaderType::New(); + shortReader->SetFileName(refFile); + shortReader->UpdateOutputInformation(); + ShortImageType::Pointer image = shortReader->GetOutput(); + + //reference image metadata + ShortImageType::RegionType region = image->GetLargestPossibleRegion(); + ShortImageType::PointType origin = image->GetOrigin(); + ShortImageType::SpacingType spacing = image->GetSpacing(); + ShortImageType::DirectionType direction = image->GetDirection(); + + //inImage componenet type + itk::ImageIOBase::IOPixelType pixelType; + itk::ImageIOBase::IOComponentType componentType; + itk::GetImageType(inFile, pixelType, componentType); + + shortReader->SetFileName(inFile); + shortReader->UpdateOutputInformation(); + image = shortReader->GetOutput(); + + bool transformIsIdentity = false; + AffineType * aTransform = dynamic_cast(transform.GetPointer()); + if (aTransform && aTransform->Metric(identity) == 0.0) + { + transformIsIdentity = true; + } + + if (region != image->GetLargestPossibleRegion() + || origin != image->GetOrigin() + || spacing != image->GetSpacing() + || direction != image->GetDirection() + || transformIsIdentity != true + || GetExtension(inFile) != GetExtension(outFile)) + { + std::cout << "Resampling " << inFile << std::endl; + std::cout << " into " << outFile << std::endl; + resampleAndWrite(inFile, outFile, + region, origin, spacing, direction, transform, + interpolationQuality, componentType); + } + else + { + std::cout << "Copying " << inFile << std::endl; + std::cout << " to " << outFile << std::endl; + itksys::SystemTools::CopyFileAlways(inFile, outFile); + } +} + +//transform the images to match the atlas root and write them +void resampleOrCopyN(const std::string & rootFilename, const std::string & outDir, + const std::vector & inFiles, + const std::vector & outFiles, + std::vector & transforms, + int interpolationQuality) +{ + for (unsigned i = 0; i < inFiles.size(); i++) + { + resampleOrCopy(inFiles[i], outDir + outFiles[i], rootFilename, transforms[i], interpolationQuality); + } +} + +#define HANDLE_IMAGE(n) \ +{ \ + if (image##n##Arg.isSet()) \ + { \ + imageFileNames.push_back(image##n); \ + if (transform##n##Arg.isSet()) \ + { \ + transformFileNames.push_back(transform##n); \ + } \ + else \ + { \ + transformFileNames.push_back(""); \ + } \ + if (segmentation##n##Arg.isSet()) \ + { \ + segmentationFileNames.push_back(segmentation##n); \ + } \ + else \ + { \ + segmentationFileNames.push_back(""); \ + } \ + } \ +} +//segmentationFileNames + +int main( int argc, char *argv[] ) +{ + PARSE_ARGS; + + int segmentationInterpolationQuality = 1; + if (segmentationInterpolationArg.isSet() && segmentationInterpolation == "Nearest") + { + segmentationInterpolationQuality = 0; + } + + int imageInterpolationQuality = 3; + if (imageInterpolationArg.isSet() && imageInterpolation == "Linear") + { + imageInterpolationQuality = -1; + } + if (imageInterpolationArg.isSet() && imageInterpolation == "WindowedSinc") + { + imageInterpolationQuality = 5; + } + + if (mode == "Create imageXML" && !imageListXMLArg.isSet()) + { + std::cerr << "Create imageXML mode requires imageListXML argument to be set!" << std::endl; + return EXIT_FAILURE; + } + + imageListXML = ReplacePathSepForOS(imageListXML); + atlasTreeXML = ReplacePathSepForOS(atlasTreeXML); + imageDir = ReplacePathSepForOS(imageDir); + + try + { + itk::MABMISImageData miData; + const std::string extension = ".nrrd"; + std::string listDir = itksys::SystemTools::GetParentDirectory(imageListXML); + listDir = itksys::SystemTools::GetRealPath(listDir); + imageDir = itksys::SystemTools::GetRealPath(imageDir); + if (!imageDirArg.isSet() || imageDir.empty() || imageDir == ".") + { + imageDir = listDir; + } + imageDir += '/'; + + std::vector transforms; + std::vector imageFileNames; + std::vector segmentationFileNames; + itk::MABMISAtlas* atlas = nullptr; + std::string rootFilename; + if (atlasTreeXMLArg.isSet() && itksys::SystemTools::FileExists(atlasTreeXML)) + { + //read atlas tree + itk::MABMISAtlasXMLFileReader::Pointer atlasReader = + itk::MABMISAtlasXMLFileReader::New(); + atlasReader->SetFilename(atlasTreeXML); + atlasReader->GenerateOutputInformation(); + atlas = atlasReader->GetOutputObject(); + + //get atlas root image + for (unsigned i = 0; i < atlas->m_NumberAllAtlases; i++) + { + if (atlas->m_Tree[i] == atlas->m_TreeRoot) + { + rootFilename = atlas->m_AtlasFilenames[i]; + break; + } + } + + const std::string& atlDir = atlas->m_AtlasDirectory; + if (atlDir.size() > 2 && atlDir[0] == '.' && (atlDir[1] == '\\' || atlDir[1] == '/')) + { + listDir = itksys::SystemTools::GetParentDirectory(atlasTreeXML) + + atlas->m_AtlasDirectory.substr(1) + '/'; + } + else + { + listDir = atlas->m_AtlasDirectory + '/'; + } + atlas->m_AtlasDirectory = listDir; + rootFilename = listDir + rootFilename; + + if (mode == "Retrain atlas") + { + std::cout << "Re-training the atlas: " << atlasTreeXML << std::endl; + imageDir = listDir; //new images should be put into the old atlas directory too + for (unsigned i = 0; i < atlas->m_NumberAllAtlases - atlas->m_NumberSimulatedAtlases; i++) + { + miData.m_ImageFileNames.push_back(atlas->m_AtlasFilenames[i]); + imageFileNames.push_back(imageDir + atlas->m_AtlasFilenames[i]); + miData.m_SegmentationFileNames.push_back(atlas->m_AtlasSegmentationFilenames[i]); + segmentationFileNames.push_back(imageDir + atlas->m_AtlasSegmentationFilenames[i]); + transforms.push_back(identity.GetPointer()); + } + } + } + else + { + if (mode == "Direct invoke" || mode == "Retrain atlas") + { + itkGenericExceptionMacro("Atlas XML file is missing " << atlasTreeXML); + } + else if (mode == "Train atlas") + { + std::cout << "Initial creation of the atlas: " << atlasTreeXML << std::endl; + } + } + unsigned numberOfPreExistingImages = miData.m_ImageFileNames.size(); + std::cout << "Using " << imageDir << " as a temporary folder." << std::endl; + itksys::SystemTools::MakeDirectory(imageDir); + + //check which parameters are present + std::vector transformFileNames; + HANDLE_IMAGE(0); + HANDLE_IMAGE(1); + HANDLE_IMAGE(2); + HANDLE_IMAGE(3); + HANDLE_IMAGE(4); + HANDLE_IMAGE(5); + HANDLE_IMAGE(6); + HANDLE_IMAGE(7); + HANDLE_IMAGE(8); + HANDLE_IMAGE(9); + + //fill miData and read array of transforms + itk::TransformFileReader::Pointer transformReader = itk::TransformFileReader::New(); + for (unsigned i = numberOfPreExistingImages; i < imageFileNames.size(); i++) + { + miData.m_ImageFileNames.push_back("image" + std::to_string(i) + extension); + miData.m_SegmentationFileNames.push_back("image" + std::to_string(i) + std::string("-label") + extension); + if (transformFileNames[i - numberOfPreExistingImages].empty()) + { + transforms.push_back(identity.GetPointer()); + } + else + { + transformReader->SetFileName(transformFileNames[i - numberOfPreExistingImages]); + transformReader->Update(); + if (transformReader->GetTransformList()->size() > 1) + { + itkGenericExceptionMacro("Only simple transforms are supported. " + << transformFileNames[i] << " contains more than one transform!"); + } + itk::TransformBaseTemplate::Pointer genericTransform = transformReader->GetTransformList()->front(); + TransformType::Pointer aTransform = static_cast(genericTransform.GetPointer()); + if (aTransform.IsNull()) + { + itkGenericExceptionMacro("Could not interpret " << transformFileNames[i] << " as a supported transform!"); + } + transforms.push_back(aTransform); + } + } + miData.m_NumberImageData = miData.m_ImageFileNames.size(); + if (miData.m_NumberImageData == 0 || (atlas && atlas->m_AtlasFilenames.size() == miData.m_NumberImageData)) + { + itkGenericExceptionMacro("Some images have to be provided"); + } + if (rootFilename.empty()) //there was no atlas + { + rootFilename = imageFileNames[0]; + } + + if (mode == "Train atlas" || mode == "Retrain atlas") + { + //check if all images have segmentations + for (unsigned i = 0; i < miData.m_SegmentationFileNames.size(); i++) + { + if (miData.m_SegmentationFileNames[i].empty()) + { + itkGenericExceptionMacro(<< miData.m_ImageFileNames[i] << " lacks a segmentation!"); + } + } + } + + resampleOrCopyN(rootFilename, imageDir, imageFileNames, miData.m_ImageFileNames, transforms, imageInterpolationQuality); + if (mode == "Create imageXML" || imageListXMLArg.isSet()) + { + itk::MABMISImageDataXMLFileWriter::Pointer miWriter = + itk::MABMISImageDataXMLFileWriter::New(); + miWriter->SetObject(&miData); + miWriter->SetFilename(imageListXML); + miWriter->WriteFile(); + } + //set path after the xml file is written, so xml uses relative paths + miData.m_DataDirectory = imageDir; + miData.m_OutputDirectory = imageDir; + + if (mode == "Direct invoke") + { + Testing(&miData, atlas, iterations, sigma); + + //inverse the transforms and resample the segmentations back into original image grids + for (unsigned i = numberOfPreExistingImages; i < segmentationFileNames.size(); i++) + { + std::string segDir = itksys::SystemTools::GetParentDirectory(segmentationFileNames[i]); + itksys::SystemTools::MakeDirectory(segDir); + + TransformType::Pointer invT = transforms[i]->GetInverseTransform(); + if (!invT) //not possible to invert + { + invT = identity; + std::cerr << "Error: it was not possible to invert tranform " << transformFileNames[i] << std::endl; + std::cerr << "Resampling segmentation back using identity transform" << std::endl; + } + + resampleOrCopy(imageDir + miData.m_SegmentationFileNames[i], segmentationFileNames[i], + imageFileNames[i], invT, segmentationInterpolationQuality); + } + } + if (mode == "Train atlas" || mode == "Retrain atlas") + { + //segmentations are now inputs which accompany the images, and need to be transformed the same way + resampleOrCopyN(rootFilename, imageDir, segmentationFileNames, miData.m_SegmentationFileNames, transforms, segmentationInterpolationQuality); + Training(&miData, atlasTreeXML, iterations, sigma); + } + } + catch (itk::ExceptionObject &exc) + { + std::cout << exc; + return EXIT_FAILURE; + } + catch (std::runtime_error &exc) + { + std::cout << "Runtime error has occurred: " << exc.what() << std::endl; + return EXIT_FAILURE; + } + catch (...) + { + std::cout << "An unknown error has occurred!" << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; //no error has occurred +} diff --git a/IGR3D_MABMIS_Direct_Invoke.xml b/IGR3D_MABMIS_Direct_Invoke.xml new file mode 100644 index 0000000..3d7d446 --- /dev/null +++ b/IGR3D_MABMIS_Direct_Invoke.xml @@ -0,0 +1,290 @@ + + + Segmentation + MABMIS Direct Invoke + Allows invoking MABMIS without pre-created XML files, as well as atlas creation and re-training. Does image resampling if needed. + + + + Dženan Zukić (Kitware Inc.) + + + + + Image list and atlas XMLs + + + m + mode + Processing mode.\n\ +\ + Direct invoke: provide atlas and test images, get segmentations for test images.\n\ +\ + Create imageXML: provide test images, get ImageXML ready for invoking Testing or Training. If atlasTreeXML is provided, test images are resampled to match its root image. Otherwise, test images are resampled to match the first provided test image.\n\ +\ + Train atlas: atlas is created using the provided images and their segmentations.\n\ +\ + Retrain atlas: atlas is expanded with additional images and their segmentations.\n + + Direct invoke + Direct invoke + Create imageXML + Train atlas + Retrain atlas + + + + atlasTreeXML + a + xml file of the multi-atlases after training + + input + + + + imageListXML + l + xml file that will contain a list of images to be segmented + + output + + + + imageDir + d + A directory where to write transformed, resampled and cropped images. Default: imageListXML's location. + + input + + + + + + Parameters used for registration + + + i + iterations + Comma separated list of iterations, for low resolution, middle resolution, and high resolution. + + 5,3,2 + + + + s + sigma + The standard deviation for smoothing deformation field + + 1.5 + + + + + + Interpolation used for re-sampling + + + imageInterpolation + Which interpolation to use for image resampling + + CubicSpline + Linear + CubicSpline + WindowedSinc + + + + segmentationInterpolation + Which interpolation to use for segmentation resampling + + Smooth + Nearest + Smooth + + + + + + Inputs and outputs + + image0 + An input image + + input + + + segmentation0 + The output segmentation + + output + + + transform0 + A transform which aligns the image to the atlas root + + input + + + image1 + An input image + + input + + + segmentation1 + The output segmentation + + output + + + transform1 + A transform which aligns the image to the atlas root + + input + + + image2 + An input image + + input + + + segmentation2 + The output segmentation + + output + + + transform2 + A transform which aligns the image to the atlas root + + input + + + image3 + An input image + + input + + + segmentation3 + The output segmentation + + output + + + transform3 + A transform which aligns the image to the atlas root + + input + + + image4 + An input image + + input + + + segmentation4 + The output segmentation + + output + + + transform4 + A transform which aligns the image to the atlas root + + input + + + image5 + An input image + + input + + + segmentation5 + The output segmentation + + output + + + transform5 + A transform which aligns the image to the atlas root + + input + + + image6 + An input image + + input + + + segmentation6 + The output segmentation + + output + + + transform6 + A transform which aligns the image to the atlas root + + input + + + image7 + An input image + + input + + + segmentation7 + The output segmentation + + output + + + transform7 + A transform which aligns the image to the atlas root + + input + + + image8 + An input image + + input + + + segmentation8 + The output segmentation + + output + + + transform8 + A transform which aligns the image to the atlas root + + input + + + image9 + An input image + + input + + + segmentation9 + The output segmentation + + output + + + transform9 + A transform which aligns the image to the atlas root + + input + + + diff --git a/IGR3D_MABMIS_Testing.cxx b/IGR3D_MABMIS_Testing.cxx index 8db7b23..2da35b5 100644 --- a/IGR3D_MABMIS_Testing.cxx +++ b/IGR3D_MABMIS_Testing.cxx @@ -8,529 +8,9 @@ PURPOSE. See the above copyright notices for more information. =========================================================================*/ -// for math -#include -#include -#include -#include -#include -#include -#include -#include -#include -// basic itk -#include "itkVector.h" - -// registration -#include "itkImageRegistrationMethod.h" -#include "itkSymmetricForcesDemonsRegistrationFilter.h" - -// interpolator -#include "itkLinearInterpolateImageFunction.h" -#include "itkNearestNeighborInterpolateImageFunction.h" - -// reader / writer -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" - -// filter -#include "itkResampleImageFilter.h" - -#include "itkCastImageFilter.h" -#include "itkWarpImageFilter.h" - -#include "itkHistogramMatchingImageFilter.h" -#include "itkAddImageFilter.h" -#include "itkDivideImageFilter.h" -#include "itkMultiplyImageFilter.h" -#include "itkWarpVectorImageFilter.h" -#include "itkInverseDeformationFieldImageFilter.h" - -// for affine transformation -#include "itkTransform.h" -#include "itkAffineTransform.h" -#include "itkImageRegistrationMethod.h" - -// for Diffeomorphic Demons -#include -#include -#include - -// including itksys::SystemTools::MakeDirectory(char*) -#include -#include - -// To include all related header files #include "IGR3D_MABMIS_TestingCLP.h" -#include "itkMABMISImageOperationFilter.h" -#include "itkMABMISDeformationFieldFilter.h" -#include "itkMABMISSimulateData.h" -#include "itkMABMISImageRegistrationFilter.h" -#include "itkMABMISTreeOperation.h" -#include "itkMABMISBasicOperationFilter.h" - -#include "itkMABMISAtlasXMLFile.h" -#include -static std::string ReplacePathSepForOS( const std::string & input ) -{ - std::string output = input; -#ifdef _WIN32 - std::replace(output.begin(), output.end(), '/', FILESEP); -#else - std::replace(output.begin(), output.end(), '\\', FILESEP); -#endif - return output; -} - -static std::string -GetExtension(const std::string & filename) -{ - std::string fileExt( itksys::SystemTools::GetFilenameLastExtension(filename) ); - //If the last extension is .gz, then need to pull off 2 extensions. - //.gz is the only valid compression extension. - if ( fileExt == std::string(".gz") ) - { - fileExt = itksys::SystemTools::GetFilenameLastExtension( - itksys::SystemTools::GetFilenameWithoutLastExtension(filename) ); - fileExt += ".gz"; - } - return ( fileExt ); -} - -static std::string -GetRootName(const std::string & filename) -{ - const std::string fileExt = GetExtension(filename); - - // Create a base filename - // i.e Image.hdr --> Image - if ( fileExt.length() > 0 //Ensure that an extension was - // found - && filename.length() > fileExt.length() //Ensure that the filename does - // not contain only the extension - ) - { - const std::string::size_type it = filename.find_last_of(fileExt); - const std::string baseName( filename, 0, it - ( fileExt.length() - 1 ) ); - return ( baseName ); - } - //Default to return same as input when the extension is nothing (Analyze) - return ( filename ); -} - -typedef double CoordinateRepType; -const unsigned int SpaceDimension = ImageDimension; - -// basic data type -typedef unsigned char CharPixelType; // for image IO usage -typedef float FloatPixelType; // for -typedef int IntPixelType; -typedef short ShortPixelType; -typedef float InternalPixelType; // for internal processing usage -typedef itk::Vector VectorPixelType; - -// basic image type -typedef itk::Image CharImageType; -typedef itk::Image IntImageType; -typedef itk::Image ShortImageType; -typedef itk::Image FloatImageType; -typedef itk::Image InternalImageType; -typedef itk::Image DeformationFieldType; - -// basic iterator type -typedef itk::ImageRegionIterator DeformationFieldIteratorType; -typedef itk::ImageRegionIterator InternalImageIteratorType; -typedef itk::ImageRegionIterator CharImageIteratorType; - -// basic image reader/writer related type -typedef itk::ImageFileReader CharImageReaderType; -typedef itk::ImageFileReader InternalImageReaderType; -typedef itk::ImageFileWriter InternalImageWriterType; - -typedef itk::WarpImageFilter InternalWarpFilterType; -typedef itk::ImageFileWriter CharImageWriterType; -typedef itk::ImageFileWriter IntImageWriterType; -typedef itk::ImageFileWriter FloatImageWriterType; -typedef itk::ImageFileWriter ShortImageWriterType; - -typedef itk::ImageFileReader DeformationFieldReaderType; -typedef itk::ImageFileWriter DeformationFieldWriterType; - -////////////////////////////////////////////////////////////////////////////// -// image filter type -typedef itk::ResampleImageFilter ResampleFilterType; -typedef itk::HistogramMatchingImageFilter InternalHistMatchFilterType; - -//////////////////////////////////////////////////////////////////////////// -// operation on deformation fields -typedef itk::WarpVectorImageFilter WarpVectorFilterType; -typedef itk::InverseDeformationFieldImageFilter - InverseDeformationFieldImageFilterType; -typedef itk::AddImageFilter AddImageFilterType; - -// global bool variables to adjust the procedure - -bool isEvaluate = false; // if false, we do not know the ground-truth of labels -bool isDebug = false; // false;//true; // if true, print out more information - -int localPatchSize = 1; // (2r+1)*(2r+1)*(2r+1) is the volume of local patch - -// demons registration parameters -// int iterInResolutions[4][3]={{5,3,2},{10,5,5},{15,10,5},{20,15,10}}; -// int itereach = 2; // -// int itereach0 = 0;int itereach1 = 1;int itereach2 = 2;int itereach3 = 3; -// double sigmaDef = 1.5; -// double sigmaDef10 = 1.0;double sigmaDef15 = 1.5;double sigmaDef20 = 2.0; -// double sigmaDef25 = 2.5;double sigmaDef30 = 3.0;double sigmaDef35 = 3.5; -bool doHistMatch = true; -bool doHistMatchTrue = true; bool doHistMatchFalse = false; - -typedef itk::Statistics::MABMISSimulateData DataSimulatorType; -DataSimulatorType::Pointer datasimulator; - -typedef itk::Statistics::MABMISImageOperationFilter ImageOperationFilterType; -ImageOperationFilterType::Pointer imgoperator; -typedef itk::Statistics::MABMISDeformationFieldFilter DeformationFieldOperationFilterType; -DeformationFieldOperationFilterType::Pointer dfoperator; -typedef itk::Statistics::MABMISImageRegistrationFilter ImageRegistrationFilterType; -ImageRegistrationFilterType::Pointer regoperator; -typedef itk::Statistics::MABMISTreeOperation TreeOperationType; -TreeOperationType::Pointer treeoperator; - -typedef itk::Statistics::MABMISBasicOperationFilter BasicOperationFilterType; -BasicOperationFilterType::Pointer basicoperator; - -// - -typedef itk::Vector ShortVectorPixelType; -typedef itk::Image ShortDeformationFieldType; -typedef itk::ImageFileWriter ShortDeformationFieldWriterType; - -void ReadImgInfo(std::vector imgfilenames); - -void HistogramMatching(InternalImageType::Pointer inputImage, InternalImageType::Pointer referenceImage, - InternalImageType::Pointer & outputImage); - -void TreeBasedRegistrationFastOniTree(vnl_vector itree, // the incremental tree - int root, // the tree root - int itree_size, // the tree size - itk::MABMISAtlas* atlasTree, // the atlas tree - itk::MABMISImageData*imageData, // the images to be segmented - std::vector iterations, // registration parameters - double sigma); - -void LabelFusion(std::string curSampleImgName, std::string outSampleSegName, std::vector allWarpedAtlasImgNames, - std::vector allDeformationFieldNames, std::vector allAtlasSegNames, int numOfAtlases); - -void GaussianWeightedLabelFusion(InternalImageType::Pointer curSampleImgPtr, InternalImageType::Pointer outSampleSegPtr, - InternalImageType::Pointer* warpedImgPtrs, InternalImageType::Pointer* warpedSegPtrs, - int numOfAtlases); - -void RegistrationOntoTreeRoot(vnl_vector itree, // the incremental tree - int root, // the tree root - int itree_size, // the tree size - itk::MABMISAtlas* atlasTree, // the atlas tree - itk::MABMISImageData*imageData, // the images to be segmented - std::vector iterations, // registration parameters - double sigma); - -// void PairwiseRegistrationOnTreeViaRoot(int root, int filenumber, int atlas_size, std::vector sub_ids); -void PairwiseRegistrationOnTreeViaRoot(int root, itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree, - std::vector iterations, double sigma); - -// void MultiAtlasBasedSegmentation(int filenumber, int atlas_size, std::vector sub_ids); -int MultiAtlasBasedSegmentation(int root, itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree); - -int Testing(itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree, - std::vector iterations, double sigma) -{ - // Validate the tree size is correct - if( atlasTree->m_TreeSize != atlasTree->m_Tree.size() || - atlasTree->m_TreeSize != atlasTree->m_NumberAllAtlases ) - { - std::cerr << "ERROR: The atlas tree size is incorrect!!" << std::endl; - return -1; - } - ///////////////////////////////////// - // the atlas tree - int tree_size = atlasTree->m_NumberAllAtlases; - - // incremental tree: the atlas tree PLUS images to be segmented - int itree_size = atlasTree->m_NumberAllAtlases + imageData->m_NumberImageData; - vnl_vector itree(itree_size); - for( int i = 0; i < itree_size; ++i ) - { - itree.put(i, 0); - } - // copy current tree into an incremental tree - for( int i = 0; i < tree_size; ++i ) - { - itree.put(i, atlasTree->m_Tree[i]); // combinative tree: all atlases including - // simulated data - } - treeoperator->FindRoot(itree, atlasTree->m_TreeSize, atlasTree->m_TreeRoot); - treeoperator->GetTreeHeight(itree, atlasTree->m_TreeSize, atlasTree->m_TreeHeight); - int root = atlasTree->m_TreeRoot; // root node - // Validate the tree is correct - for( int n = 0; n < atlasTree->m_TreeSize; ++n ) - { - bool valid = true; - if( atlasTree->m_Tree[n] >= atlasTree->m_TreeSize ) - { - valid = false; - } - if( atlasTree->m_Tree[n] < 0 && n != atlasTree->m_TreeRoot ) - { - valid = false; - } - if( !valid ) - { - std::cerr << "ERROR: The atlas tree size is incorrect!!" << std::endl; - return -1; - } - } - - if( isDebug ) - { - treeoperator->SaveTreeWithInfo(itree, itree_size, imageData->m_DataDirectory + "itree.txt"); - } - - // build the incremental tree - std::cout << "---------------------------------------" << std::endl; - std::cout << "Build incremental tree ... " << std::endl; - - int totalNumAtlases = atlasTree->m_NumberAllAtlases; - int totalNumFiles = totalNumAtlases + imageData->m_NumberImageData; - std::vector allfilenames(totalNumFiles); - for( int i = 0; i < totalNumFiles; ++i ) - { - - if( i < totalNumAtlases ) - { - allfilenames[i] = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[i]); - } - else - { - allfilenames[i] = ReplacePathSepForOS(imageData->m_DataDirectory + imageData->m_ImageFileNames[i - totalNumAtlases]); - } - } - - std::cout << "---------------------------------------" << std::endl; - std::cout << "1. Calculating pairwise distance cross training and test images ..." << std::endl; - vnl_matrix cross_dist(imageData->m_NumberImageData, totalNumFiles); - for( int i = 0; i < imageData->m_NumberImageData; ++i ) - { - for( int j = 0; j < totalNumFiles; ++j ) - { - cross_dist.put(i, j, imgoperator->calculateDistanceMSD(allfilenames[totalNumAtlases + i], allfilenames[j]) ); - } - } - - if( isDebug ) - { - basicoperator->SaveMatrix2File(cross_dist, imageData->m_NumberImageData, totalNumFiles, - imageData->m_DataDirectory + "sample2atlas_dist.txt"); - } - - // indices to for atlas files and image files to be segmented. - int* atlas_cur = new int[totalNumFiles]; - for( int i = 0; i < totalNumFiles; ++i ) - { - atlas_cur[i] = 0; - } - for( int i = 0; i < totalNumAtlases; ++i ) - { - atlas_cur[i] = i; - } - int* images_cur = new int[imageData->m_NumberImageData]; - for( int i = 0; i < imageData->m_NumberImageData; ++i ) - { - images_cur[i] = i + totalNumAtlases; - } - - // build the incremental tree - treeoperator->SetAllAtlasNumber(totalNumAtlases); - std::cout << "---------------------------------------" << std::endl; - std::cout << "2. Build the tree... " << std::endl; - itree = treeoperator->BuildIncrementalTree(imageData->m_NumberImageData, // #images to be segmented - totalNumAtlases, // #atlases, including original ones and - // simulated ones - totalNumFiles, // total #image files, = - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // imageData->m_NumberImageData+totalNumAtlases - images_cur, // file indices for images to be segmented - atlas_cur, // file indices for atlas images - cross_dist, // pairwise distance between images - itree); // the output incremental tree - - // delete - delete[] atlas_cur; - delete[] images_cur; - - // check if new tree is actually a tree - const bool iTreeValid = treeoperator->ValidateTree(itree, itree_size); - if( !iTreeValid ) - { - std::cout << "The new built tree is NOT valid!" << std::endl; - - return 105; // EXIT_FAILURE; - } - // save new tree - if( isDebug ) - { - treeoperator->SaveTreeWithInfo(itree, itree_size, imageData->m_DataDirectory + "itree.txt"); - // std::cout << "Tree is built up!" << std::endl; - } - std::cout << "Done! " << std::endl; - - std::cout << "---------------------------------------" << std::endl; - std::cout << "3. Register all images to the root... " << std::endl; - // register all images to the root, if not done yet. - RegistrationOntoTreeRoot(itree, // the incremental tree - root, // the tree root - itree_size, // the tree size - atlasTree, - imageData, - iterations, sigma // registration parameters - ); - std::cout << "Done. " << std::endl; - - std::cout << "---------------------------------------" << std::endl; - std::cout << "4. Pairwise registration via the tree root... " << std::endl; - PairwiseRegistrationOnTreeViaRoot(root, imageData, atlasTree, iterations, sigma); - std::cout << "Done. " << std::endl; - - //////////////////////////////// - // do all multi-atlas based segmentation with weighted label fusion - std::cout << "---------------------------------------" << std::endl; - std::cout << "5. Segment images using label fusion... " << std::endl; - int numIter = MultiAtlasBasedSegmentation(root, imageData, atlasTree); - - std::cout << "---------------------------------------" << std::endl; - std::cout << "Done. " << std::endl; - - // do all multi-atlas based segmentation with weighted label fusion - ///////////////////////////////// - // - // basicoperator->RemoveFile("*_*_cbq_000*"); - // basicoperator->RemoveFile("*_*_deform_000*"); - // - ///////////////////////////////// - - // remove intermediate results - // deformedImageFileName and deformation fields - if( !isDebug ) - { - std::cout << "---------------------------------------" << std::endl; - std::cout << "Removing intermediate registration files..." << std::endl; - for( int n = 0; n < imageData->m_NumberImageData; ++n ) - { - char n_str[10], m_str[10]; - sprintf(n_str, "%03d", n); - for( int m = 0; m < atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; ++m ) - { - sprintf(m_str, "%03d", m); - - std::string imHdr = std::string(m_str) + "_to_" + "testImage" + n_str + "_cbq_000.nii.gz"; - std::string dfMha = std::string(m_str) + "_to_" + "testImage" + n_str + "_deform_000.nii.gz"; - imHdr = imageData->m_DataDirectory + imHdr; - dfMha = imageData->m_DataDirectory + dfMha; - basicoperator->RemoveFile(imHdr.c_str() ); - basicoperator->RemoveFile(dfMha.c_str() ); - - if( m == root ) - { - imHdr = std::string("testImage") + n_str + "_to_" + m_str + "_cbq_reg.nii.gz"; - dfMha = std::string("testImage") + n_str + "_to_" + m_str + "_deform_000.nii.gz"; - imHdr = imageData->m_DataDirectory + imHdr; - dfMha = imageData->m_DataDirectory + dfMha; - basicoperator->RemoveFile(imHdr.c_str() ); - basicoperator->RemoveFile(dfMha.c_str() ); - } - } - for( int m = 0; m < imageData->m_NumberImageData; ++m ) - { - if( m == n ) - { - continue; - } - sprintf(m_str, "%03d", m); - std::string imHdr = std::string("testImage") + n_str + "_to_testImage" + m_str + "_cbq_000.nii.gz"; - std::string dfMha = std::string("testImage") + n_str + "_to_testImage" + m_str + "_deform_000.nii.gz"; - imHdr = imageData->m_DataDirectory + imHdr; - dfMha = imageData->m_DataDirectory + dfMha; - basicoperator->RemoveFile(imHdr.c_str() ); - basicoperator->RemoveFile(dfMha.c_str() ); - } - } - } - // remove segmentation results in iterations. - for( int n = 0; n < imageData->m_NumberImageData; ++n ) - { - const std::string imageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[n]; - const std::string baseFileName = GetRootName(imageFileName); - char i_str[10]; - bool copied = false; - for( int i = numIter - 1; i >= 0; i-- ) - { - sprintf(i_str, "%03d", i); - const std::string segHdr = baseFileName + "_seg_" + i_str + ".nii.gz"; - if( itksys::SystemTools::FileExists(segHdr.c_str(), true) ) - { - if( !copied ) - { - std::string segHdr_save = baseFileName + "_seg.nii.gz"; - - const size_t dir_sep = segHdr_save.find_last_of(FILESEP); - if( dir_sep != std::string::npos ) - { - segHdr_save = imageData->m_OutputDirectory + segHdr_save.substr(dir_sep + 1, std::string::npos); - } - else - { - segHdr_save = imageData->m_OutputDirectory + segHdr_save; - } - - itksys::SystemTools::CopyFileAlways(segHdr.c_str(), segHdr_save.c_str() ); - copied = true; - } - basicoperator->RemoveFile(segHdr.c_str() ); - } - } - } - - std::cout << "Done. " << std::endl; - return EXIT_SUCCESS; -} - -template -int DoIt( itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree, - std::vector iterations, double sigma) -{ - return Testing(imageData, atlasTree, iterations, sigma); -} +#include "Testing.h" int main( int argc, char *argv[] ) { @@ -590,6 +70,7 @@ int main( int argc, char *argv[] ) } inputImageData->m_OutputDirectory = OutputFolder; } + itksys::SystemTools::MakeDirectory(inputImageData->m_OutputDirectory); // load the tree-structured atlases that is generated in the training step itk::MABMISAtlasXMLFileReader::Pointer treeAtlasXMLReader = itk::MABMISAtlasXMLFileReader::New(); @@ -623,1065 +104,9 @@ int main( int argc, char *argv[] ) atlasTree->m_AtlasDirectory = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + FILESEP ); } - // Will look into getting rid of these global variables later ---Xiaofeng - datasimulator = DataSimulatorType::New(); - imgoperator = ImageOperationFilterType::New(); - dfoperator = DeformationFieldOperationFilterType::New(); - regoperator = ImageRegistrationFilterType::New(); - treeoperator = TreeOperationType::New(); - basicoperator = BasicOperationFilterType::New(); - - int retVal = DoIt( inputImageData, atlasTree, iterations, SmoothingKernelSize); - + int retVal = Testing( inputImageData, atlasTree, iterations, SmoothingKernelSize); delete inputImageData; delete atlasTree; return retVal; } -void LabelFusion(std::string curSampleImgName, std::string outSampleSegName, std::vector allWarpedAtlasImgNames, - std::vector allDeformationFieldNames, std::vector allAtlasSegNames, int numOfAtlases) -{ - // create pointers for each images including: cur_image, our_seg, - // warpedImg(after histogram matching), warpedSeg (after warping) - - // cur sample image pointer - InternalImageType::Pointer curSampleImgPtr = 0; - - imgoperator->ReadImage(curSampleImgName, curSampleImgPtr); - - // output sample label image pointer - InternalImageType::Pointer outSampleSegPtr = 0; - imgoperator->ReadImage(curSampleImgName, outSampleSegPtr); - - // atlas related images pointers - InternalImageType::Pointer* warpedImgPtrs = new InternalImageType::Pointer[numOfAtlases]; - InternalImageType::Pointer* warpedSegPtrs = new InternalImageType::Pointer[numOfAtlases]; - InternalImageType::IndexType index; - index[0] = 19; index[1] = 19; index[2] = 8; - for( int i = 0; i < numOfAtlases; ++i ) - { - // histogram match each warped atlas to the current sample - warpedImgPtrs[i] = 0; - InternalImageType::Pointer curWarpedAtlasPtr = 0; - imgoperator->ReadImage(allWarpedAtlasImgNames[i], curWarpedAtlasPtr); - HistogramMatching(curWarpedAtlasPtr, curSampleImgPtr, warpedImgPtrs[i]); - - // load each deformation field - DeformationFieldType::Pointer deformFieldPtr = 0; - dfoperator->ReadDeformationField(allDeformationFieldNames[i], deformFieldPtr); - - // warp each atlas label image to the current sample space - InternalImageType::Pointer curAtlasSegPtr = 0; - imgoperator->ReadImage(allAtlasSegNames[i], curAtlasSegPtr); - dfoperator->ApplyDeformationField(curAtlasSegPtr, deformFieldPtr, warpedSegPtrs[i], false); - } - - // weighted label fusion - GaussianWeightedLabelFusion(curSampleImgPtr, outSampleSegPtr, warpedImgPtrs, warpedSegPtrs, numOfAtlases); - - // output segmentation image - imgoperator->WriteImage(outSampleSegName, outSampleSegPtr); - - delete[] warpedImgPtrs; - delete[] warpedSegPtrs; - - return; -} - -void GaussianWeightedLabelFusion(InternalImageType::Pointer curSampleImgPtr, InternalImageType::Pointer outSampleSegPtr, - InternalImageType::Pointer* warpedImgPtrs, InternalImageType::Pointer* warpedSegPtrs, - int numOfAtlases) -{ - // introduce the iterators of each image - // InternalImageIteratorType sampleImgIt(curSampleImgPtr, - // curSampleImgPtr->GetLargestPossibleRegion());//GetRequestedRegion()); - InternalImageIteratorType sampleSegIt(outSampleSegPtr, outSampleSegPtr->GetLargestPossibleRegion() ); // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // GetRequestedRegion()); - - InternalImageType::SizeType sampleImSize = outSampleSegPtr->GetBufferedRegion().GetSize(); - - // InternalImageIteratorType* warpedImgIts = new InternalImageIteratorType[numOfAtlases]; - InternalImageIteratorType* warpedSegIts = new InternalImageIteratorType[numOfAtlases]; - - typedef itk::ConstNeighborhoodIterator InternalImageNeighborhoodIteratorType; - InternalImageNeighborhoodIteratorType* warpedImgNeighborhoodIts = - new InternalImageNeighborhoodIteratorType[numOfAtlases]; - InternalImageType::SizeType radius; - radius[0] = localPatchSize; radius[1] = localPatchSize; radius[2] = localPatchSize; - - InternalImageNeighborhoodIteratorType sampleImgNeighborhoodIt(radius, curSampleImgPtr, - curSampleImgPtr->GetRequestedRegion() ); - - InternalImageNeighborhoodIteratorType sampleImgNeighborIt(radius, curSampleImgPtr, - curSampleImgPtr->GetLargestPossibleRegion() ); - for( int i = 0; i < numOfAtlases; ++i ) - { - // InternalImageIteratorType it1( warpedImgPtrs[i], - // warpedImgPtrs[i]->GetLargestPossibleRegion());//GetRequestedRegion() ); - // warpedImgIts[i] = it1; - - InternalImageIteratorType it2( warpedSegPtrs[i], warpedSegPtrs[i]->GetLargestPossibleRegion() ); // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // GetRequestedRegion() - // ); - warpedSegIts[i] = it2; - - InternalImageNeighborhoodIteratorType neighborIt(radius, warpedImgPtrs[i], - warpedImgPtrs[i]->GetLargestPossibleRegion() ); - warpedImgNeighborhoodIts[i] = neighborIt; - } - int maxLabel = 0; // the maximum labels - // find out the maximum label from the label images - for( int i = 0; i < numOfAtlases; ++i ) - { - warpedSegIts[i].GoToBegin(); - } - while( !warpedSegIts[0].IsAtEnd() ) - { - for( int i = 0; i < numOfAtlases; ++i ) - { - if( warpedSegIts[i].Get() > maxLabel ) - { - maxLabel = warpedSegIts[i].Get(); - } - ++warpedSegIts[i]; - } - } - - InternalImageIteratorType::IndexType idx; - InternalImageIteratorType::IndexType idx1; - - int temp1, temp2; - double kernel_sigma = 1; // std of Gaussian approximation, set to the median of all distances - int* index = new int[numOfAtlases]; - double* sort1 = new double[numOfAtlases]; - int* label_index = new int[maxLabel + 1]; - - int* label_pool = new int[numOfAtlases]; - double* mse = new double[numOfAtlases]; - double* weight_sum = new double[maxLabel + 1]; - // for each voxel do label fusion,(vx, vy, vz) is the current location - double kernel_sigma_square = kernel_sigma * kernel_sigma; - for( int i = 0; i < numOfAtlases; ++i ) - { - warpedImgNeighborhoodIts[i].GoToBegin(); - warpedSegIts[i].GoToBegin(); - } - for( sampleImgNeighborIt.GoToBegin(), sampleSegIt.GoToBegin(); - !sampleImgNeighborIt.IsAtEnd(); - ++sampleImgNeighborIt, ++sampleSegIt ) - { -// idx = sampleImgNeighborIt.GetIndex(); -// if (idx[0] == 30 && idx[1]==20 && idx[2] ==8) -// idx[0] = idx[0]; - double msec = 0.0; - for( int i = 0; i < numOfAtlases; ++i ) - { - label_pool[i] = warpedSegIts[i].Get(); - for( int n = 0; n < sampleImgNeighborIt.Size(); ++n ) - { - temp1 = sampleImgNeighborIt.GetPixel(n); - temp2 = warpedImgNeighborhoodIts[i].GetPixel(n); - // add together differences, squared sum - msec += (temp1 - temp2) * (temp1 - temp2); - } - mse[i] = msec; - } - // sort the mean square difference - for( int i = 0; i < numOfAtlases; ++i ) - { - index[i] = i; - sort1[i] = mse[i]; - } - basicoperator->bubbleSort(sort1, index, numOfAtlases); - kernel_sigma_square = sort1[(int)( (numOfAtlases - 1) / 2)] + 0.0001; - // weight each label - for( int i = 0; i < maxLabel + 1; ++i ) - { - weight_sum[i] = 0.0; - } - for( int i = 0; i < numOfAtlases; ++i ) - { - // add-up all weights - if ((mse[i]) / (2.0 * kernel_sigma_square) < 50) - weight_sum[label_pool[i]] += exp(0.0 - (mse[i]) / (2.0 * kernel_sigma_square) ); - - } - - // weighted label fusion - for( int i = 0; i < maxLabel + 1; ++i ) - { - label_index[i] = i; - } - - // determine the label with the maximum weight - // basicoperator->bubbleSort(weight_sum, label_index, maxLabel+1); - // int label_final; - // label_final = label_index[maxLabel]; - int label_final = 0; - double weight_final = 0.0; - for( int i = 0; i < maxLabel + 1; ++i ) - { - if( weight_sum[i] > weight_final ) - { - label_final = i; - weight_final = weight_sum[i]; - } - } - - // assign label to the segmentations - sampleSegIt.Set(label_final); - for( int i = 0; i < numOfAtlases; ++i ) - { - ++warpedImgNeighborhoodIts[i]; - ++warpedSegIts[i]; - } - } - /* - for (int vz = 0; vz < sampleImSize[2]; vz++) - { - //if (isDebug) - // std::cout << vz << ", "; - for (int vy = 0; vy < sampleImSize[1]; vy++) - { - for (int vx = 0; vx < sampleImSize[0]; vx++) - { - // process each voxel - idx.SetElement(0,vx); - idx.SetElement(1,vy); - idx.SetElement(2,vz); - - // get labels from all atlases on this current location - - for (int i = 0; i < numOfAtlases; ++i) - { - - warpedSegIts[i].SetIndex(idx); - label_pool[i] = warpedSegIts[i].Get(); - } - - // get the local patch differences - - for (int i = 0; i < numOfAtlases; ++i) mse[i] = 0.0; - for (int i = 0; i < numOfAtlases; ++i) - { - double msec = 0.0; - // compare with all other images at the local patch of current voxel - // (nx,ny,nz) is the incremental position on (vx,vy,vz) - for(int nz = 0-localPatchSize; nz <= localPatchSize; nz++) - { - for(int ny = 0-localPatchSize; ny <= localPatchSize; ny++) - { - for(int nx = 0-localPatchSize; nx <= localPatchSize; nx++) - { - if ((vx+nx>=0)&&(vx+nx=0)&&(vy+ny=0)&&(vz+nzbubbleSort(sort1, index, numOfAtlases); - kernel_sigma = sort1[(int)((numOfAtlases-1)/2)]+0.0001; - - // weight each label - - for (int i = 0; i < maxLabel+1; ++i) weight_sum[i] = 0.0; - for (int i = 0; i < numOfAtlases; ++i) - { - // add-up all weights - weight_sum[label_pool[i]]+= exp(0-(mse[i]*mse[i])/(2*kernel_sigma*kernel_sigma)); - } - - // weighted label fusion - for (int i = 0; i < maxLabel+1; ++i) label_index[i] = i; - - basicoperator->bubbleSort(weight_sum, label_index, maxLabel+1); - int label_final; - label_final = label_index[maxLabel]; - - // assign label to the segmentations - sampleSegIt.SetIndex(idx); - sampleSegIt.Set(label_final); - - } - } - }// end of for vz - -*/ - delete[] label_pool; - delete[] mse; - delete[] weight_sum; - // - // delete[] warpedImgIts; - delete[] warpedSegIts; - delete[] warpedImgNeighborhoodIts; - delete[] sort1; - delete[] index; - delete[] label_index; - - return; -} - -void TreeBasedRegistrationFastOniTree(vnl_vector itree, // the incremental tree - int root, // the tree root - int itree_size, // the tree size - itk::MABMISAtlas* atlasTree, // the atlas tree - itk::MABMISImageData*imageData, // the images to be segmented - std::vector iterations, - double sigma - ) -// (vnl_vector itree_t, int itree_size_t, int atlas_size_t, int simulate_size_t, int sample_size_t, std::vector -// originalIntensityImageFileNames, std::vector deformedImgFileNames, std::vector deformationFieldFileNames) -{ - // int root = -1; // - // treeoperator->FindRoot(itree_t, itree_size_t, root); - - // calculate the depth of each node - int* node_depth = new int[itree_size]; - - for( int i = 0; i < itree_size; ++i ) - { - int curnode = i; - node_depth[i] = 0; - bool isRoot = false; - - while( !isRoot ) - { - if( itree[curnode] >= 0 ) - { - node_depth[i]++; - curnode = itree[curnode]; - } - else - { - isRoot = true; - } - } - } - - // sort all nodes by its depth - int* index = new int[itree_size]; - { - // - double* arr = new double[itree_size]; - for( int i = 0; i < itree_size; ++i ) - { - arr[i] = node_depth[i]; - index[i] = i; - } - basicoperator->bubbleSort(arr, index, itree_size); // index[0] = root - delete[] arr; - } - - int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; - int atlas_total_size = atlasTree->m_NumberAllAtlases; - - // start to register each image to the root node step by step - for( int ii = 1; ii < itree_size; ++ii ) // starting from 1, since index[0] showing the root - { - // if (isDebug) - // std::cout << ii << ", "; - int curnode = index[ii]; - int parentnode = itree[curnode]; - - std::cout << "Between nodes: [" << curnode << ", " << parentnode << "]" << std::endl; - - if( (curnode >= atlas_image_size) && (curnode < atlas_total_size) ) - { - // this is a simulated image, then pass - continue; - } - - // first check whether registrations from atlas images to root exist. If exist, no registration is needed - bool isRegistered = false; - char i_str[10], root_str[10], p_str[10]; - sprintf(i_str, "%03d", curnode); - sprintf(root_str, "%03d", root); - sprintf(p_str, "%03d", parentnode); - - std::string dfFileName, df_ImageFileName; - if( curnode == root ) - { - isRegistered = true; - } - else if( curnode < atlas_image_size ) - { - dfFileName = std::string(i_str) + "_to_" + root_str + "_deform_000.nii.gz"; - df_ImageFileName = std::string(i_str) + "_to_" + root_str + "_cbq_reg.nii.gz"; - dfFileName = atlasTree->m_AtlasDirectory + dfFileName; - df_ImageFileName = atlasTree->m_AtlasDirectory + df_ImageFileName; - - if( (itksys::SystemTools::FileExists(dfFileName.c_str(), true) && - itksys::SystemTools::FileExists(df_ImageFileName.c_str(), true) ) - ) - { - isRegistered = true; - } - } - else if( curnode >= atlas_total_size ) - { - int n = curnode - atlas_total_size; - char n_str[10]; sprintf(n_str, "%03d", n); - dfFileName = std::string("testImage") + n_str + "_to_" + root_str + "_deform_000.nii.gz"; - df_ImageFileName = std::string("testImage") + n_str + "_to_" + root_str + "_cbq_reg.nii.gz"; - - dfFileName = imageData->m_DataDirectory + dfFileName; - df_ImageFileName = imageData->m_DataDirectory + df_ImageFileName; - - isRegistered = false; - } - // parent deformation fields - std::string df_ParentFileName; - if( parentnode < atlas_image_size ) - { - df_ParentFileName = std::string(p_str) + "_to_" + root_str + "_deform_000.nii.gz"; - df_ParentFileName = atlasTree->m_AtlasDirectory + df_ParentFileName; - } - else if( parentnode >= atlas_total_size ) - { - int n = parentnode - atlas_total_size; - char n_str[10]; sprintf(n_str, "%03d", n); - df_ParentFileName = std::string("testImage") + n_str + "_to_" + root_str + "_deform_000.nii.gz"; - df_ParentFileName = imageData->m_DataDirectory + df_ParentFileName; - } - - // do registration - if( isRegistered ) - { - continue; // goto next one - } - - std::string rootImageFile = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[root]); - std::string testImageFile; - if( curnode < atlas_image_size ) - { - testImageFile = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[curnode]); // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // atlas - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // - // image - } - else - { - testImageFile = imageData->m_DataDirectory + imageData->m_ImageFileNames[curnode - atlas_total_size]; // test - // image - } - if( parentnode == root ) // direct registration without initial deformation - { - regoperator->DiffeoDemonsRegistrationWithParameters( - rootImageFile, - testImageFile, - df_ImageFileName, - dfFileName, - sigma, doHistMatch, iterations); - } - else // registration with initial deformation - { - if( (parentnode < atlas_image_size) || (parentnode >= atlas_total_size) ) - { - // parent is a real image - regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( - rootImageFile, - testImageFile, - df_ParentFileName, - df_ImageFileName, - dfFileName, - sigma, doHistMatch, iterations); - } - else - { - // parent is a simulated image - - std::string index_string; basicoperator->myitoa( parentnode - atlas_image_size, index_string, 3 ); - std::string sim_df_FileName = "simulated_deform_" + index_string + ".nii.gz"; - sim_df_FileName = atlasTree->m_AtlasDirectory + sim_df_FileName; - - regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( - rootImageFile, - testImageFile, - sim_df_FileName, - df_ImageFileName, - dfFileName, - sigma, doHistMatch, iterations); - } - } - } - std::cout << "Done." << std::endl; - - delete[] node_depth; - delete[] index; -} - -// other filter -// histogram matching -void HistogramMatching(InternalImageType::Pointer inputImage, - InternalImageType::Pointer referenceImage, - InternalImageType::Pointer & outputImage) -{ - InternalHistMatchFilterType::Pointer matcher = InternalHistMatchFilterType::New(); - - matcher->SetInput( inputImage ); - matcher->SetReferenceImage( referenceImage ); - matcher->SetNumberOfHistogramLevels( 1024 ); // 1024 - matcher->SetNumberOfMatchPoints( 7 ); - matcher->ThresholdAtMeanIntensityOn(); - try - { - matcher->Update(); - } - catch( itk::ExceptionObject & err ) - { - std::cout << err << std::endl; - return; - } - outputImage = matcher->GetOutput(); - return; -} - -// void RegistrationOntoTreeRoot(vnl_vector itree, int root, int filenumber, int itree_size, int atlas_size, int -// simulate_size, int sample_size, std::vector sub_ids, std::vector segfilenames) -void RegistrationOntoTreeRoot(vnl_vector itree, // the incremental tree - int root, // the tree root - int itree_size, // the tree size - itk::MABMISAtlas* atlasTree, // the atlas tree - itk::MABMISImageData*imageData, // the images to be segmented - std::vector iterations, // registration parameters - double sigma - ) -{ - // do registration on the tree - std::cout << "Start registration ... " << std::endl; - - // to get the numbers original defined in the data; - const int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; - const int test_image_size = imageData->m_NumberImageData; - -// do tree-based registration - TreeBasedRegistrationFastOniTree(itree, // the incremental tree - root, // the tree root - itree_size, // the tree size - atlasTree, // the atlas tree - imageData, // the images to be segmented - iterations, sigma // registration parameters - ); - - // end of do registration to the root - - // - // basicoperator->RemoveFile("simulated_*"); - - std::cout << "Reverse deformation field ... " << std::endl;; - // reverse each deformation field (from the root) - std::string rootImageTag; - char root_str[10], i_str[10]; - sprintf(root_str, "%03d", root); - rootImageTag = root_str; - std::string fixedImageTag; - - for( int i = 0; i < atlas_image_size + test_image_size; i++ ) - { - //if( isDebug ) - { - std::cout << i << ", "; - } - if( i == root ) - { - continue; - } - - // prepare file names - - // std::string originalImgImageFileName = sub_ids[root] + "_cbq_000.nii.gz"; - // std::string originalSegImageFileName = sub_ids[root] + "_seg_000.nii.gz"; - std::string rootImageFileName = ReplacePathSepForOS( - atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[root] ); - std::string rootSegmentFileName = ReplacePathSepForOS( - atlasTree->m_AtlasDirectory + atlasTree->m_AtlasSegmentationFilenames[root] ); - - std::string fixedImageFileName; - std::string deformedImageFileName; - std::string deformedImageFileNameImg; - std::string deformedSegmentFileName; - - std::string deformationFileName; - std::string invDeformationFileName; - - if( i < atlas_image_size ) - { - fixedImageFileName = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[i]); - sprintf(i_str, "%03d", i); - fixedImageTag = i_str; - - deformedImageFileName = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + +"_cbq_000.nii.gz"; - deformedImageFileNameImg = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + +"_cbq_000.nii.gz"; - - deformedSegmentFileName = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + +"_seg_000.nii.gz"; - - deformationFileName = atlasTree->m_AtlasDirectory + fixedImageTag + "_to_" + rootImageTag + "_deform_000.nii.gz"; - invDeformationFileName = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + "_deform_000.nii.gz"; - } - else - { - fixedImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[i - atlas_image_size]; - sprintf(i_str, "%03d", i - atlas_image_size); - fixedImageTag = std::string("testImage") + i_str; - - deformedImageFileName = imageData->m_DataDirectory + rootImageTag + "_to_" + fixedImageTag + +"_cbq_000.nii.gz"; - deformedImageFileNameImg = imageData->m_DataDirectory + rootImageTag + "_to_" + fixedImageTag + +"_cbq_000.nii.gz"; - - deformedSegmentFileName = imageData->m_DataDirectory + rootImageTag + "_to_" + fixedImageTag + +"_seg_000.nii.gz"; - - deformationFileName = imageData->m_DataDirectory + fixedImageTag + "_to_" + rootImageTag + "_deform_000.nii.gz"; - invDeformationFileName = imageData->m_DataDirectory + rootImageTag + "_to_" + fixedImageTag + "_deform_000.nii.gz"; - } - - // if exist?? - bool isReversedExist = true; - if( !itksys::SystemTools::FileExists(deformedImageFileName.c_str(), true) ) - { - isReversedExist = false; - } - if( !itksys::SystemTools::FileExists(deformedImageFileNameImg.c_str(), true) ) - { - isReversedExist = false; - } - if( !itksys::SystemTools::FileExists(invDeformationFileName.c_str(), true) ) - { - isReversedExist = false; - } - if( isReversedExist ) - { - continue; - } - - // reverse deformation field - DeformationFieldType::Pointer deformationField = 0; - dfoperator->ReadDeformationField(deformationFileName, deformationField); - - DeformationFieldType::Pointer inversedDeformationField = DeformationFieldType::New(); - inversedDeformationField->SetRegions(deformationField->GetLargestPossibleRegion() ); - inversedDeformationField->SetSpacing(deformationField->GetSpacing() ); - inversedDeformationField->SetDirection(deformationField->GetDirection() ); - inversedDeformationField->SetOrigin(deformationField->GetOrigin() ); - inversedDeformationField->Allocate(); - - dfoperator->InverseDeformationField3D(deformationField, inversedDeformationField); - // output reversed deformation field - - dfoperator->WriteDeformationField(invDeformationFileName, inversedDeformationField); - // //update - regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( - fixedImageFileName, rootImageFileName, - invDeformationFileName, - deformedImageFileName, invDeformationFileName, - sigma, doHistMatch, iterations); - - // CompressDeformationField2Short(inversedDeformationFieldFileName); - // apply deformation field on seg file - if( isEvaluate ) - { - dfoperator->ApplyDeformationFieldAndWriteWithFileNames( - rootSegmentFileName, - invDeformationFileName, - deformedSegmentFileName, false); - } - - std::cout << "Done!" << std::endl; - } // end of reverse each deformation field (from the root) - // std::cout << "done." << std::endl; - std::cout << std::endl; -} - -void PairwiseRegistrationOnTreeViaRoot(int root, - itk::MABMISImageData* imageData, - itk::MABMISAtlas* atlasTree, - std::vector iterations, - double sigma) -{ - ///////////////////////////////// - // with the help of root node, obtain the pairwise registration among all images - std::cout << "Generate all pairwise registration ..." << std::endl;; - - int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; - int test_image_size = imageData->m_NumberImageData; - - // do for all real images, including both atlases and test images - for( int all_index = 0; all_index < atlas_image_size + test_image_size; all_index++ ) - { - if( isDebug ) - { - std::cout << all_index << ": "; - } - if( all_index == root ) - { - continue; - } - - char a_str[10], s_str[10], root_str[10]; - sprintf(root_str, "%03d", root); - sprintf(a_str, "%03d", all_index); - - std::string movingImageFileName; - std::string movingSegmentFileName; - // get the moving filenames - std::string movingImageTag; - if( all_index < atlas_image_size ) - { - movingImageFileName = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[all_index]); - movingSegmentFileName = atlasTree->m_AtlasDirectory + atlasTree->m_AtlasSegmentationFilenames[all_index]; - movingImageTag = std::string(a_str); - } - else - { - movingImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[all_index - atlas_image_size]; - movingSegmentFileName = movingImageFileName; - const std::string movingSegmentBaseFileName = GetRootName(movingSegmentFileName); - - movingSegmentFileName = movingSegmentBaseFileName + "_seg_000.nii.gz"; - sprintf(a_str, "%03d", all_index - atlas_image_size); - movingImageTag = std::string("testImage") + a_str; - } - // for (int sample_index = atlas_size; sample_index < filenumber; sample_index++) - for( int sample_index = 0; sample_index < imageData->m_NumberImageData; sample_index++ ) - { - if( isDebug ) - { - std::cout << "[" << all_index << "," << sample_index << "], "; - } - - if( sample_index == all_index - atlas_image_size ) - { - continue; - } - - std::string fixedImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[sample_index]; - - sprintf(s_str, "%03d", sample_index); - std::string fixedImageTag = std::string("testImage") + s_str; - - std::string deformedImageFileName = movingImageTag + "_to_" + fixedImageTag + "_cbq_000.nii.gz"; - std::string deformedImageFileNameImg = movingImageTag + "_to_" + fixedImageTag + "_cbq_000.nii.gz"; - std::string deformedSegmentFileName = movingImageTag + "_to_" + "testImage" + s_str + "_seg_000.nii.gz"; - std::string deformedSegmentFileNameImg = movingImageTag + "_to_" + "testImage" + s_str + "_seg_000.nii.gz"; - std::string dfFileName = movingImageTag + "_to_" + fixedImageTag + "_deform_000.nii.gz"; - - deformedImageFileName = imageData->m_DataDirectory + deformedImageFileName; - deformedImageFileNameImg = imageData->m_DataDirectory + deformedImageFileNameImg; - dfFileName = imageData->m_DataDirectory + dfFileName; - - deformedSegmentFileName = imageData->m_DataDirectory + deformedSegmentFileName; - deformedSegmentFileNameImg = imageData->m_DataDirectory + deformedSegmentFileNameImg; - - // std::string fixedImgImageFileName = sub_ids[sample_index] + "_cbq_000.nii.gz"; - // std::string movingImgImageFileName = sub_ids[all_index] + "_cbq_000.nii.gz"; - // std::string movingSegImageFileName = sub_ids[all_index] + "_seg_000.nii.gz"; - - // std::string deformedImgImageFileName = sub_ids[sample_index] + "_" + sub_ids[all_index] + "_cbq_000.nii.gz"; - // std::string deformedSegImageFileName = sub_ids[sample_index] + "_" + sub_ids[all_index] + "_seg_000.nii.gz"; - // std::string deformedImgImageFileNameImg = sub_ids[sample_index] + "_" + sub_ids[all_index] + "_cbq_000.nii.gz"; - // std::string deformedSegImageFileNameImg = sub_ids[sample_index] + "_" + sub_ids[all_index] + "_seg_000.nii.gz"; - // std::string deformationFieldFileName = sub_ids[sample_index] + "_" + sub_ids[all_index] + "_deform_000.nii.gz"; - - // if exist - bool isComposeExist = true; - if( !itksys::SystemTools::FileExists(deformedImageFileName.c_str(), true) ) - { - isComposeExist = false; - } - if( !itksys::SystemTools::FileExists(deformedImageFileNameImg.c_str(), true) ) - { - isComposeExist = false; - } - - if( !itksys::SystemTools::FileExists(dfFileName.c_str(), true) ) - { - isComposeExist = false; - } - if( isComposeExist ) - { - continue; - } - - // compose - // std::string inputDeformationFieldFileName = sub_ids[root] + "_" + sub_ids[all_index] + "_deform_000.nii.gz"; - // std::string secondDeformationFieldFileName = sub_ids[sample_index] + "_" + sub_ids[root] + "_deform_000.nii.gz"; - std::string inputDeformationFieldFileName; - if( all_index < atlas_image_size ) - { - inputDeformationFieldFileName = atlasTree->m_AtlasDirectory + movingImageTag + "_to_" + root_str - + "_deform_000.nii.gz"; - } - else - { - inputDeformationFieldFileName = imageData->m_DataDirectory + movingImageTag + "_to_" + root_str - + "_deform_000.nii.gz"; - } - std::string secondDeformationFieldFileName = imageData->m_DataDirectory + root_str + "_to_" + fixedImageTag - + "_deform_000.nii.gz"; - // output composed deformation field - dfoperator->ComposeDeformationFieldsAndSave(inputDeformationFieldFileName, - secondDeformationFieldFileName, - dfFileName); - - // update: refine the deformation - regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( - fixedImageFileName, - movingImageFileName, - dfFileName, - deformedImageFileName, - dfFileName, - sigma, doHistMatch, iterations); - - if( isEvaluate ) - { - dfoperator->ApplyDeformationFieldAndWriteWithFileNames( - movingSegmentFileName, - dfFileName, - deformedSegmentFileName, false); - } - } // end of sample_index - } // end of all_index - - std::cout << std::endl; - std::cout << "Done!" << std::endl; -} - -// void MultiAtlasBasedSegmentation(int filenumber, int atlas_size, std::vector sub_ids) -int MultiAtlasBasedSegmentation(int root, - itk::MABMISImageData* imageData, - itk::MABMISAtlas* atlasTree) -{ - std::cout << "Start to process segmentation..." << std::endl; - int iter = 1; - - bool isConverged = false; // the iterative update can stop - - int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; - - int test_image_size = imageData->m_NumberImageData; - - while( !isConverged ) - { - // if (isDebug) - std::cout << "iteration " << iter << ": "; - // iter == 1, only using atlases; otherwise using all images - // number of other image to be used - int numOtherAtlases = 0; - if( iter == 1 ) - { - numOtherAtlases = atlas_image_size; - } - else - { - numOtherAtlases = atlas_image_size + test_image_size - 1; - } - // do label fusion for each new sample - for( int sample_index = 0; sample_index < test_image_size; sample_index++ ) - { - // if (isDebug) - // std::cout << sample_index << ":" << sub_ids[sample_index] << ", "; - std::cout << sample_index << ", "; - - char s_str[10]; - sprintf(s_str, "%03d", sample_index); - std::string testImageTag = std::string("testImage") + s_str; - - // the sample to be processed - - std::string testImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[sample_index]; - - // std::string curSampleImgName = sub_ids[sample_index] +"_cbq_000.nii.gz"; - // the output label image of current sample - std::string index_string; - basicoperator->myitoa( iter, index_string, 3 ); - std::string testSementFileName = testImageFileName; - const std::string testSegmentBaseFileName = GetRootName(testSementFileName); - testSementFileName = testSegmentBaseFileName + "_seg_" + index_string + ".nii.gz"; - - // std::string outSampleSegName = sub_ids[sample_index] + "_seg_" + index_string + ".nii.gz"; - - // prepare all current atlases names - std::vector allWarpedAtlasImgNames; - std::vector allDeformationFieldNames; - std::vector allAtlasSegNames; - - int atlas_index = 0; - for( int i = 0; i < atlas_image_size + test_image_size; ++i ) - { - if( (iter == 1) && (i >= atlas_image_size) ) - { - continue; - } - if( i - atlas_image_size == sample_index ) - { - continue; - } - - std::string movingImageTag; - char i_str[10]; - - if( i < atlas_image_size ) - { - sprintf(i_str, "%03d", i); - movingImageTag = std::string(i_str); - } - else - { - sprintf(i_str, "%03d", i - atlas_image_size); - movingImageTag = std::string("testImage") + i_str; - } - - // assign names - // std::string allWarpedAtlasImgName = sub_ids[sample_index] + "_" + sub_ids[i] + "_cbq_000.nii.gz"; - // std::string allDeformationFieldName = sub_ids[sample_index] + "_" + sub_ids[i] + "_deform_000.nii.gz"; - std::string warpedImageFileName = movingImageTag + "_to_" + testImageTag + "_cbq_000.nii.gz"; - std::string dfFileName = movingImageTag + "_to_" + testImageTag + "_deform_000.nii.gz"; - - warpedImageFileName = imageData->m_DataDirectory + warpedImageFileName; - dfFileName = imageData->m_DataDirectory + dfFileName; - - std::string atlasSegName; - if( i < atlas_image_size ) - { - atlasSegName = atlasTree->m_AtlasDirectory + atlasTree->m_AtlasSegmentationFilenames[i]; - } - else - { - atlasSegName = imageData->m_DataDirectory + imageData->m_ImageFileNames[i - atlas_image_size]; - const std::string atlasSegBaseName = GetRootName( atlasSegName ); - - basicoperator->myitoa( iter - 1, index_string, 3 ); - atlasSegName = atlasSegBaseName + "_seg_" + index_string + ".nii.gz"; - } - - allWarpedAtlasImgNames.push_back(warpedImageFileName); - allDeformationFieldNames.push_back(dfFileName); - allAtlasSegNames.push_back(atlasSegName); - - // go to next atlas - atlas_index++; - } - // do weighted label fusion - - LabelFusion(testImageFileName, // the input image to be segmented - testSementFileName, // the output segmented image - allWarpedAtlasImgNames, - allDeformationFieldNames, - allAtlasSegNames, - numOtherAtlases); - } - std::cout << std::endl; - - // to decide if converged - if( iter >= 3 ) - { - isConverged = true; // - } - // calculate overlap rate - if( isEvaluate ) - { - // pairwise segmentation accuracy - - // pairwise segmentation consistency - } - iter++; - } - - return iter; -} - diff --git a/IGR3D_MABMIS_Training.cxx b/IGR3D_MABMIS_Training.cxx index 0c39772..2ec1ae7 100644 --- a/IGR3D_MABMIS_Training.cxx +++ b/IGR3D_MABMIS_Training.cxx @@ -8,499 +8,9 @@ PURPOSE. See the above copyright notices for more information. =========================================================================*/ -// for math -#include -#include -#include -#include -#include -#include -#include -#include -#include -// basic itk -#include "itkVector.h" - -// registration -#include "itkImageRegistrationMethod.h" -#include "itkSymmetricForcesDemonsRegistrationFilter.h" - -// interpolators -#include "itkLinearInterpolateImageFunction.h" -#include "itkNearestNeighborInterpolateImageFunction.h" - -// reader / writer -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" - -// filter -#include "itkResampleImageFilter.h" - -#include "itkCastImageFilter.h" -#include "itkWarpImageFilter.h" - -#include "itkHistogramMatchingImageFilter.h" -#include "itkAddImageFilter.h" -#include "itkDivideByConstantImageFilter.h" -#include "itkMultiplyByConstantImageFilter.h" -#include "itkWarpVectorImageFilter.h" -#include "itkInverseDeformationFieldImageFilter.h" - -// for affine transformation -#include "itkTransform.h" -#include "itkAffineTransform.h" -#include "itkImageRegistrationMethod.h" - -// for Diffeomorphic Demons -#include -#include -#include - -// including itksys::SystemTools::MakeDirectory(char*) -#include -#include - -// To include all related header files #include "IGR3D_MABMIS_TrainingCLP.h" -#include "itkMABMISImageOperationFilter.h" -#include "itkMABMISDeformationFieldFilter.h" -#include "itkMABMISSimulateData.h" -#include "itkMABMISImageRegistrationFilter.h" -#include "itkMABMISTreeOperation.h" -#include "itkMABMISBasicOperationFilter.h" - -#include "itkMABMISAtlasXMLFile.h" - -static std::string ReplacePathSepForOS( const std::string & input ) -{ - std::string output = input; -#ifdef _WIN32 - std::replace(output.begin(), output.end(), '/', FILESEP); -#else - std::replace(output.begin(), output.end(), '\\', FILESEP); -#endif - return output; -} - -typedef double CoordinateRepType; -const unsigned int SpaceDimension = ImageDimension; - -// basic data type -typedef unsigned char CharPixelType; // for image IO usage -typedef float FloatPixelType; // for -typedef int IntPixelType; -typedef short ShortPixelType; -typedef float InternalPixelType; // for internal processing usage -typedef itk::Vector VectorPixelType; - -// basic image type -typedef itk::Image CharImageType; -typedef itk::Image IntImageType; -typedef itk::Image ShortImageType; -typedef itk::Image FloatImageType; -typedef itk::Image InternalImageType; -typedef itk::Image DeformationFieldType; - -// basic iterator type -typedef itk::ImageRegionIterator DeformationFieldIteratorType; -typedef itk::ImageRegionIterator InternalImageIteratorType; -typedef itk::ImageRegionIterator CharImageIteratorType; - -// basic image reader/writer related type -typedef itk::ImageFileReader CharImageReaderType; -typedef itk::ImageFileReader InternalImageReaderType; -typedef itk::ImageFileWriter InternalImageWriterType; - -typedef itk::WarpImageFilter InternalWarpFilterType; -typedef itk::ImageFileWriter CharImageWriterType; -typedef itk::ImageFileWriter IntImageWriterType; -typedef itk::ImageFileWriter FloatImageWriterType; -typedef itk::ImageFileWriter ShortImageWriterType; - -typedef itk::ImageFileReader DeformationFieldReaderType; -typedef itk::ImageFileWriter DeformationFieldWriterType; - -////////////////////////////////////////////////////////////////////////////// -// image filter type -typedef itk::ResampleImageFilter ResampleFilterType; -typedef itk::HistogramMatchingImageFilter InternalHistMatchFilterType; - -//////////////////////////////////////////////////////////////////////////// -// operation on deformation fields -typedef itk::WarpVectorImageFilter WarpVectorFilterType; -typedef itk::InverseDeformationFieldImageFilter - InverseDeformationFieldImageFilterType; -typedef itk::AddImageFilter AddImageFilterType; - -// global bool variables to adjust the procedure - -bool isEvaluate = false; // if false, we do not know the ground-truth of labels -bool isDebug = false; // false;//true; // if true, print out more information - -// ----------------------------------------------------------------------------- -// global variables - -// int localPatchSize = 1; //(2r+1)*(2r+1)*(2r+1) is the volume of local patch - -// demons registration parameters -// int iterInResolutions[4][3]={{5,3,2},{10,5,5},{15,10,5},{20,15,10}}; -// int itereach = 2; // -// int itereach0 = 0;int itereach1 = 1;int itereach2 = 2;int itereach3 = 3; -// double sigmaDef = 1.5; -// double sigmaDef10 = 1.0;double sigmaDef15 = 1.5;double sigmaDef20 = 2.0; -// double sigmaDef25 = 2.5;double sigmaDef30 = 3.0;double sigmaDef35 = 3.5; -bool doHistMatch = true; -bool doHistMatchTrue = true; bool doHistMatchFalse = false; - -// ----------------------------------------------------------------------------- - -typedef itk::Statistics::MABMISSimulateData DataSimulatorType; -DataSimulatorType::Pointer datasimulator; - -typedef itk::Statistics::MABMISImageOperationFilter ImageOperationFilterType; -ImageOperationFilterType::Pointer imgoperator; -typedef itk::Statistics::MABMISDeformationFieldFilter DeformationFieldOperationFilterType; -DeformationFieldOperationFilterType::Pointer dfoperator; -typedef itk::Statistics::MABMISImageRegistrationFilter ImageRegistrationFilterType; -ImageRegistrationFilterType::Pointer regoperator; -typedef itk::Statistics::MABMISTreeOperation TreeOperationType; -TreeOperationType::Pointer treeoperator; - -typedef itk::Statistics::MABMISBasicOperationFilter BasicOperationFilterType; -BasicOperationFilterType::Pointer basicoperator; - -// -DeformationFieldType::SpacingType df_spacing; -DeformationFieldType::DirectionType df_direction; -DeformationFieldType::PointType df_origin; - -typedef itk::Vector ShortVectorPixelType; -typedef itk::Image ShortDeformationFieldType; -typedef itk::ImageFileWriter ShortDeformationFieldWriterType; - -void SearchRootAmongAtlases(std::vector imgfilenames, int & root); - -int RegistrationBetweenRootandAtlases(int root, std::vector imageFileNames, std::vector iterations, - double sigma); - -int BuildStatisticalDeformationModel(int root, std::vector imageFileNames, int simulatedAtlasSize); - -std::vector GenerateSimulatedData(int root, std::vector imgfilenames, int simulatedAtlasSize); - -void strtrim(std::string& str) -{ - std::string::size_type pos = str.find_last_not_of(' '); - - if( pos != std::string::npos ) - { - str.erase(pos + 1); - pos = str.find_first_not_of(' '); - if( pos != std::string::npos ) - { - str.erase(0, pos); - } - } - else - { - str.erase(str.begin(), str.end() ); - } -} - -int Training( itk::MABMISImageData* trainingData, std::string outputFile, - std::vector iterations, double sigma) -{ - // sanity check - std::cout << "m_DataDirectory: " << trainingData->m_DataDirectory << std::endl; - if( trainingData->m_SegmentationFileNames.size() != trainingData->m_ImageFileNames.size() ) - { - std::cerr << "The numbers of image files and segmentation files are NOT equal!!!" << std::endl; - return -1; - } - if( trainingData->m_ImageFileNames.size() == 0 ) - { - std::cerr << "No Atlas is provided. QUIT. " << std::endl; - return -1; - } - if( trainingData->m_NumberImageData != trainingData->m_ImageFileNames.size() ) - { - trainingData->m_NumberImageData = trainingData->m_ImageFileNames.size(); - } - - // step 0: files and folders - // get output path from the output file name - size_t sep = outputFile.find_last_of(FILESEP); - std::string outputFolder = ""; - if( sep != std::string::npos ) - { - outputFolder = outputFile.substr(0, sep) + FILESEP; - } - else - { - sep = -1; - } - // get output file name without extension, and use it to create a directory to save trained atlas - size_t fnSep = outputFile.find_last_of('.'); - std::string outputxmlname = ""; - if( fnSep != std::string::npos ) - { - outputxmlname = outputFile.substr(sep + 1, fnSep - sep - 1); - } - else - { - outputxmlname = outputFile.substr(sep + 1, std::string::npos); - } - - std::string atlasFolder; - if( outputFolder.empty() ) - { - atlasFolder = outputxmlname; - } - else - { - atlasFolder = outputFolder + outputxmlname; - } - atlasFolder = ReplacePathSepForOS(atlasFolder); - if( !atlasFolder.empty() ) - { - atlasFolder = atlasFolder + FILESEP; - } - itksys::SystemTools::MakeDirectory(atlasFolder.c_str() ); - - // copy training data into the atlas folder - std::vector imageFiles(trainingData->m_NumberImageData); - std::vector segmentFiles(trainingData->m_NumberImageData); - for( int i = 0; i < imageFiles.size(); ++i ) - { - imageFiles[i] = ReplacePathSepForOS(atlasFolder + trainingData->m_ImageFileNames[i]); - segmentFiles[i] = ReplacePathSepForOS(atlasFolder + trainingData->m_SegmentationFileNames[i]); - - // copy - std::string fromImagefile = trainingData->m_ImageFileNames[i]; - if( !trainingData->m_DataDirectory.empty() ) - { - fromImagefile = ReplacePathSepForOS(trainingData->m_DataDirectory + trainingData->m_ImageFileNames[i]); - } - - bool ret1=itksys::SystemTools::CopyFileAlways(fromImagefile.c_str(), imageFiles[i].c_str() ); - if (!ret1) - { - std::cerr << "ERROR: Cannot copy atlas image file to trained atlas folder!" << std::endl; - return -1; - } - - // if file extension is '.nii.gz', also copy the 'img file' - if( itksys::SystemTools::Strucmp(fromImagefile.substr(fromImagefile.size() - 3, 3).c_str(), "hdr") == 0 ) - { - std::string fromfileImg = fromImagefile.substr(0, fromImagefile.size() - 3) + "img"; - std::string tofileImg = imageFiles[i].substr(0, imageFiles[i].size() - 3) + "img"; - itksys::SystemTools::CopyFileAlways(fromfileImg.c_str(), tofileImg.c_str() ); - } - - std::string fromSegfile = trainingData->m_DataDirectory + trainingData->m_SegmentationFileNames[i]; - - bool ret2 = itksys::SystemTools::CopyFileAlways(fromSegfile.c_str(), segmentFiles[i].c_str() ); - if (!ret1) - { - std::cerr << "ERROR: Cannot copy atlas image file to trained atlas folder!" << std::endl; - return -1; - } - // if file extension is '.nii.gz', also copy the 'img file' - if( itksys::SystemTools::Strucmp(fromSegfile.substr(fromSegfile.size() - 3, 3).c_str(), "hdr") == 0 ) - { - std::string fromfileImg = fromSegfile.substr(0, fromSegfile.size() - 3) + "img"; - std::string tofileImg = segmentFiles[i].substr(0, segmentFiles[i].size() - 3) + "img"; - itksys::SystemTools::CopyFileAlways(fromfileImg.c_str(), tofileImg.c_str() ); - } - } - - // step 1: find root of the atlas - std::cout << "---------------------------------------" << std::endl; - std::cout << "1. Find root atlas ... " << std::endl;; - //////////////////////////////// - // build tree on atlases to find the root - int root = -1; // root node of tree - SearchRootAmongAtlases(imageFiles, root); - std::cout << "Done." << std::endl;; - - // step 2.1: register all atlas to the root - ////////////////////////////// - // registration between root and other atlases - std::cout << "---------------------------------------" << std::endl; - std::cout << "2. Run coarse registration ... " << std::endl; - int val = RegistrationBetweenRootandAtlases(root, imageFiles, iterations, sigma); - if( val != 0 ) - { - return -1; - } - std::cout << "Done!" << std::endl; - - int simulatedAtlasSize = 2 * imageFiles.size(); - - std::cout << "---------------------------------------" << std::endl; - std::cout << "3. Build statistical deformation model..." << std::endl; - if( BuildStatisticalDeformationModel(root, imageFiles, simulatedAtlasSize) != 0 ) - { - return -1; - } - std::cout << "Done. " << std::endl; - - // step 2.5: generate templaet images - ///////////////////////// - // generate simulated template with deformation field - std::cout << "---------------------------------------" << std::endl; - std::cout << "4. Generate simulated templates ... " << std::endl; - std::vector simulatedImageFiles = GenerateSimulatedData(root, imageFiles, simulatedAtlasSize); - std::cout << "Done. " << std::endl; - - ///////////////////////// - // delete subsampled deformation field - // if (!isDebug) if (false) - { - for( int i = 0; i < imageFiles.size(); ++i ) - { - if( i == root ) - { - continue; - } - char i_str[10], root_str[10]; - sprintf(i_str, "%03d", i); - sprintf(root_str, "%03d", root); - std::string subdeformationFieldFileName = std::string(i_str) + "_to_" + std::string(root_str) - + "deform_000_sub.nii.gz"; - subdeformationFieldFileName = outputFolder + subdeformationFieldFileName; - basicoperator->RemoveFile(subdeformationFieldFileName); - } - } - - ////////////////////////////// - // build the combinative tree - std::cout << "---------------------------------------" << std::endl; - std::cout << "5. Build combinative tree ... " << std::endl; - int totalAtlasSize = imageFiles.size() + simulatedAtlasSize; - int tree_size = totalAtlasSize; - vnl_vector tree(tree_size); // each tree - - // fill distance matrix - std::vector atlasfilenames; - for( int i = 0; i < tree_size; ++i ) - { - std::string atlasfilename; - - if( i < imageFiles.size() ) - { - atlasfilename.append(imageFiles[i]); - } - else - { - atlasfilename.append(simulatedImageFiles[i - imageFiles.size()]); - } - - // atlasfilename = outputFolder + atlasfilename; - atlasfilenames.push_back(atlasfilename); - } - - tree = treeoperator->BuildCombinativeTree(root, atlasfilenames, tree_size, tree); // cout << "pass: - // BuildCombinativeTree " << endl; - - std::cout << "Done. " << std::endl; - - // write trained atlas to xml file - std::cout << "---------------------------------------" << std::endl; - std::cout << "6. Write to output file..." << std::endl; - itk::MABMISAtlasXMLFileWriter::Pointer atlasWriter = itk::MABMISAtlasXMLFileWriter::New(); - itk::MABMISAtlas atlas; - atlas.m_NumberAllAtlases = totalAtlasSize; - atlas.m_NumberSimulatedAtlases = simulatedAtlasSize; - atlas.m_AtlasDirectory = std::string(".") + FILESEP + outputxmlname; - atlas.m_Tree.resize(tree_size); - atlas.m_IsSimulatedImage.resize(tree_size); - atlas.m_RealImageIDs.resize(totalAtlasSize - simulatedAtlasSize); - atlas.m_SimulatedImageIDs.resize(simulatedAtlasSize); - for( int i = 0; i < tree_size; ++i ) - { - atlas.m_Tree[i] = tree[i]; - } - atlas.m_AtlasFilenames.resize(atlasfilenames.size() ); - int countRealImages = 0, countSimulatedImages = 0; - for( int i = 0; i < atlasfilenames.size(); ++i ) - { - const size_t sep = atlasfilenames[i].find_last_of(FILESEP); - std::string fname = atlasfilenames[i]; - if( sep != std::string::npos ) - { - fname = atlasfilenames[i].substr(sep + 1, std::string::npos); -/* fname = atlasfilenames[i].substr(0, sep); - sep = fname.find_last_of(FILESEP); - if (sep != std::string::npos) - fname = atlasfilenames[i].substr(sep+1, std::string::npos); - else - fname = atlasfilenames[i]; - */ - } - else - { - fname = atlasfilenames[i]; - } - - atlas.m_AtlasFilenames[i] = fname; - - if( i < totalAtlasSize - simulatedAtlasSize ) - { - atlas.m_IsSimulatedImage[i] = false; - atlas.m_RealImageIDs[countRealImages++] = i; - } - else - { - atlas.m_IsSimulatedImage[i] = true; - atlas.m_SimulatedImageIDs[countSimulatedImages++] = i; - } - } - // atlas segmentation files - atlas.m_AtlasSegmentationFilenames.resize(segmentFiles.size() ); - for( int i = 0; i < segmentFiles.size(); ++i ) - { - const size_t sep = segmentFiles[i].find_last_of(FILESEP); - std::string fname = segmentFiles[i]; - if( sep != std::string::npos ) - { - fname = segmentFiles[i].substr(sep + 1, std::string::npos); -/* fname = segmentFiles[i].substr(0, sep); - sep = fname.find_last_of(FILESEP); - if (sep != std::string::npos) - fname = segmentFiles[i].substr(sep+1, std::string::npos); - else - fname = segmentFiles[i]; -*/ - } - else - { - fname = segmentFiles[i]; - } - atlas.m_AtlasSegmentationFilenames[i] = fname; - } - - treeoperator->FindRoot(tree, tree_size, atlas.m_TreeRoot); - treeoperator->GetTreeHeight(tree, tree_size, atlas.m_TreeHeight); - - atlasWriter->SetObject(&atlas); - atlasWriter->SetFilename(outputFile); - atlasWriter->WriteFile(); - - std::cout << "Done." << std::endl; - - return EXIT_SUCCESS; -} - -template -int DoIt( itk::MABMISImageData* trainingData, std::string outputFile, - std::vector iterations, double sigma) -{ - return Training(trainingData, outputFile, iterations, sigma); -} +#include "Training.h" int main( int argc, char *argv[] ) { @@ -544,260 +54,8 @@ int main( int argc, char *argv[] ) trainingData->m_DataDirectory = trainingData->m_DataDirectory + FILESEP; } } - // Will look into getting rid of these global variables later ---Xiaofeng - datasimulator = DataSimulatorType::New(); - imgoperator = ImageOperationFilterType::New(); - dfoperator = DeformationFieldOperationFilterType::New(); - regoperator = ImageRegistrationFilterType::New(); - treeoperator = TreeOperationType::New(); - basicoperator = BasicOperationFilterType::New(); - - // get number of iterations - int retVal = DoIt(trainingData, TrainingOutputFile, iterations, SmoothingKernelSize); + int retVal = Training(trainingData, TrainingOutputFile, iterations, SmoothingKernelSize); delete trainingData; - return retVal; - - return EXIT_SUCCESS; -} - -void SearchRootAmongAtlases(std::vector imgfilenames, int & root) -{ - int tree_size = imgfilenames.size(); - - vnl_matrix distanceMatrix(tree_size, tree_size); - for( int i = 0; i < tree_size; ++i ) - { - for( int j = 0; j < tree_size; ++j ) - { - distanceMatrix[i][j] = 0.0; - } - } - - vnl_vector tree(tree_size); // each tree - - // fill distance matrix - imgoperator->PairwiseDistanceAmongImages(imgfilenames, tree_size, distanceMatrix); - if( isDebug ) - { - std::string distFileName = "dist_all_atlases.txt"; - basicoperator->SaveMatrix2File(distanceMatrix, tree_size, tree_size, distFileName); - } - - // build MST - tree = treeoperator->generateMSTFromMatrix(distanceMatrix, tree_size, tree); - - treeoperator->FindRoot(tree, tree_size, root); - - // save distance matrix and tree - if( isDebug ) - { - std::string treeFileName = "tree_all_atlases.txt"; - treeoperator->SaveTreeWithInfo(tree, tree_size, treeFileName); - } - - std::cout << "The root is " << root << ": " << imgfilenames[root] << ". " << std::endl; -} - -int RegistrationBetweenRootandAtlases(int root, std::vector imageFileNames, - std::vector iterations, double sigma) -{ - // get the path from the filename - const size_t sep = imageFileNames[0].find_last_of(FILESEP); - std::string outputFolder = ""; - - if( sep != std::string::npos ) - { - outputFolder = imageFileNames[0].substr(0, sep); - } - if( !outputFolder.empty() ) - { - outputFolder = outputFolder + FILESEP; - } - - const int atlas_size = imageFileNames.size(); - std::cout << "Register between root and atlases ... " << std::endl; - for( int i = 0; i < atlas_size; ++i ) - { - std::cout << i << ", "; - if( i == root ) - { - continue; - } - // if not the root image do registration - std::string fixedImageFileName = imageFileNames[root]; - std::string movingImageFileName = imageFileNames[i]; - - char i_str[10], root_str[10]; - sprintf(i_str, "%03d", i); - sprintf(root_str, "%03d", root); - std::string deformedImageFileName = std::string(i_str) + "_to_" + std::string(root_str) + "_cbq_reg.nii.gz"; - std::string deformationFieldFileName = std::string(i_str) + "_to_" + std::string(root_str) + "_deform_000.nii.gz"; - - deformedImageFileName = outputFolder + deformedImageFileName; - deformationFieldFileName = outputFolder + deformationFieldFileName; - - if( !itksys::SystemTools::FileExists(deformationFieldFileName.c_str(), true) ) - { - // rough registration - - int val = regoperator->DiffeoDemonsRegistrationWithParameters( - fixedImageFileName, movingImageFileName, - deformedImageFileName, deformationFieldFileName, - sigma + 0.5, doHistMatchTrue, iterations); - if( val != 0 ) - { - std::cout << "Cannot perform registration between :" << std::endl; - std::cout << fixedImageFileName << " and " << std::endl; - std::cout << movingImageFileName << "." << std::endl; - std::cout << "Please verify the input images are correct." << std::endl; - } - } - - // --------------------------------------------------------------------- - // ------------- Generate the inverse deformation field - - DeformationFieldType::Pointer deformationField = 0; - dfoperator->ReadDeformationField(deformationFieldFileName, deformationField); - - DeformationFieldType::Pointer inversedDeformationField = DeformationFieldType::New(); - inversedDeformationField->SetRegions(deformationField->GetLargestPossibleRegion() ); - inversedDeformationField->SetSpacing(deformationField->GetSpacing() ); - inversedDeformationField->SetDirection(deformationField->GetDirection() ); - inversedDeformationField->SetOrigin(deformationField->GetOrigin() ); - inversedDeformationField->Allocate(); - - dfoperator->InverseDeformationField3D(deformationField, inversedDeformationField); - - std::string invDeformedImageFileName = std::string(root_str) + "_to_" + std::string(i_str) + "_cbq_reg.nii.gz"; - std::string invDeformationFileName = std::string(root_str) + "_to_" + std::string(i_str) + "_deform_000.nii.gz"; - invDeformedImageFileName = outputFolder + invDeformedImageFileName; - invDeformationFileName = outputFolder + invDeformationFileName; - - std::string invDeformedSegmentFileName = std::string(root_str) + "_to_" + std::string(i_str) + "seg_000.nii.gz"; - invDeformedSegmentFileName = outputFolder + invDeformedSegmentFileName; - - dfoperator->WriteDeformationField(invDeformationFileName, inversedDeformationField); - - // update - regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( - movingImageFileName, fixedImageFileName, - invDeformationFileName, - invDeformedImageFileName, invDeformationFileName, - sigma, doHistMatch, iterations); - - if( isEvaluate ) - { - dfoperator->ApplyDeformationFieldAndWriteWithFileNames( - fixedImageFileName, - invDeformationFileName, - invDeformedSegmentFileName, false); - } - } - std::cout << std::endl; - std::cout << "Done. " << std::endl; - - return 0; } - -int BuildStatisticalDeformationModel(int root, std::vector imageFileNames, int simulatedAtlasSize) -{ - /////////////////////////////// - // do PCA simulation - // std::cout << "Build deformation field model ... "; - std::vector allDeformationFieldFileNames; - - // get the path from the filename - const size_t sep = imageFileNames[0].find_last_of(FILESEP); - std::string outputFolder = ""; - if( sep != std::string::npos ) - { - outputFolder = imageFileNames[0].substr(0, sep); - } - if( !outputFolder.empty() ) - { - outputFolder = outputFolder + FILESEP; - } - - int atlas_size = imageFileNames.size(); - - datasimulator->SetRoot(root); - datasimulator->SetAtlasSize(atlas_size); - datasimulator->SetSimulateSize(simulatedAtlasSize); - for( int i = 0; i < atlas_size; ++i ) - { - if( i == root ) - { - continue; - } - - char i_str[10], root_str[10]; - sprintf(i_str, "%03d", i); - sprintf(root_str, "%03d", root); - - std::string deformationFieldFileName = std::string(i_str) + "_to_" + std::string(root_str) + "_deform_000.nii.gz"; - deformationFieldFileName = outputFolder + deformationFieldFileName; - allDeformationFieldFileNames.push_back(deformationFieldFileName); - } - - return datasimulator->DoPCATraining(allDeformationFieldFileNames, atlas_size - 1, imageFileNames, root); -} - -// return: simulated template file names -std::vector GenerateSimulatedData(int root, std::vector imageFiles, int simulatedAtlasSize) -{ - // output folder - const size_t sep = imageFiles[0].find_last_of(FILESEP); - std::string outputFolder = ""; - - if( sep != std::string::npos ) - { - outputFolder = imageFiles[0].substr(0, sep); - } - if( !outputFolder.empty() ) - { - outputFolder = outputFolder + FILESEP; - } - - std::vector simulateDeformationFieldFileNames(0); - std::vector simulateTemplateFileNames(0); - std::cout << "Generating simulated images: " << std::endl; - for( int i = 0; i < simulatedAtlasSize; ++i ) - { - std::cout << i << ", "; - - std::string index_string; - basicoperator->myitoa( i, index_string, 3 ); - - std::string simulateDeformationFieldFileName = "simulated_deform_" + index_string + ".nii.gz"; - std::string simulateTemplateFileName = "simulated_cbq_" + index_string + ".nii.gz"; - simulateDeformationFieldFileName = outputFolder + simulateDeformationFieldFileName; - simulateTemplateFileName = outputFolder + simulateTemplateFileName; - - simulateDeformationFieldFileNames.push_back(simulateDeformationFieldFileName); - simulateTemplateFileNames.push_back(simulateTemplateFileName); - - // load simulated deformation field - DeformationFieldType::Pointer df = 0; - dfoperator->ReadDeformationField(simulateDeformationFieldFileNames[i], df); - - DeformationFieldType::Pointer invdf = DeformationFieldType::New(); - invdf->SetRegions(df->GetLargestPossibleRegion() ); - invdf->SetSpacing(df->GetSpacing() ); - invdf->SetDirection(df->GetDirection() ); - invdf->SetOrigin(df->GetOrigin() ); - invdf->Allocate(); - - dfoperator->InverseDeformationField3D(df, invdf); - InternalImageType::Pointer rootImg = 0; - imgoperator->ReadImage(imageFiles[root], rootImg); - InternalImageType::Pointer simImg = 0; - dfoperator->ApplyDeformationField(rootImg, invdf, simImg, true); - imgoperator->WriteImage(simulateTemplateFileNames[i], simImg); - } - std::cout << std::endl; - - return simulateTemplateFileNames; -} - diff --git a/Testing.cxx b/Testing.cxx new file mode 100644 index 0000000..c37277a --- /dev/null +++ b/Testing.cxx @@ -0,0 +1,1262 @@ +/*========================================================================= + Copyright (c) IDEA LAB, UNC-Chapel Hill, 2013. + + MABMIS (Multi-Atlas-Based Multi-Image Segmentation) + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#include "Testing.h" + +static std::string +GetDefaultSegmentationFilename(itk::MABMISImageData* imageData, size_t imageIndex, std::string iterationString) +{ + std::string prefix = iterationString; + if (!iterationString.empty()) + { + prefix += "_"; + } + + if (imageData->m_SegmentationFileNames.size() >= imageIndex && !imageData->m_SegmentationFileNames[imageIndex].empty()) + { + return imageData->m_OutputDirectory + prefix + imageData->m_SegmentationFileNames[imageIndex]; + } + else + { + std::string imageFileName = imageData->m_OutputDirectory + prefix + imageData->m_ImageFileNames[imageIndex]; + return GetRootName(imageFileName) + "-label" + GetExtension(imageFileName); + } +} + + +void HistogramMatching(InternalImageType::Pointer inputImage, InternalImageType::Pointer referenceImage, + InternalImageType::Pointer & outputImage); + +void TreeBasedRegistrationFastOniTree(vnl_vector itree, // the incremental tree + int root, // the tree root + int itree_size, // the tree size + itk::MABMISAtlas* atlasTree, // the atlas tree + itk::MABMISImageData*imageData, // the images to be segmented + std::vector iterations, // registration parameters + double sigma); + +void LabelFusion(std::string curSampleImgName, std::string outSampleSegName, std::vector allWarpedAtlasImgNames, + std::vector allDeformationFieldNames, std::vector allAtlasSegNames, int numOfAtlases); + +void GaussianWeightedLabelFusion(InternalImageType::Pointer curSampleImgPtr, InternalImageType::Pointer outSampleSegPtr, + InternalImageType::Pointer* warpedImgPtrs, InternalImageType::Pointer* warpedSegPtrs, + unsigned numOfAtlases); + +void RegistrationOntoTreeRoot(vnl_vector itree, // the incremental tree + int root, // the tree root + int itree_size, // the tree size + itk::MABMISAtlas* atlasTree, // the atlas tree + itk::MABMISImageData*imageData, // the images to be segmented + std::vector iterations, // registration parameters + double sigma); + +// void PairwiseRegistrationOnTreeViaRoot(int root, int filenumber, int atlas_size, std::vector sub_ids); +void PairwiseRegistrationOnTreeViaRoot(int root, itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree, + std::vector iterations, double sigma); + +// void MultiAtlasBasedSegmentation(int filenumber, int atlas_size, std::vector sub_ids); +int MultiAtlasBasedSegmentation(itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree); + +int Testing(itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree, + std::vector iterations, double sigma) +{ + // Validate the tree size is correct + if( atlasTree->m_TreeSize != atlasTree->m_Tree.size() || + atlasTree->m_TreeSize != atlasTree->m_NumberAllAtlases ) + { + std::cerr << "ERROR: The atlas tree size is incorrect!!" << std::endl; + return -1; + } + ///////////////////////////////////// + // the atlas tree + unsigned tree_size = atlasTree->m_NumberAllAtlases; + + // incremental tree: the atlas tree PLUS images to be segmented + unsigned itree_size = atlasTree->m_NumberAllAtlases + imageData->m_NumberImageData; + vnl_vector itree(itree_size); + for( unsigned i = 0; i < itree_size; ++i ) + { + itree.put(i, 0); + } + // copy current tree into an incremental tree + for( unsigned i = 0; i < tree_size; ++i ) + { + itree.put(i, atlasTree->m_Tree[i]); // combinative tree: all atlases including + // simulated data + } + treeoperator->FindRoot(itree, atlasTree->m_TreeSize, atlasTree->m_TreeRoot); + treeoperator->GetTreeHeight(itree, atlasTree->m_TreeSize, atlasTree->m_TreeHeight); + int root = atlasTree->m_TreeRoot; // root node + // Validate the tree is correct + for( unsigned n = 0; n < atlasTree->m_TreeSize; ++n ) + { + bool valid = true; + if( atlasTree->m_Tree[n] >= int(atlasTree->m_TreeSize) ) + { + valid = false; + } + if( atlasTree->m_Tree[n] < 0 && int(n) != atlasTree->m_TreeRoot ) + { + valid = false; + } + if( !valid ) + { + std::cerr << "ERROR: The atlas tree size is incorrect!!" << std::endl; + return -1; + } + } + + if( isDebug ) + { + treeoperator->SaveTreeWithInfo(itree, itree_size, imageData->m_OutputDirectory + "itree.txt"); + } + + // build the incremental tree + std::cout << "---------------------------------------" << std::endl; + std::cout << "Build incremental tree ... " << std::endl; + + unsigned totalNumAtlases = atlasTree->m_NumberAllAtlases; + unsigned totalNumFiles = totalNumAtlases + imageData->m_NumberImageData; + std::vector allfilenames(totalNumFiles); + for( unsigned i = 0; i < totalNumFiles; ++i ) + { + + if( i < totalNumAtlases ) + { + allfilenames[i] = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[i]); + } + else + { + allfilenames[i] = ReplacePathSepForOS(imageData->m_DataDirectory + imageData->m_ImageFileNames[i - totalNumAtlases]); + } + } + + std::cout << "---------------------------------------" << std::endl; + std::cout << "1. Calculating pairwise distance cross training and test images ..." << std::endl; + vnl_matrix cross_dist(imageData->m_NumberImageData, totalNumFiles); + for( unsigned i = 0; i < imageData->m_NumberImageData; ++i ) + { + for( unsigned j = 0; j < totalNumFiles; ++j ) + { + cross_dist.put(i, j, imgoperator->calculateDistanceMSD(allfilenames[totalNumAtlases + i], allfilenames[j]) ); + } + } + + if( isDebug ) + { + basicoperator->SaveMatrix2File(cross_dist, imageData->m_NumberImageData, totalNumFiles, + imageData->m_OutputDirectory + "sample2atlas_dist.txt"); + } + + // indices to for atlas files and image files to be segmented. + int* atlas_cur = new int[totalNumFiles]; + for( unsigned i = 0; i < totalNumFiles; ++i ) + { + atlas_cur[i] = 0; + } + for( unsigned i = 0; i < totalNumAtlases; ++i ) + { + atlas_cur[i] = i; + } + int* images_cur = new int[imageData->m_NumberImageData]; + for( unsigned i = 0; i < imageData->m_NumberImageData; ++i ) + { + images_cur[i] = i + totalNumAtlases; + } + + // build the incremental tree + treeoperator->SetAllAtlasNumber(totalNumAtlases); + std::cout << "---------------------------------------" << std::endl; + std::cout << "2. Build the tree... " << std::endl; + itree = treeoperator->BuildIncrementalTree(imageData->m_NumberImageData, // #images to be segmented + totalNumAtlases, // #atlases, including original ones and simulated ones + totalNumFiles, // total #image files, = imageData->m_NumberImageData+totalNumAtlases + images_cur, // file indices for images to be segmented + atlas_cur, // file indices for atlas images + cross_dist, // pairwise distance between images + itree); // the output incremental tree + + // delete + delete[] atlas_cur; + delete[] images_cur; + + // check if new tree is actually a tree + const bool iTreeValid = treeoperator->ValidateTree(itree, itree_size); + if( !iTreeValid ) + { + std::cout << "The new built tree is NOT valid!" << std::endl; + + return 105; // EXIT_FAILURE; + } + // save new tree + if( isDebug ) + { + treeoperator->SaveTreeWithInfo(itree, itree_size, imageData->m_OutputDirectory + "itree.txt"); + // std::cout << "Tree is built up!" << std::endl; + } + std::cout << "Done! " << std::endl; + + std::cout << "---------------------------------------" << std::endl; + std::cout << "3. Register all images to the root... " << std::endl; + // register all images to the root, if not done yet. + RegistrationOntoTreeRoot(itree, // the incremental tree + root, // the tree root + itree_size, // the tree size + atlasTree, + imageData, + iterations, sigma // registration parameters + ); + std::cout << "Done. " << std::endl; + + std::cout << "---------------------------------------" << std::endl; + std::cout << "4. Pairwise registration via the tree root... " << std::endl; + PairwiseRegistrationOnTreeViaRoot(root, imageData, atlasTree, iterations, sigma); + std::cout << "Done. " << std::endl; + + //////////////////////////////// + // do all multi-atlas based segmentation with weighted label fusion + std::cout << "---------------------------------------" << std::endl; + std::cout << "5. Segment images using label fusion... " << std::endl; + int numIter = MultiAtlasBasedSegmentation(imageData, atlasTree); + + std::cout << "---------------------------------------" << std::endl; + std::cout << "Done. " << std::endl; + + // do all multi-atlas based segmentation with weighted label fusion + ///////////////////////////////// + // + // basicoperator->RemoveFile("*_*_cbq_000*"); + // basicoperator->RemoveFile("*_*_deform_000*"); + // + ///////////////////////////////// + + // remove intermediate results + // deformedImageFileName and deformation fields + if( !isDebug ) + { + std::cout << "---------------------------------------" << std::endl; + std::cout << "Removing intermediate registration files..." << std::endl; + for( unsigned n = 0; n < imageData->m_NumberImageData; ++n ) + { + char n_str[10], m_str[10]; + sprintf(n_str, "%03d", n); + for( int m = 0; m < int(atlasTree->m_NumberAllAtlases) - int(atlasTree->m_NumberSimulatedAtlases); ++m ) + { + sprintf(m_str, "%03d", m); + + std::string imHdr = std::string(m_str) + "_to_" + "testImage" + n_str + ImageSuffix; + std::string dfMha = std::string(m_str) + "_to_" + "testImage" + n_str + DeformationSuffix; + imHdr = imageData->m_OutputDirectory + imHdr; + dfMha = imageData->m_OutputDirectory + dfMha; + basicoperator->RemoveFile(imHdr.c_str() ); + basicoperator->RemoveFile(dfMha.c_str() ); + + if( m == root ) + { + imHdr = std::string("testImage") + n_str + "_to_" + m_str + RegisteredSuffix; + dfMha = std::string("testImage") + n_str + "_to_" + m_str + DeformationSuffix; + imHdr = imageData->m_OutputDirectory + imHdr; + dfMha = imageData->m_OutputDirectory + dfMha; + basicoperator->RemoveFile(imHdr.c_str() ); + basicoperator->RemoveFile(dfMha.c_str() ); + } + } + for( unsigned m = 0; m < imageData->m_NumberImageData; ++m ) + { + if( m == n ) + { + continue; + } + sprintf(m_str, "%03d", m); + std::string imHdr = std::string("testImage") + n_str + "_to_testImage" + m_str + ImageSuffix; + std::string dfMha = std::string("testImage") + n_str + "_to_testImage" + m_str + DeformationSuffix; + imHdr = imageData->m_OutputDirectory + imHdr; + dfMha = imageData->m_OutputDirectory + dfMha; + basicoperator->RemoveFile(imHdr.c_str() ); + basicoperator->RemoveFile(dfMha.c_str() ); + } + } + } + // remove segmentation results in iterations. + for( unsigned n = 0; n < imageData->m_NumberImageData; ++n ) + { + const std::string imageFileName = imageData->m_OutputDirectory + imageData->m_ImageFileNames[n]; + char i_str[10]; + bool copied = false; + for( int i = numIter - 1; i >= 0; i-- ) + { + sprintf(i_str, "%03d", i); + const std::string segHdr = GetDefaultSegmentationFilename(imageData, n, i_str); + if( itksys::SystemTools::FileExists(segHdr.c_str(), true) ) + { + if( !copied ) + { + std::string segHdr_save = GetDefaultSegmentationFilename(imageData, n, ""); + itksys::SystemTools::CopyFileAlways(segHdr.c_str(), segHdr_save.c_str()); + std::cout << "Segmentation " << i_str << ": " + segHdr_save << std::endl; + copied = true; + } + basicoperator->RemoveFile(segHdr.c_str() ); + } + } + } + + std::cout << "Done. " << std::endl; + return EXIT_SUCCESS; +} + +void LabelFusion(std::string curSampleImgName, std::string outSampleSegName, std::vector allWarpedAtlasImgNames, + std::vector allDeformationFieldNames, std::vector allAtlasSegNames, int numOfAtlases) +{ + // create pointers for each images including: cur_image, our_seg, + // warpedImg(after histogram matching), warpedSeg (after warping) + + // cur sample image pointer + InternalImageType::Pointer curSampleImgPtr = nullptr; + + imgoperator->ReadImage(curSampleImgName, curSampleImgPtr); + const itk::ImageIOBase::IOComponentType ioType = imgoperator->GetIOPixelType(curSampleImgName); + + // output sample label image pointer + InternalImageType::Pointer outSampleSegPtr = nullptr; + imgoperator->ReadImage(curSampleImgName, outSampleSegPtr); + + // atlas related images pointers + InternalImageType::Pointer* warpedImgPtrs = new InternalImageType::Pointer[numOfAtlases]; + InternalImageType::Pointer* warpedSegPtrs = new InternalImageType::Pointer[numOfAtlases]; + InternalImageType::IndexType index; + index[0] = 19; index[1] = 19; index[2] = 8; + for( int i = 0; i < numOfAtlases; ++i ) + { + // histogram match each warped atlas to the current sample + warpedImgPtrs[i] = 0; + InternalImageType::Pointer curWarpedAtlasPtr = nullptr; + imgoperator->ReadImage(allWarpedAtlasImgNames[i], curWarpedAtlasPtr); + HistogramMatching(curWarpedAtlasPtr, curSampleImgPtr, warpedImgPtrs[i]); + + // load each deformation field + DeformationFieldType::Pointer deformFieldPtr = nullptr; + dfoperator->ReadDeformationField(allDeformationFieldNames[i], deformFieldPtr); + + // warp each atlas label image to the current sample space + InternalImageType::Pointer curAtlasSegPtr = nullptr; + imgoperator->ReadImage(allAtlasSegNames[i], curAtlasSegPtr); + dfoperator->ApplyDeformationField(curAtlasSegPtr, deformFieldPtr, warpedSegPtrs[i], false); + } + + // weighted label fusion + GaussianWeightedLabelFusion(curSampleImgPtr, outSampleSegPtr, warpedImgPtrs, warpedSegPtrs, numOfAtlases); + + // output segmentation image + imgoperator->WriteImage(outSampleSegName, outSampleSegPtr, ioType); + + delete[] warpedImgPtrs; + delete[] warpedSegPtrs; + + return; +} + +void GaussianWeightedLabelFusion(InternalImageType::Pointer curSampleImgPtr, InternalImageType::Pointer outSampleSegPtr, + InternalImageType::Pointer* warpedImgPtrs, InternalImageType::Pointer* warpedSegPtrs, + unsigned numOfAtlases) +{ + // introduce the iterators of each image + // InternalImageIteratorType sampleImgIt(curSampleImgPtr, + // curSampleImgPtr->GetLargestPossibleRegion());//GetRequestedRegion()); + InternalImageIteratorType sampleSegIt(outSampleSegPtr, outSampleSegPtr->GetLargestPossibleRegion() ); // GetRequestedRegion()); + + // InternalImageIteratorType* warpedImgIts = new InternalImageIteratorType[numOfAtlases]; + InternalImageIteratorType* warpedSegIts = new InternalImageIteratorType[numOfAtlases]; + + typedef itk::ConstNeighborhoodIterator InternalImageNeighborhoodIteratorType; + InternalImageNeighborhoodIteratorType* warpedImgNeighborhoodIts = + new InternalImageNeighborhoodIteratorType[numOfAtlases]; + InternalImageType::SizeType radius; + radius[0] = localPatchSize; radius[1] = localPatchSize; radius[2] = localPatchSize; + + InternalImageNeighborhoodIteratorType sampleImgNeighborhoodIt(radius, curSampleImgPtr, + curSampleImgPtr->GetRequestedRegion() ); + + InternalImageNeighborhoodIteratorType sampleImgNeighborIt(radius, curSampleImgPtr, + curSampleImgPtr->GetLargestPossibleRegion() ); + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + // InternalImageIteratorType it1( warpedImgPtrs[i], + // warpedImgPtrs[i]->GetLargestPossibleRegion());//GetRequestedRegion() ); + // warpedImgIts[i] = it1; + + InternalImageIteratorType it2( warpedSegPtrs[i], warpedSegPtrs[i]->GetLargestPossibleRegion() ); // GetRequestedRegion()); + warpedSegIts[i] = it2; + + InternalImageNeighborhoodIteratorType neighborIt(radius, warpedImgPtrs[i], + warpedImgPtrs[i]->GetLargestPossibleRegion() ); + warpedImgNeighborhoodIts[i] = neighborIt; + } + unsigned maxLabel = 0; // the maximum labels + // find out the maximum label from the label images + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + warpedSegIts[i].GoToBegin(); + } + while( !warpedSegIts[0].IsAtEnd() ) + { + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + if( warpedSegIts[i].Get() > maxLabel ) + { + maxLabel = warpedSegIts[i].Get(); + } + ++warpedSegIts[i]; + } + } + + int temp1, temp2; + double kernel_sigma = 1; // std of Gaussian approximation, set to the median of all distances + int* index = new int[numOfAtlases]; + double* sort1 = new double[numOfAtlases]; + int* label_index = new int[maxLabel + 1]; + + int* label_pool = new int[numOfAtlases]; + double* mse = new double[numOfAtlases]; + double* weight_sum = new double[maxLabel + 1]; + // for each voxel do label fusion,(vx, vy, vz) is the current location + double kernel_sigma_square = kernel_sigma * kernel_sigma; + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + warpedImgNeighborhoodIts[i].GoToBegin(); + warpedSegIts[i].GoToBegin(); + } + for( sampleImgNeighborIt.GoToBegin(), sampleSegIt.GoToBegin(); + !sampleImgNeighborIt.IsAtEnd(); + ++sampleImgNeighborIt, ++sampleSegIt ) + { +// idx = sampleImgNeighborIt.GetIndex(); +// if (idx[0] == 30 && idx[1]==20 && idx[2] ==8) +// idx[0] = idx[0]; + double msec = 0.0; + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + label_pool[i] = warpedSegIts[i].Get(); + for( unsigned n = 0; n < sampleImgNeighborIt.Size(); ++n ) + { + temp1 = sampleImgNeighborIt.GetPixel(n); + temp2 = warpedImgNeighborhoodIts[i].GetPixel(n); + // add together differences, squared sum + msec += (temp1 - temp2) * (temp1 - temp2); + } + mse[i] = msec; + } + // sort the mean square difference + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + index[i] = i; + sort1[i] = mse[i]; + } + basicoperator->bubbleSort(sort1, index, numOfAtlases); + kernel_sigma_square = sort1[(int(numOfAtlases) - 1) / 2] + 0.0001; + // weight each label + for( unsigned i = 0; i < maxLabel + 1; ++i ) + { + weight_sum[i] = 0.0; + } + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + // add-up all weights + if ((mse[i]) / (2.0 * kernel_sigma_square) < 50) + weight_sum[label_pool[i]] += exp(0.0 - (mse[i]) / (2.0 * kernel_sigma_square) ); + + } + + // weighted label fusion + for( unsigned i = 0; i < maxLabel + 1; ++i ) + { + label_index[i] = i; + } + + // determine the label with the maximum weight + // basicoperator->bubbleSort(weight_sum, label_index, maxLabel+1); + // int label_final; + // label_final = label_index[maxLabel]; + unsigned label_final = 0; + double weight_final = 0.0; + for( unsigned i = 0; i < maxLabel + 1; ++i ) + { + if( weight_sum[i] > weight_final ) + { + label_final = i; + weight_final = weight_sum[i]; + } + } + + // assign label to the segmentations + sampleSegIt.Set(label_final); + for( unsigned i = 0; i < numOfAtlases; ++i ) + { + ++warpedImgNeighborhoodIts[i]; + ++warpedSegIts[i]; + } + } + /* + for (int vz = 0; vz < sampleImSize[2]; vz++) + { + //if (isDebug) + // std::cout << vz << ", "; + for (int vy = 0; vy < sampleImSize[1]; vy++) + { + for (int vx = 0; vx < sampleImSize[0]; vx++) + { + // process each voxel + idx.SetElement(0,vx); + idx.SetElement(1,vy); + idx.SetElement(2,vz); + + // get labels from all atlases on this current location + + for (int i = 0; i < numOfAtlases; ++i) + { + + warpedSegIts[i].SetIndex(idx); + label_pool[i] = warpedSegIts[i].Get(); + } + + // get the local patch differences + + for (int i = 0; i < numOfAtlases; ++i) mse[i] = 0.0; + for (int i = 0; i < numOfAtlases; ++i) + { + double msec = 0.0; + // compare with all other images at the local patch of current voxel + // (nx,ny,nz) is the incremental position on (vx,vy,vz) + for(int nz = 0-localPatchSize; nz <= localPatchSize; nz++) + { + for(int ny = 0-localPatchSize; ny <= localPatchSize; ny++) + { + for(int nx = 0-localPatchSize; nx <= localPatchSize; nx++) + { + if ((vx+nx>=0)&&(vx+nx=0)&&(vy+ny=0)&&(vz+nzbubbleSort(sort1, index, numOfAtlases); + kernel_sigma = sort1[(int)((numOfAtlases-1)/2)]+0.0001; + + // weight each label + + for (int i = 0; i < maxLabel+1; ++i) weight_sum[i] = 0.0; + for (int i = 0; i < numOfAtlases; ++i) + { + // add-up all weights + weight_sum[label_pool[i]]+= exp(0-(mse[i]*mse[i])/(2*kernel_sigma*kernel_sigma)); + } + + // weighted label fusion + for (int i = 0; i < maxLabel+1; ++i) label_index[i] = i; + + basicoperator->bubbleSort(weight_sum, label_index, maxLabel+1); + int label_final; + label_final = label_index[maxLabel]; + + // assign label to the segmentations + sampleSegIt.SetIndex(idx); + sampleSegIt.Set(label_final); + + } + } + }// end of for vz + +*/ + delete[] label_pool; + delete[] mse; + delete[] weight_sum; + // + // delete[] warpedImgIts; + delete[] warpedSegIts; + delete[] warpedImgNeighborhoodIts; + delete[] sort1; + delete[] index; + delete[] label_index; + + return; +} + +void TreeBasedRegistrationFastOniTree(vnl_vector itree, // the incremental tree + int root, // the tree root + int itree_size, // the tree size + itk::MABMISAtlas* atlasTree, // the atlas tree + itk::MABMISImageData*imageData, // the images to be segmented + std::vector iterations, + double sigma + ) +// (vnl_vector itree_t, int itree_size_t, int atlas_size_t, int simulate_size_t, int sample_size_t, std::vector +// originalIntensityImageFileNames, std::vector deformedImgFileNames, std::vector deformationFieldFileNames) +{ + // int root = -1; // + // treeoperator->FindRoot(itree_t, itree_size_t, root); + + // calculate the depth of each node + int* node_depth = new int[itree_size]; + + for( int i = 0; i < itree_size; ++i ) + { + int curnode = i; + node_depth[i] = 0; + bool isRoot = false; + + while( !isRoot ) + { + if( itree[curnode] >= 0 ) + { + node_depth[i]++; + curnode = itree[curnode]; + } + else + { + isRoot = true; + } + } + } + + // sort all nodes by its depth + int* index = new int[itree_size]; + { + // + double* arr = new double[itree_size]; + for( int i = 0; i < itree_size; ++i ) + { + arr[i] = node_depth[i]; + index[i] = i; + } + basicoperator->bubbleSort(arr, index, itree_size); // index[0] = root + delete[] arr; + } + + int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; + int atlas_total_size = atlasTree->m_NumberAllAtlases; + + // start to register each image to the root node step by step + for( int ii = 1; ii < itree_size; ++ii ) // starting from 1, since index[0] showing the root + { + // if (isDebug) + // std::cout << ii << ", "; + int curnode = index[ii]; + int parentnode = itree[curnode]; + + std::cout << "Between nodes: [" << curnode << ", " << parentnode << "]" << std::endl; + + if( (curnode >= atlas_image_size) && (curnode < atlas_total_size) ) + { + // this is a simulated image, then pass + continue; + } + + // first check whether registrations from atlas images to root exist. If exist, no registration is needed + bool isRegistered = false; + char i_str[10], root_str[10], p_str[10]; + sprintf(i_str, "%03d", curnode); + sprintf(root_str, "%03d", root); + sprintf(p_str, "%03d", parentnode); + + std::string dfFileName, df_ImageFileName; + if( curnode == root ) + { + isRegistered = true; + } + else if( curnode < atlas_image_size ) + { + dfFileName = std::string(i_str) + "_to_" + root_str + DeformationSuffix; + df_ImageFileName = std::string(i_str) + "_to_" + root_str + RegisteredSuffix; + dfFileName = atlasTree->m_AtlasDirectory + dfFileName; + df_ImageFileName = atlasTree->m_AtlasDirectory + df_ImageFileName; + + if( (itksys::SystemTools::FileExists(dfFileName.c_str(), true) && + itksys::SystemTools::FileExists(df_ImageFileName.c_str(), true) ) + ) + { + isRegistered = true; + } + } + else if( curnode >= atlas_total_size ) + { + int n = curnode - atlas_total_size; + char n_str[10]; sprintf(n_str, "%03d", n); + dfFileName = std::string("testImage") + n_str + "_to_" + root_str + DeformationSuffix; + df_ImageFileName = std::string("testImage") + n_str + "_to_" + root_str + RegisteredSuffix; + + dfFileName = imageData->m_OutputDirectory + dfFileName; + df_ImageFileName = imageData->m_OutputDirectory + df_ImageFileName; + + isRegistered = false; + } + // parent deformation fields + std::string df_ParentFileName; + if( parentnode < atlas_image_size ) + { + df_ParentFileName = std::string(p_str) + "_to_" + root_str + DeformationSuffix; + df_ParentFileName = atlasTree->m_AtlasDirectory + df_ParentFileName; + } + else if( parentnode >= atlas_total_size ) + { + int n = parentnode - atlas_total_size; + char n_str[10]; sprintf(n_str, "%03d", n); + df_ParentFileName = std::string("testImage") + n_str + "_to_" + root_str + DeformationSuffix; + df_ParentFileName = imageData->m_OutputDirectory + df_ParentFileName; + } + + // do registration + if( isRegistered ) + { + continue; // goto next one + } + + std::string rootImageFile = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[root]); + std::string testImageFile; + if( curnode < atlas_image_size ) + { + testImageFile = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[curnode]); // atlas image + } + else + { + testImageFile = imageData->m_DataDirectory + imageData->m_ImageFileNames[curnode - atlas_total_size]; // test image + } + if( parentnode == root ) // direct registration without initial deformation + { + regoperator->DiffeoDemonsRegistrationWithParameters( + rootImageFile, + testImageFile, + df_ImageFileName, + dfFileName, + sigma, doHistMatch, iterations); + } + else // registration with initial deformation + { + if( (parentnode < atlas_image_size) || (parentnode >= atlas_total_size) ) + { + // parent is a real image + regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( + rootImageFile, + testImageFile, + df_ParentFileName, + df_ImageFileName, + dfFileName, + sigma, doHistMatch, iterations); + } + else + { + // parent is a simulated image + + std::string index_string = basicoperator->myitoa( parentnode - atlas_image_size, 3 ); + std::string sim_df_FileName = "simulated_deform_" + index_string + ".nii.gz"; + sim_df_FileName = atlasTree->m_AtlasDirectory + sim_df_FileName; + + regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( + rootImageFile, + testImageFile, + sim_df_FileName, + df_ImageFileName, + dfFileName, + sigma, doHistMatch, iterations); + } + } + } + std::cout << "Done." << std::endl; + + delete[] node_depth; + delete[] index; +} + +// other filter +// histogram matching +void HistogramMatching(InternalImageType::Pointer inputImage, + InternalImageType::Pointer referenceImage, + InternalImageType::Pointer & outputImage) +{ + InternalHistMatchFilterType::Pointer matcher = InternalHistMatchFilterType::New(); + + matcher->SetInput( inputImage ); + matcher->SetReferenceImage( referenceImage ); + matcher->SetNumberOfHistogramLevels( 1024 ); // 1024 + matcher->SetNumberOfMatchPoints( 7 ); + matcher->ThresholdAtMeanIntensityOn(); + try + { + matcher->Update(); + } + catch( itk::ExceptionObject & err ) + { + std::cout << err << std::endl; + return; + } + outputImage = matcher->GetOutput(); + return; +} + +// void RegistrationOntoTreeRoot(vnl_vector itree, int root, int filenumber, int itree_size, int atlas_size, int +// simulate_size, int sample_size, std::vector sub_ids, std::vector segfilenames) +void RegistrationOntoTreeRoot(vnl_vector itree, // the incremental tree + int root, // the tree root + int itree_size, // the tree size + itk::MABMISAtlas* atlasTree, // the atlas tree + itk::MABMISImageData*imageData, // the images to be segmented + std::vector iterations, // registration parameters + double sigma + ) +{ + // do registration on the tree + std::cout << "Start registration ... " << std::endl; + + // to get the numbers original defined in the data; + const int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; + const int test_image_size = imageData->m_NumberImageData; + +// do tree-based registration + TreeBasedRegistrationFastOniTree(itree, // the incremental tree + root, // the tree root + itree_size, // the tree size + atlasTree, // the atlas tree + imageData, // the images to be segmented + iterations, sigma // registration parameters + ); + + // end of do registration to the root + + // + // basicoperator->RemoveFile("simulated_*"); + + std::cout << "Reverse deformation field ... " << std::endl;; + // reverse each deformation field (from the root) + std::string rootImageTag; + char root_str[10], i_str[10]; + sprintf(root_str, "%03d", root); + rootImageTag = root_str; + std::string fixedImageTag; + + for( int i = 0; i < atlas_image_size + test_image_size; i++ ) + { + //if( isDebug ) + { + std::cout << i << ", "; + } + if( i == root ) + { + continue; + } + + // prepare file names + + // std::string originalImgImageFileName = sub_ids[root] + ImageSuffix; + // std::string originalSegImageFileName = sub_ids[root] + SegmentationSuffix; + std::string rootImageFileName = ReplacePathSepForOS( + atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[root] ); + std::string rootSegmentFileName = ReplacePathSepForOS( + atlasTree->m_AtlasDirectory + atlasTree->m_AtlasSegmentationFilenames[root] ); + + std::string fixedImageFileName; + std::string deformedImageFileName; + std::string deformedImageFileNameImg; + std::string deformedSegmentFileName; + + std::string deformationFileName; + std::string invDeformationFileName; + + if( i < atlas_image_size ) + { + fixedImageFileName = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[i]); + sprintf(i_str, "%03d", i); + fixedImageTag = i_str; + + deformedImageFileName = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + ImageSuffix; + deformedImageFileNameImg = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + ImageSuffix; + + deformedSegmentFileName = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + SegmentationSuffix; + + deformationFileName = atlasTree->m_AtlasDirectory + fixedImageTag + "_to_" + rootImageTag + DeformationSuffix; + invDeformationFileName = atlasTree->m_AtlasDirectory + rootImageTag + "_to_" + fixedImageTag + DeformationSuffix; + } + else + { + fixedImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[i - atlas_image_size]; + sprintf(i_str, "%03d", i - atlas_image_size); + fixedImageTag = std::string("testImage") + i_str; + + deformedImageFileName = imageData->m_OutputDirectory + rootImageTag + "_to_" + fixedImageTag + ImageSuffix; + deformedImageFileNameImg = imageData->m_OutputDirectory + rootImageTag + "_to_" + fixedImageTag + ImageSuffix; + + deformedSegmentFileName = imageData->m_OutputDirectory + rootImageTag + "_to_" + fixedImageTag + SegmentationSuffix; + + deformationFileName = imageData->m_OutputDirectory + fixedImageTag + "_to_" + rootImageTag + DeformationSuffix; + invDeformationFileName = imageData->m_OutputDirectory + rootImageTag + "_to_" + fixedImageTag + DeformationSuffix; + } + + // if exist?? + bool isReversedExist = true; + if( !itksys::SystemTools::FileExists(deformedImageFileName.c_str(), true) ) + { + isReversedExist = false; + } + if( !itksys::SystemTools::FileExists(deformedImageFileNameImg.c_str(), true) ) + { + isReversedExist = false; + } + if( !itksys::SystemTools::FileExists(invDeformationFileName.c_str(), true) ) + { + isReversedExist = false; + } + if( isReversedExist ) + { + continue; + } + + // reverse deformation field + DeformationFieldType::Pointer deformationField = nullptr; + dfoperator->ReadDeformationField(deformationFileName, deformationField); + + DeformationFieldType::Pointer inversedDeformationField = DeformationFieldType::New(); + inversedDeformationField->SetRegions(deformationField->GetLargestPossibleRegion() ); + inversedDeformationField->SetSpacing(deformationField->GetSpacing() ); + inversedDeformationField->SetDirection(deformationField->GetDirection() ); + inversedDeformationField->SetOrigin(deformationField->GetOrigin() ); + inversedDeformationField->Allocate(); + + dfoperator->InverseDeformationField3D(deformationField, inversedDeformationField); + // output reversed deformation field + + dfoperator->WriteDeformationField(invDeformationFileName, inversedDeformationField); + // //update + regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( + fixedImageFileName, rootImageFileName, + invDeformationFileName, + deformedImageFileName, invDeformationFileName, + sigma, doHistMatch, iterations); + + // CompressDeformationField2Short(inversedDeformationFieldFileName); + // apply deformation field on seg file + if( isEvaluate ) + { + dfoperator->ApplyDeformationFieldAndWriteWithTypeWithFileNames( + rootSegmentFileName, + invDeformationFileName, + deformedSegmentFileName, false); + } + + std::cout << "Done!" << std::endl; + } // end of reverse each deformation field (from the root) + // std::cout << "done." << std::endl; + std::cout << std::endl; +} + +void PairwiseRegistrationOnTreeViaRoot(int root, + itk::MABMISImageData* imageData, + itk::MABMISAtlas* atlasTree, + std::vector iterations, + double sigma) +{ + ///////////////////////////////// + // with the help of root node, obtain the pairwise registration among all images + std::cout << "Generate all pairwise registration ..." << std::endl;; + + int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; + int test_image_size = imageData->m_NumberImageData; + + // do for all real images, including both atlases and test images + for( int all_index = 0; all_index < atlas_image_size + test_image_size; all_index++ ) + { + if( isDebug ) + { + std::cout << all_index << ": "; + } + if( all_index == root ) + { + continue; + } + + char a_str[10], s_str[10], root_str[10]; + sprintf(root_str, "%03d", root); + sprintf(a_str, "%03d", all_index); + + std::string movingImageFileName; + std::string movingSegmentFileName; + // get the moving filenames + std::string movingImageTag; + if( all_index < atlas_image_size ) + { + movingImageFileName = ReplacePathSepForOS(atlasTree->m_AtlasDirectory + atlasTree->m_AtlasFilenames[all_index]); + movingSegmentFileName = atlasTree->m_AtlasDirectory + atlasTree->m_AtlasSegmentationFilenames[all_index]; + movingImageTag = std::string(a_str); + } + else + { + movingImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[all_index - atlas_image_size]; + movingSegmentFileName = movingImageFileName; + const std::string movingSegmentBaseFileName = GetRootName(movingSegmentFileName); + + movingSegmentFileName = movingSegmentBaseFileName + SegmentationSuffix; + sprintf(a_str, "%03d", all_index - atlas_image_size); + movingImageTag = std::string("testImage") + a_str; + } + // for (int sample_index = atlas_size; sample_index < filenumber; sample_index++) + for( int sample_index = 0; sample_index < int(imageData->m_NumberImageData); sample_index++ ) + { + if( isDebug ) + { + std::cout << "[" << all_index << "," << sample_index << "], "; + } + + if( sample_index == all_index - atlas_image_size ) + { + continue; + } + + std::string fixedImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[sample_index]; + + sprintf(s_str, "%03d", sample_index); + std::string fixedImageTag = std::string("testImage") + s_str; + + std::string deformedImageFileName = movingImageTag + "_to_" + fixedImageTag + ImageSuffix; + std::string deformedImageFileNameImg = movingImageTag + "_to_" + fixedImageTag + ImageSuffix; + std::string deformedSegmentFileName = movingImageTag + "_to_" + "testImage" + s_str + SegmentationSuffix; + std::string deformedSegmentFileNameImg = movingImageTag + "_to_" + "testImage" + s_str + SegmentationSuffix; + std::string dfFileName = movingImageTag + "_to_" + fixedImageTag + DeformationSuffix; + + deformedImageFileName = imageData->m_OutputDirectory + deformedImageFileName; + deformedImageFileNameImg = imageData->m_OutputDirectory + deformedImageFileNameImg; + dfFileName = imageData->m_OutputDirectory + dfFileName; + + deformedSegmentFileName = imageData->m_OutputDirectory + deformedSegmentFileName; + deformedSegmentFileNameImg = imageData->m_OutputDirectory + deformedSegmentFileNameImg; + + // std::string fixedImgImageFileName = sub_ids[sample_index] + ImageSuffix; + // std::string movingImgImageFileName = sub_ids[all_index] + ImageSuffix; + // std::string movingSegImageFileName = sub_ids[all_index] + SegmentationSuffix; + + // std::string deformedImgImageFileName = sub_ids[sample_index] + "_" + sub_ids[all_index] + ImageSuffix; + // std::string deformedSegImageFileName = sub_ids[sample_index] + "_" + sub_ids[all_index] + SegmentationSuffix; + // std::string deformedImgImageFileNameImg = sub_ids[sample_index] + "_" + sub_ids[all_index] + ImageSuffix; + // std::string deformedSegImageFileNameImg = sub_ids[sample_index] + "_" + sub_ids[all_index] + SegmentationSuffix; + // std::string deformationFieldFileName = sub_ids[sample_index] + "_" + sub_ids[all_index] + DeformationSuffix; + + // if exist + bool isComposeExist = true; + if( !itksys::SystemTools::FileExists(deformedImageFileName.c_str(), true) ) + { + isComposeExist = false; + } + if( !itksys::SystemTools::FileExists(deformedImageFileNameImg.c_str(), true) ) + { + isComposeExist = false; + } + + if( !itksys::SystemTools::FileExists(dfFileName.c_str(), true) ) + { + isComposeExist = false; + } + if( isComposeExist ) + { + continue; + } + + // compose + // std::string inputDeformationFieldFileName = sub_ids[root] + "_" + sub_ids[all_index] + DeformationSuffix; + // std::string secondDeformationFieldFileName = sub_ids[sample_index] + "_" + sub_ids[root] + DeformationSuffix; + std::string inputDeformationFieldFileName; + if( all_index < atlas_image_size ) + { + inputDeformationFieldFileName = atlasTree->m_AtlasDirectory + movingImageTag + "_to_" + root_str + DeformationSuffix; + } + else + { + inputDeformationFieldFileName = imageData->m_OutputDirectory + movingImageTag + "_to_" + root_str + DeformationSuffix; + } + std::string secondDeformationFieldFileName = imageData->m_OutputDirectory + root_str + "_to_" + fixedImageTag + DeformationSuffix; + // output composed deformation field + dfoperator->ComposeDeformationFieldsAndSave(inputDeformationFieldFileName, + secondDeformationFieldFileName, + dfFileName); + + // update: refine the deformation + regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( + fixedImageFileName, + movingImageFileName, + dfFileName, + deformedImageFileName, + dfFileName, + sigma, doHistMatch, iterations); + + if( isEvaluate ) + { + dfoperator->ApplyDeformationFieldAndWriteWithTypeWithFileNames( + movingSegmentFileName, + dfFileName, + deformedSegmentFileName, false); + } + } // end of sample_index + } // end of all_index + + std::cout << std::endl; + std::cout << "Done!" << std::endl; +} + +int MultiAtlasBasedSegmentation(itk::MABMISImageData* imageData, + itk::MABMISAtlas* atlasTree) +{ + std::cout << "Start to process segmentation..." << std::endl; + int iter = 1; + + bool isConverged = false; // the iterative update can stop + + int atlas_image_size = atlasTree->m_NumberAllAtlases - atlasTree->m_NumberSimulatedAtlases; + + int test_image_size = imageData->m_NumberImageData; + + while( !isConverged ) + { + // if (isDebug) + std::cout << "iteration " << iter << ": "; + // iter == 1, only using atlases; otherwise using all images + // number of other image to be used + int numOtherAtlases = 0; + if( iter == 1 ) + { + numOtherAtlases = atlas_image_size; + } + else + { + numOtherAtlases = atlas_image_size + test_image_size - 1; + } + // do label fusion for each new sample + for( int sample_index = 0; sample_index < test_image_size; sample_index++ ) + { + // if (isDebug) + // std::cout << sample_index << ":" << sub_ids[sample_index] << ", "; + std::cout << sample_index << ", "; + + char s_str[10]; + sprintf(s_str, "%03d", sample_index); + std::string testImageTag = std::string("testImage") + s_str; + + // the sample to be processed + + std::string testImageFileName = imageData->m_DataDirectory + imageData->m_ImageFileNames[sample_index]; + + // std::string curSampleImgName = sub_ids[sample_index] +ImageSuffix; + // the output label image of current sample + std::string index_string = basicoperator->myitoa( iter, 3 ); + + const std::string testSegmentFileName = GetDefaultSegmentationFilename(imageData, sample_index, index_string); + + // prepare all current atlases names + std::vector allWarpedAtlasImgNames; + std::vector allDeformationFieldNames; + std::vector allAtlasSegNames; + + int atlas_index = 0; + for( int i = 0; i < atlas_image_size + test_image_size; ++i ) + { + if( (iter == 1) && (i >= atlas_image_size) ) + { + continue; + } + if( i - atlas_image_size == sample_index ) + { + continue; + } + + std::string movingImageTag; + char i_str[10]; + + if( i < atlas_image_size ) + { + sprintf(i_str, "%03d", i); + movingImageTag = std::string(i_str); + } + else + { + sprintf(i_str, "%03d", i - atlas_image_size); + movingImageTag = std::string("testImage") + i_str; + } + + // assign names + std::string warpedImageFileName = movingImageTag + "_to_" + testImageTag + ImageSuffix; + std::string dfFileName = movingImageTag + "_to_" + testImageTag + DeformationSuffix; + + warpedImageFileName = imageData->m_OutputDirectory + warpedImageFileName; + dfFileName = imageData->m_OutputDirectory + dfFileName; + + std::string atlasSegName; + if( i < atlas_image_size ) + { + atlasSegName = atlasTree->m_AtlasDirectory + atlasTree->m_AtlasSegmentationFilenames[i]; + } + else + { + index_string = basicoperator->myitoa( iter-1, 3 ); + atlasSegName = GetDefaultSegmentationFilename(imageData, i - atlas_image_size, index_string); + } + + allWarpedAtlasImgNames.push_back(warpedImageFileName); + allDeformationFieldNames.push_back(dfFileName); + allAtlasSegNames.push_back(atlasSegName); + + // go to next atlas + atlas_index++; + } + // do weighted label fusion + + LabelFusion(testImageFileName, + testSegmentFileName, + allWarpedAtlasImgNames, + allDeformationFieldNames, + allAtlasSegNames, + numOtherAtlases); + } + std::cout << std::endl; + + // to decide if converged + if( iter >= 3 ) + { + isConverged = true; // + } + // calculate overlap rate + if( isEvaluate ) + { + // pairwise segmentation accuracy + // pairwise segmentation consistency + } + iter++; + } + + return iter; +} diff --git a/Testing.h b/Testing.h new file mode 100644 index 0000000..0c417eb --- /dev/null +++ b/Testing.h @@ -0,0 +1,9 @@ +#ifndef __Testing_h +#define __Testing_h + +#include "commonMABMIS.h" + +int Testing(itk::MABMISImageData* imageData, itk::MABMISAtlas* atlasTree, + std::vector iterations, double sigma); + +#endif //__Testing_h diff --git a/Testing/CMakeLists.txt b/Testing/CMakeLists.txt index 078322b..fddde79 100644 --- a/Testing/CMakeLists.txt +++ b/Testing/CMakeLists.txt @@ -3,80 +3,95 @@ set(TEMP ${CMAKE_CURRENT_BINARY_DIR}/../Testing/Temporary) set(INPUT ${CMAKE_CURRENT_SOURCE_DIR}/../Data/Input) set(BASELINE ${CMAKE_CURRENT_SOURCE_DIR}/../Data/Baseline) -message (STATUS ${CMAKE_CURRENT_BINARY_DIR}) -message (STATUS ${CMAKE_CURRENT_SOURCE_DIR}) - - -macro(configure_files srcDir destDir) - message(STATUS "Configuring directory ${destDir}") - make_directory(${destDir}) - - file(GLOB templateFiles RELATIVE ${srcDir} ${srcDir}/*) - foreach(templateFile ${templateFiles}) - set(srcTemplatePath ${srcDir}/${templateFile}) - if(NOT IS_DIRECTORY ${srcTemplatePath}) - message(STATUS "Configuring file ${templateFile}") - configure_file( - ${srcTemplatePath} - ${destDir}/${templateFile} - COPYONLY) - endif(NOT IS_DIRECTORY ${srcTemplatePath}) - endforeach(templateFile) -endmacro(configure_files) - #----------------------------------------------------------------------------- # Training module set(Training_CLP IGR3D_MABMIS_Training) #----------------------------------------------------------------------------- -add_executable(${Training_CLP}Test ${Training_CLP}Test.cxx) -target_link_libraries(${Training_CLP}Test ${Training_CLP}Lib) -set_target_properties(${Training_CLP}Test PROPERTIES LABELS ${Training_CLP}) - - -configure_file(${INPUT}/TrainingData/TrainingData.xml ${TEMP}/TrainingData.xml) -configure_files(${INPUT}/TrainingData/ ${TEMP}/) +add_executable(${Training_CLP}Test Test.cxx) +target_link_libraries(${Training_CLP}Test ${Training_CLP}Lib ${SlicerExecutionModel_EXTRA_EXECUTABLE_TARGET_LIBRARIES}) #----------------------------------------------------------------------------- set(testname ${Training_CLP}Test) ExternalData_add_test(${Training_CLP}Data NAME ${testname} COMMAND ${Slicer_LAUNCH_COMMAND} $ ModuleEntryPoint -i 5,3,2 -s 1.5 - --trainingXML ${TEMP}/TrainingData.xml --atlasTreeXML ${TEMP}/TrainedAtlas.xml + --trainingXML ${INPUT}/TrainingData/TrainingData.xml --atlasTreeXML ${TEMP}/TrainedAtlas.xml ) -set_property(TEST ${testname} PROPERTY LABELS ${Training_CLP}) +set_property(TEST ${testname} PROPERTY LABELS MABMIS) +set(testname ${Training_CLP}_compareTraingingXMLs) if (WIN32) - set(testname compareTraingingXMLs) ExternalData_add_test(${Training_CLP}Data NAME ${testname} COMMAND diff ${INPUT}/AtlasTree/TrainedAtlas-win.xml ${TEMP}/TrainedAtlas.xml) else() - set(testname compareTraingingXMLs) ExternalData_add_test(${Training_CLP}Data NAME ${testname} COMMAND diff ${INPUT}/AtlasTree/TrainedAtlas.xml ${TEMP}/TrainedAtlas.xml) endif() -set_property(TEST compareTraingingXMLs APPEND PROPERTY DEPENDS ${Training_CLP}Test) +set_property(TEST ${testname} PROPERTY LABELS MABMIS) +set_property(TEST ${testname} APPEND PROPERTY DEPENDS ${Training_CLP}Test) #----------------------------------------------------------------------------- # Testing module set(Testing_CLP IGR3D_MABMIS_Testing) #----------------------------------------------------------------------------- -add_executable(${Testing_CLP}Test ${Testing_CLP}Test.cxx) -target_link_libraries(${Testing_CLP}Test ${Testing_CLP}Lib) -set_target_properties(${Testing_CLP}Test PROPERTIES LABELS ${Testing_CLP}) +add_executable(${Testing_CLP}Test Test.cxx) +target_link_libraries(${Testing_CLP}Test ${Testing_CLP}Lib ${SlicerExecutionModel_EXTRA_EXECUTABLE_TARGET_LIBRARIES}) #----------------------------------------------------------------------------- -configure_file(${INPUT}/AtlasTree/TrainedAtlas.xml ${TEMP}/TrainedAtlas.xml) -configure_files(${INPUT}/AtlasTree/TrainedAtlas ${TEMP}/TrainedAtlas/) - -configure_files(${INPUT}/TestData ${TEMP}/TestData/) - set(testname ${Testing_CLP}Test) ExternalData_add_test(${Testing_CLP}Data NAME ${testname} COMMAND ${Slicer_LAUNCH_COMMAND} $ --compareRadiusTolerance 3 --compare DATA{${BASELINE}/UA_cbq_000_seg.nii.gz} - ${TEMP}/UA_cbq_000_seg.nii.gz + ${TEMP}/TestResults/UA_cbq_000_seg.nii.gz + ModuleEntryPoint + -i 5,3,2 -s 1.5 -o ${TEMP}/TestResults + --atlasTreeXML ${INPUT}/AtlasTree/TrainedAtlas.xml --imageListXML ${INPUT}/TestData/TestData.xml + ) +set_property(TEST ${testname} PROPERTY LABELS MABMIS) + + +#----------------------------------------------------------------------------- +# Direct invoke module +set(Testing_CLP IGR3D_MABMIS_Direct_Invoke) + +add_executable(${Testing_CLP}Test Test.cxx) +target_link_libraries(${Testing_CLP}Test ${Testing_CLP}Lib ${SlicerExecutionModel_EXTRA_EXECUTABLE_TARGET_LIBRARIES}) + +# Test creating the atlas +set(testname ${Testing_CLP}TrainingTest) +ExternalData_add_test(${Testing_CLP}Data NAME ${testname} COMMAND ${Slicer_LAUNCH_COMMAND} $ + ModuleEntryPoint + -i 5,3,2 -s 1.5 --mode "Train atlas" --imageDir ${TEMP}/TrainingTemp + --atlasTreeXML ${TEMP}/Atlas.xml + --image0 ${INPUT}/TrainingData/AA_cbq_000.nii.gz --segmentation0 ${INPUT}/TrainingData/AA_seg_000.nii.gz + --image1 ${INPUT}/TrainingData/AB_cbq_000.nii.gz --segmentation1 ${INPUT}/TrainingData/AB_seg_000.nii.gz + --image2 ${INPUT}/TrainingData/AC_cbq_000.nii.gz --segmentation2 ${INPUT}/TrainingData/AC_seg_000.nii.gz + --image3 ${INPUT}/TrainingData/AD_cbq_000.nii.gz --segmentation3 ${INPUT}/TrainingData/AD_seg_000.nii.gz + ) +set_property(TEST ${testname} PROPERTY LABELS MABMIS) + +# Test using the atlas +set(testname ${Testing_CLP}SegmentationTest) +ExternalData_add_test(${Testing_CLP}Data NAME ${testname} COMMAND ${Slicer_LAUNCH_COMMAND} $ + --compareRadiusTolerance 3 --compare DATA{${BASELINE}/UA_cbq_000_seg.nii.gz} ${TEMP}/SegmentationResults/UA-label.nrrd + ModuleEntryPoint + -i 5,3,2 -s 1.5 --mode "Direct invoke" --imageDir ${TEMP}/SegTemp --segmentationInterpolation Nearest + --atlasTreeXML ${TEMP}/Atlas.xml --imageListXML ${TEMP}/List.xml + --image0 ${INPUT}/TestData/UA_cbq_000.nii.gz --segmentation0 ${TEMP}/SegmentationResults/UA-label.nrrd + --image1 ${INPUT}/TestData/UB_cbq_000.nii.gz --segmentation1 ${TEMP}/SegmentationResults/UB-label.nrrd + --image2 ${INPUT}/TestData/VA_cbq_000.nii.gz --segmentation2 ${TEMP}/SegmentationResults/VA-label.nrrd + --image3 ${INPUT}/TestData/VB_cbq_000.nii.gz --segmentation3 ${TEMP}/SegmentationResults/VB-label.nrrd + ) +set_property(TEST ${testname} PROPERTY LABELS MABMIS) +set_property(TEST ${testname} APPEND PROPERTY DEPENDS ${Testing_CLP}TrainingTest) + +set(testname ${Testing_CLP}SegmentationTest2) +ExternalData_add_test(${Testing_CLP}Data NAME ${testname} COMMAND ${Slicer_LAUNCH_COMMAND} $ + --compareRadiusTolerance 4 --compare DATA{${BASELINE}/UA_cbq_000_seg.nii.gz} ${TEMP}/Segmentation2Result/UA-label.nrrd ModuleEntryPoint - -i 5,3,2 -s 1.5 -o ${TEMP} - --atlasTreeXML ${TEMP}/TrainedAtlas.xml --imageListXML ${TEMP}/TestData/TestData.xml + -i 5,3,2 -s 1.5 --mode "Direct invoke" --imageDir ${TEMP}/Seg2Temp --segmentationInterpolation Nearest + --atlasTreeXML ${TEMP}/Atlas.xml --imageListXML ${TEMP}/List.xml #just one image and higher tolerance + --image0 ${INPUT}/TestData/UA_cbq_000.nii.gz --segmentation0 ${TEMP}/Segmentation2Result/UA-label.nrrd ) -set_property(TEST ${testname} PROPERTY LABELS ${Testing_CLP}) +set_property(TEST ${testname} PROPERTY LABELS MABMIS) +set_property(TEST ${testname} APPEND PROPERTY DEPENDS ${Testing_CLP}TrainingTest) diff --git a/Testing/IGR3D_MABMIS_TrainingTest.cxx b/Testing/IGR3D_MABMIS_TrainingTest.cxx deleted file mode 100644 index d2f1e5f..0000000 --- a/Testing/IGR3D_MABMIS_TrainingTest.cxx +++ /dev/null @@ -1,14 +0,0 @@ -#include "itkTestMain.h" - -#ifdef WIN32 -#define MODULE_IMPORT __declspec(dllimport) -#else -#define MODULE_IMPORT -#endif - -extern "C" MODULE_IMPORT int ModuleEntryPoint(int, char * []); - -void RegisterTests() -{ - StringToTestFunctionMap["ModuleEntryPoint"] = ModuleEntryPoint; -} diff --git a/Testing/IGR3D_MABMIS_TestingTest.cxx b/Testing/Test.cxx similarity index 100% rename from Testing/IGR3D_MABMIS_TestingTest.cxx rename to Testing/Test.cxx diff --git a/Training.cxx b/Training.cxx new file mode 100644 index 0000000..575b960 --- /dev/null +++ b/Training.cxx @@ -0,0 +1,571 @@ +/*========================================================================= + Copyright (c) IDEA LAB, UNC-Chapel Hill, 2013. + + MABMIS (Multi-Atlas-Based Multi-Image Segmentation) + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#include "Training.h" + +DeformationFieldType::SpacingType df_spacing; +DeformationFieldType::DirectionType df_direction; +DeformationFieldType::PointType df_origin; + +typedef itk::Vector ShortVectorPixelType; +typedef itk::Image ShortDeformationFieldType; +typedef itk::ImageFileWriter ShortDeformationFieldWriterType; + +void SearchRootAmongAtlases(std::vector imgfilenames, int & root); + +int RegistrationBetweenRootandAtlases(int root, std::vector imageFileNames, std::vector iterations, double sigma); + +int BuildStatisticalDeformationModel(int root, std::vector imageFileNames, int simulatedAtlasSize); + +std::vector GenerateSimulatedData(int root, std::vector imgfilenames, int simulatedAtlasSize); + +void strtrim(std::string& str) +{ + std::string::size_type pos = str.find_last_not_of(' '); + + if( pos != std::string::npos ) + { + str.erase(pos + 1); + pos = str.find_first_not_of(' '); + if( pos != std::string::npos ) + { + str.erase(0, pos); + } + } + else + { + str.erase(str.begin(), str.end() ); + } +} + +int Training( itk::MABMISImageData* trainingData, std::string outputFile, + std::vector iterations, double sigma) +{ + // sanity check + std::cout << "m_DataDirectory: " << trainingData->m_DataDirectory << std::endl; + if( trainingData->m_SegmentationFileNames.size() != trainingData->m_ImageFileNames.size() ) + { + std::cerr << "The numbers of image files and segmentation files are NOT equal!!!" << std::endl; + return -1; + } + if( trainingData->m_ImageFileNames.size() == 0 ) + { + std::cerr << "No Atlas is provided. QUIT. " << std::endl; + return -1; + } + if( trainingData->m_NumberImageData != trainingData->m_ImageFileNames.size() ) + { + trainingData->m_NumberImageData = trainingData->m_ImageFileNames.size(); + } + + // step 0: files and folders + // get output path from the output file name + size_t sep = outputFile.find_last_of(FILESEP); + std::string outputFolder = ""; + if( sep != std::string::npos ) + { + outputFolder = outputFile.substr(0, sep) + FILESEP; + } + else + { + sep = -1; + } + // get output file name without extension, and use it to create a directory to save trained atlas + size_t fnSep = outputFile.find_last_of('.'); + std::string outputxmlname = ""; + if( fnSep != std::string::npos ) + { + outputxmlname = outputFile.substr(sep + 1, fnSep - sep - 1); + } + else + { + outputxmlname = outputFile.substr(sep + 1, std::string::npos); + } + + std::string atlasFolder; + if( outputFolder.empty() ) + { + atlasFolder = outputxmlname; + } + else + { + atlasFolder = outputFolder + outputxmlname; + } + atlasFolder = ReplacePathSepForOS(atlasFolder); + if( !atlasFolder.empty() ) + { + atlasFolder = atlasFolder + FILESEP; + } + itksys::SystemTools::MakeDirectory(atlasFolder.c_str() ); + + // copy training data into the atlas folder + std::vector imageFiles(trainingData->m_NumberImageData); + std::vector segmentFiles(trainingData->m_NumberImageData); + for( unsigned i = 0; i < imageFiles.size(); ++i ) + { + imageFiles[i] = ReplacePathSepForOS(atlasFolder + trainingData->m_ImageFileNames[i]); + segmentFiles[i] = ReplacePathSepForOS(atlasFolder + trainingData->m_SegmentationFileNames[i]); + + // copy + std::string fromImagefile = trainingData->m_ImageFileNames[i]; + if( !trainingData->m_DataDirectory.empty() ) + { + fromImagefile = ReplacePathSepForOS(trainingData->m_DataDirectory + trainingData->m_ImageFileNames[i]); + } + + bool ret1=itksys::SystemTools::CopyFileAlways(fromImagefile.c_str(), imageFiles[i].c_str()); + if (!ret1) + { + std::cerr << "ERROR: Cannot copy atlas image file to trained atlas folder!" << std::endl; + return -1; + } + + // if file extension is '.nii.gz', also copy the 'img file' + if( itksys::SystemTools::Strucmp(fromImagefile.substr(fromImagefile.size() - 3, 3).c_str(), "hdr") == 0 ) + { + std::string fromfileImg = fromImagefile.substr(0, fromImagefile.size() - 3) + "img"; + std::string tofileImg = imageFiles[i].substr(0, imageFiles[i].size() - 3) + "img"; + itksys::SystemTools::CopyFileAlways(fromfileImg.c_str(), tofileImg.c_str()); + } + + std::string fromSegfile = trainingData->m_DataDirectory + trainingData->m_SegmentationFileNames[i]; + + bool ret2 = itksys::SystemTools::CopyFileAlways(fromSegfile.c_str(), segmentFiles[i].c_str()); + if (!ret2) + { + std::cerr << "ERROR: Cannot copy atlas segmentation file to trained atlas folder!" << std::endl; + return -1; + } + // if file extension is '.nii.gz', also copy the 'img file' + if( itksys::SystemTools::Strucmp(fromSegfile.substr(fromSegfile.size() - 3, 3).c_str(), "hdr") == 0 ) + { + std::string fromfileImg = fromSegfile.substr(0, fromSegfile.size() - 3) + "img"; + std::string tofileImg = segmentFiles[i].substr(0, segmentFiles[i].size() - 3) + "img"; + itksys::SystemTools::CopyFileAlways(fromfileImg.c_str(), tofileImg.c_str()); + } + } + + // step 1: find root of the atlas + std::cout << "---------------------------------------" << std::endl; + std::cout << "1. Find root atlas ... " << std::endl;; + //////////////////////////////// + // build tree on atlases to find the root + int root = -1; // root node of tree + SearchRootAmongAtlases(imageFiles, root); + std::cout << "Done." << std::endl;; + + // step 2.1: register all atlas to the root + ////////////////////////////// + // registration between root and other atlases + std::cout << "---------------------------------------" << std::endl; + std::cout << "2. Run coarse registration ... " << std::endl; + int val = RegistrationBetweenRootandAtlases(root, imageFiles, iterations, sigma); + if( val != 0 ) + { + return -1; + } + std::cout << "Done!" << std::endl; + + unsigned simulatedAtlasSize = 2 * imageFiles.size(); + + std::cout << "---------------------------------------" << std::endl; + std::cout << "3. Build statistical deformation model..." << std::endl; + if( BuildStatisticalDeformationModel(root, imageFiles, simulatedAtlasSize) != 0 ) + { + return -1; + } + std::cout << "Done. " << std::endl; + + // step 2.5: generate templaet images + ///////////////////////// + // generate simulated template with deformation field + std::cout << "---------------------------------------" << std::endl; + std::cout << "4. Generate simulated templates ... " << std::endl; + std::vector simulatedImageFiles = GenerateSimulatedData(root, imageFiles, simulatedAtlasSize); + std::cout << "Done. " << std::endl; + + ///////////////////////// + // delete subsampled deformation field + // if (!isDebug) if (false) + { + for( unsigned i = 0; i < imageFiles.size(); ++i ) + { + if( int(i) == root ) + { + continue; + } + char i_str[10], root_str[10]; + sprintf(i_str, "%03d", i); + sprintf(root_str, "%03d", root); + std::string subdeformationFieldFileName = std::string(i_str) + "_to_" + std::string(root_str) + + "deform_000_sub.nii.gz"; + subdeformationFieldFileName = outputFolder + subdeformationFieldFileName; + basicoperator->RemoveFile(subdeformationFieldFileName); + } + } + + ////////////////////////////// + // build the combinative tree + std::cout << "---------------------------------------" << std::endl; + std::cout << "5. Build combinative tree ... " << std::endl; + unsigned totalAtlasSize = imageFiles.size() + simulatedAtlasSize; + unsigned tree_size = totalAtlasSize; + vnl_vector tree(tree_size); // each tree + + // fill distance matrix + std::vector atlasfilenames; + for( unsigned i = 0; i < tree_size; ++i ) + { + std::string atlasfilename; + + if( i < imageFiles.size() ) + { + atlasfilename.append(imageFiles[i]); + } + else + { + atlasfilename.append(simulatedImageFiles[i - imageFiles.size()]); + } + + // atlasfilename = outputFolder + atlasfilename; + atlasfilenames.push_back(atlasfilename); + } + + tree = treeoperator->BuildCombinativeTree(root, atlasfilenames, tree_size, tree); // cout << "pass: + // BuildCombinativeTree " << endl; + + std::cout << "Done. " << std::endl; + + // write trained atlas to xml file + std::cout << "---------------------------------------" << std::endl; + std::cout << "6. Write to output file..." << std::endl; + itk::MABMISAtlasXMLFileWriter::Pointer atlasWriter = itk::MABMISAtlasXMLFileWriter::New(); + itk::MABMISAtlas atlas; + atlas.m_NumberAllAtlases = totalAtlasSize; + atlas.m_NumberSimulatedAtlases = simulatedAtlasSize; + atlas.m_AtlasDirectory = std::string(".") + FILESEP + outputxmlname; + atlas.m_Tree.resize(tree_size); + atlas.m_IsSimulatedImage.resize(tree_size); + atlas.m_RealImageIDs.resize(totalAtlasSize - simulatedAtlasSize); + atlas.m_SimulatedImageIDs.resize(simulatedAtlasSize); + for( unsigned i = 0; i < tree_size; ++i ) + { + atlas.m_Tree[i] = tree[i]; + } + atlas.m_AtlasFilenames.resize(atlasfilenames.size() ); + unsigned countRealImages = 0, countSimulatedImages = 0; + for( unsigned i = 0; i < atlasfilenames.size(); ++i ) + { + const size_t sep = atlasfilenames[i].find_last_of(FILESEP); + std::string fname = atlasfilenames[i]; + if( sep != std::string::npos ) + { + fname = atlasfilenames[i].substr(sep + 1, std::string::npos); +/* fname = atlasfilenames[i].substr(0, sep); + sep = fname.find_last_of(FILESEP); + if (sep != std::string::npos) + fname = atlasfilenames[i].substr(sep+1, std::string::npos); + else + fname = atlasfilenames[i]; + */ + } + else + { + fname = atlasfilenames[i]; + } + + atlas.m_AtlasFilenames[i] = fname; + + if( i < totalAtlasSize - simulatedAtlasSize ) + { + atlas.m_IsSimulatedImage[i] = false; + atlas.m_RealImageIDs[countRealImages++] = i; + } + else + { + atlas.m_IsSimulatedImage[i] = true; + atlas.m_SimulatedImageIDs[countSimulatedImages++] = i; + } + } + // atlas segmentation files + atlas.m_AtlasSegmentationFilenames.resize(segmentFiles.size() ); + for( unsigned i = 0; i < segmentFiles.size(); ++i ) + { + const size_t sep = segmentFiles[i].find_last_of(FILESEP); + std::string fname = segmentFiles[i]; + if( sep != std::string::npos ) + { + fname = segmentFiles[i].substr(sep + 1, std::string::npos); +/* fname = segmentFiles[i].substr(0, sep); + sep = fname.find_last_of(FILESEP); + if (sep != std::string::npos) + fname = segmentFiles[i].substr(sep+1, std::string::npos); + else + fname = segmentFiles[i]; +*/ + } + else + { + fname = segmentFiles[i]; + } + atlas.m_AtlasSegmentationFilenames[i] = fname; + } + + treeoperator->FindRoot(tree, tree_size, atlas.m_TreeRoot); + treeoperator->GetTreeHeight(tree, tree_size, atlas.m_TreeHeight); + + atlasWriter->SetObject(&atlas); + atlasWriter->SetFilename(outputFile); + atlasWriter->WriteFile(); + + std::cout << "Done." << std::endl; + + return EXIT_SUCCESS; +} + +void SearchRootAmongAtlases(std::vector imgfilenames, int & root) +{ + int tree_size = imgfilenames.size(); + + vnl_matrix distanceMatrix(tree_size, tree_size); + for( int i = 0; i < tree_size; ++i ) + { + for( int j = 0; j < tree_size; ++j ) + { + distanceMatrix[i][j] = 0.0; + } + } + + vnl_vector tree(tree_size); // each tree + + // fill distance matrix + imgoperator->PairwiseDistanceAmongImages(imgfilenames, tree_size, distanceMatrix); + if( isDebug ) + { + std::string distFileName = "dist_all_atlases.txt"; + basicoperator->SaveMatrix2File(distanceMatrix, tree_size, tree_size, distFileName); + } + + // build MST + tree = treeoperator->generateMSTFromMatrix(distanceMatrix, tree_size, tree); + + treeoperator->FindRoot(tree, tree_size, root); + + // save distance matrix and tree + if( isDebug ) + { + std::string treeFileName = "tree_all_atlases.txt"; + treeoperator->SaveTreeWithInfo(tree, tree_size, treeFileName); + } + + std::cout << "The root is " << root << ": " << imgfilenames[root] << ". " << std::endl; +} + +int RegistrationBetweenRootandAtlases(int root, std::vector imageFileNames, + std::vector iterations, double sigma) +{ + // get the path from the filename + const size_t sep = imageFileNames[0].find_last_of(FILESEP); + std::string outputFolder = ""; + + if( sep != std::string::npos ) + { + outputFolder = imageFileNames[0].substr(0, sep); + } + if( !outputFolder.empty() ) + { + outputFolder = outputFolder + FILESEP; + } + + const int atlas_size = imageFileNames.size(); + std::cout << "Register between root and atlases ... " << std::endl; + for( int i = 0; i < atlas_size; ++i ) + { + std::cout << i << ", "; + if( i == root ) + { + continue; + } + // if not the root image do registration + std::string fixedImageFileName = imageFileNames[root]; + std::string movingImageFileName = imageFileNames[i]; + + char i_str[10], root_str[10]; + sprintf(i_str, "%03d", i); + sprintf(root_str, "%03d", root); + std::string deformedImageFileName = std::string(i_str) + "_to_" + std::string(root_str) + RegisteredSuffix; + std::string deformationFieldFileName = std::string(i_str) + "_to_" + std::string(root_str) + DeformationSuffix; + + deformedImageFileName = outputFolder + deformedImageFileName; + deformationFieldFileName = outputFolder + deformationFieldFileName; + + if( !itksys::SystemTools::FileExists(deformationFieldFileName.c_str(), true) ) + { + // rough registration + + int val = regoperator->DiffeoDemonsRegistrationWithParameters( + fixedImageFileName, movingImageFileName, + deformedImageFileName, deformationFieldFileName, + sigma + 0.5, doHistMatchTrue, iterations); + if( val != 0 ) + { + std::cout << "Cannot perform registration between :" << std::endl; + std::cout << fixedImageFileName << " and " << std::endl; + std::cout << movingImageFileName << "." << std::endl; + std::cout << "Please verify the input images are correct." << std::endl; + } + } + + // --------------------------------------------------------------------- + // ------------- Generate the inverse deformation field + + DeformationFieldType::Pointer deformationField = nullptr; + dfoperator->ReadDeformationField(deformationFieldFileName, deformationField); + + DeformationFieldType::Pointer inversedDeformationField = DeformationFieldType::New(); + inversedDeformationField->SetRegions(deformationField->GetLargestPossibleRegion() ); + inversedDeformationField->SetSpacing(deformationField->GetSpacing() ); + inversedDeformationField->SetDirection(deformationField->GetDirection() ); + inversedDeformationField->SetOrigin(deformationField->GetOrigin() ); + inversedDeformationField->Allocate(); + + dfoperator->InverseDeformationField3D(deformationField, inversedDeformationField); + + std::string invDeformedImageFileName = std::string(root_str) + "_to_" + std::string(i_str) + RegisteredSuffix; + std::string invDeformationFileName = std::string(root_str) + "_to_" + std::string(i_str) + DeformationSuffix; + invDeformedImageFileName = outputFolder + invDeformedImageFileName; + invDeformationFileName = outputFolder + invDeformationFileName; + + std::string invDeformedSegmentFileName = std::string(root_str) + "_to_" + std::string(i_str) + SegmentationSuffix; + invDeformedSegmentFileName = outputFolder + invDeformedSegmentFileName; + + dfoperator->WriteDeformationField(invDeformationFileName, inversedDeformationField); + + // update + regoperator->DiffeoDemonsRegistrationWithInitialWithParameters( + movingImageFileName, fixedImageFileName, + invDeformationFileName, + invDeformedImageFileName, invDeformationFileName, + sigma, doHistMatch, iterations); + + if( isEvaluate ) + { + dfoperator->ApplyDeformationFieldAndWriteWithTypeWithFileNames( + fixedImageFileName, + invDeformationFileName, + invDeformedSegmentFileName, false); + } + } + std::cout << std::endl; + std::cout << "Done. " << std::endl; + + return 0; +} + +int BuildStatisticalDeformationModel(int root, std::vector imageFileNames, int simulatedAtlasSize) +{ + /////////////////////////////// + // do PCA simulation + // std::cout << "Build deformation field model ... "; + std::vector allDeformationFieldFileNames; + + // get the path from the filename + const size_t sep = imageFileNames[0].find_last_of(FILESEP); + std::string outputFolder = ""; + if( sep != std::string::npos ) + { + outputFolder = imageFileNames[0].substr(0, sep); + } + if( !outputFolder.empty() ) + { + outputFolder = outputFolder + FILESEP; + } + + int atlas_size = imageFileNames.size(); + + datasimulator->SetRoot(root); + datasimulator->SetAtlasSize(atlas_size); + datasimulator->SetSimulateSize(simulatedAtlasSize); + for( int i = 0; i < atlas_size; ++i ) + { + if( i == root ) + { + continue; + } + + char i_str[10], root_str[10]; + sprintf(i_str, "%03d", i); + sprintf(root_str, "%03d", root); + + std::string deformationFieldFileName = std::string(i_str) + "_to_" + std::string(root_str) + DeformationSuffix; + deformationFieldFileName = outputFolder + deformationFieldFileName; + allDeformationFieldFileNames.push_back(deformationFieldFileName); + } + + return datasimulator->DoPCATraining(allDeformationFieldFileNames, atlas_size - 1, imageFileNames, root); +} + +// return: simulated template file names +std::vector GenerateSimulatedData(int root, std::vector imageFiles, int simulatedAtlasSize) +{ + // output folder + const size_t sep = imageFiles[0].find_last_of(FILESEP); + std::string outputFolder = ""; + + if( sep != std::string::npos ) + { + outputFolder = imageFiles[0].substr(0, sep); + } + if( !outputFolder.empty() ) + { + outputFolder = outputFolder + FILESEP; + } + + std::vector simulateDeformationFieldFileNames(0); + std::vector simulateTemplateFileNames(0); + std::cout << "Generating simulated images: " << std::endl; + for( int i = 0; i < simulatedAtlasSize; ++i ) + { + std::cout << i << ", "; + + std::string index_string = basicoperator->myitoa( i, 3 ); + + std::string simulateDeformationFieldFileName = "simulated_deform_" + index_string + ".nii.gz"; + std::string simulateTemplateFileName = "simulated_cbq_" + index_string + ".nii.gz"; + simulateDeformationFieldFileName = outputFolder + simulateDeformationFieldFileName; + simulateTemplateFileName = outputFolder + simulateTemplateFileName; + + simulateDeformationFieldFileNames.push_back(simulateDeformationFieldFileName); + simulateTemplateFileNames.push_back(simulateTemplateFileName); + + // load simulated deformation field + DeformationFieldType::Pointer df = nullptr; + dfoperator->ReadDeformationField(simulateDeformationFieldFileNames[i], df); + + DeformationFieldType::Pointer invdf = DeformationFieldType::New(); + invdf->SetRegions(df->GetLargestPossibleRegion() ); + invdf->SetSpacing(df->GetSpacing() ); + invdf->SetDirection(df->GetDirection() ); + invdf->SetOrigin(df->GetOrigin() ); + invdf->Allocate(); + + dfoperator->InverseDeformationField3D(df, invdf); + InternalImageType::Pointer rootImg = nullptr; + imgoperator->ReadImage(imageFiles[root], rootImg); + const itk::ImageIOBase::IOComponentType ioType = imgoperator->GetIOPixelType(imageFiles[root]); + InternalImageType::Pointer simImg = nullptr; + dfoperator->ApplyDeformationField(rootImg, invdf, simImg, true); + imgoperator->WriteImage(simulateTemplateFileNames[i], simImg, ioType); + } + std::cout << std::endl; + + return simulateTemplateFileNames; +} diff --git a/Training.h b/Training.h new file mode 100644 index 0000000..2f6392c --- /dev/null +++ b/Training.h @@ -0,0 +1,9 @@ +#ifndef __Training_h +#define __Training_h + +#include "commonMABMIS.h" + +int Training(itk::MABMISImageData* trainingData, std::string outputFile, + std::vector iterations, double sigma); + +#endif //__Training_h diff --git a/commonMABMIS.cxx b/commonMABMIS.cxx new file mode 100644 index 0000000..71b0f61 --- /dev/null +++ b/commonMABMIS.cxx @@ -0,0 +1,71 @@ +#include "commonMABMIS.h" + +#include + +std::string ReplacePathSepForOS( const std::string & input ) +{ + std::string output = input; +#ifdef _WIN32 + std::replace(output.begin(), output.end(), '/', FILESEP); +#else + std::replace(output.begin(), output.end(), '\\', FILESEP); +#endif + return output; +} + +std::string +GetExtension(const std::string & filename) +{ + std::string fileExt( itksys::SystemTools::GetFilenameLastExtension(filename) ); + //If the last extension is .gz, then need to pull off 2 extensions. + //.gz is the only valid compression extension. + if ( fileExt == std::string(".gz") ) + { + fileExt = itksys::SystemTools::GetFilenameLastExtension( + itksys::SystemTools::GetFilenameWithoutLastExtension(filename) ); + fileExt += ".gz"; + } + return ( fileExt ); +} + +std::string +GetRootName(const std::string & filename) +{ + const std::string fileExt = GetExtension(filename); + + // Create a base filename + // i.e Image.hdr --> Image + if ( fileExt.length() > 0 //Ensure that an extension was found + && filename.length() > fileExt.length() //Ensure that the filename does + // not contain only the extension + ) + { + const std::string::size_type it = filename.find_last_of(fileExt); + const std::string baseName( filename, 0, it - ( fileExt.length() - 1 ) ); + return ( baseName ); + } + //Default to return same as input when the extension is nothing (Analyze) + return ( filename ); +} + +// demons registration parameters +// int iterInResolutions[4][3]={{5,3,2},{10,5,5},{15,10,5},{20,15,10}}; +// int itereach = 2; +// int itereach0 = 0; +// int itereach1 = 1; +// int itereach2 = 2; +// int itereach3 = 3; +// double sigmaDef = 1.5; +// double sigmaDef10 = 1.0; +// double sigmaDef15 = 1.5; +// double sigmaDef20 = 2.0; +// double sigmaDef25 = 2.5; +// double sigmaDef30 = 3.0; +// double sigmaDef35 = 3.5; + +DataSimulatorType::Pointer datasimulator = DataSimulatorType::New(); +ImageOperationFilterType::Pointer imgoperator = ImageOperationFilterType::New(); +DeformationFieldOperationFilterType::Pointer dfoperator = DeformationFieldOperationFilterType::New(); +ImageRegistrationFilterType::Pointer regoperator = ImageRegistrationFilterType::New(); +TreeOperationType::Pointer treeoperator = TreeOperationType::New(); +BasicOperationFilterType::Pointer basicoperator = BasicOperationFilterType::New(); diff --git a/commonMABMIS.h b/commonMABMIS.h new file mode 100644 index 0000000..b5594d0 --- /dev/null +++ b/commonMABMIS.h @@ -0,0 +1,172 @@ +#ifndef __commonMABMIS_h +#define __commonMABMIS_h + +// for math +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef WIN32 +constexpr char FILESEP = '\\'; +#else +constexpr char FILESEP = '/'; +#endif + +// including itksys::SystemTools::MakeDirectory(char*) +#include +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkImageToImageFilter.h" +#include "itkImageRegionIterator.h" + +constexpr unsigned int ImageDimension = 3; + +// basic itk +#include "itkVector.h" + +// registration +#include "itkImageRegistrationMethod.h" +#include "itkSymmetricForcesDemonsRegistrationFilter.h" + +// interpolator +#include "itkLinearInterpolateImageFunction.h" +#include "itkNearestNeighborInterpolateImageFunction.h" +#include "itkNearestNeighborExtrapolateImageFunction.h" + +// filter +#include "itkResampleImageFilter.h" + +#include "itkCastImageFilter.h" +#include "itkWarpImageFilter.h" + +#include "itkHistogramMatchingImageFilter.h" +#include "itkAddImageFilter.h" +#include "itkDivideImageFilter.h" +#include "itkMultiplyImageFilter.h" +#include "itkWarpVectorImageFilter.h" +#include "itkInverseDisplacementFieldImageFilter.h" + +// for affine transformation +#include "itkTransform.h" +#include "itkAffineTransform.h" +#include "itkImageRegistrationMethod.h" + +// for Diffeomorphic Demons +#include +#include +#include + +// To include all related header files +#include "itkMABMISBasicOperationFilter.h" +#include "itkMABMISImageOperationFilter.h" +#include "itkMABMISDeformationFieldFilter.h" +#include "itkMABMISSimulateData.h" +#include "itkMABMISImageRegistrationFilter.h" +#include "itkMABMISTreeOperation.h" +#include "itkMABMISAtlasXMLFile.h" + +typedef double CoordinateRepType; + +// basic data type +typedef unsigned char CharPixelType; // for image IO usage +typedef float FloatPixelType; // for +typedef int IntPixelType; +typedef short ShortPixelType; +typedef float InternalPixelType; // for internal processing usage +typedef itk::Vector VectorPixelType; + +// basic image type +typedef itk::Image CharImageType; +typedef itk::Image IntImageType; +typedef itk::Image ShortImageType; +typedef itk::Image FloatImageType; +typedef itk::Image InternalImageType; +typedef itk::Image DeformationFieldType; + +// basic iterator type +typedef itk::ImageRegionIterator DeformationFieldIteratorType; +typedef itk::ImageRegionIterator InternalImageIteratorType; +typedef itk::ImageRegionIterator CharImageIteratorType; + +// basic image reader/writer related type +typedef itk::ImageFileReader CharImageReaderType; +typedef itk::ImageFileReader InternalImageReaderType; +typedef itk::ImageFileWriter InternalImageWriterType; + +typedef itk::WarpImageFilter InternalWarpFilterType; +typedef itk::ImageFileWriter CharImageWriterType; +typedef itk::ImageFileWriter IntImageWriterType; +typedef itk::ImageFileWriter FloatImageWriterType; +typedef itk::ImageFileWriter ShortImageWriterType; + +typedef itk::ImageFileReader DeformationFieldReaderType; +typedef itk::ImageFileWriter DeformationFieldWriterType; + +////////////////////////////////////////////////////////////////////////////// +// image filter type +typedef itk::HistogramMatchingImageFilter InternalHistMatchFilterType; + +//////////////////////////////////////////////////////////////////////////// +// operation on deformation fields +typedef itk::WarpVectorImageFilter WarpVectorFilterType; +typedef itk::AddImageFilter AddImageFilterType; + +typedef itk::Vector ShortVectorPixelType; +typedef itk::Image ShortDeformationFieldType; +typedef itk::ImageFileWriter ShortDeformationFieldWriterType; + +std::string ReplacePathSepForOS(const std::string & input); + +std::string GetRootName(const std::string & filename); + +std::string GetExtension(const std::string & filename); + + +// global bool variables to adjust the procedure +constexpr bool isEvaluate = false; // if false, we do not know the ground-truth of labels +constexpr bool isDebug = false; // false;//true; // if true, print out more information +constexpr int localPatchSize = 1; // (2r+1)*(2r+1)*(2r+1) is the volume of local patch +constexpr bool doHistMatch = true; +constexpr bool doHistMatchTrue = true; +constexpr bool doHistMatchFalse = false; + +//demons registration parameters +//extern int iterInResolutions[4][3]; +//extern int itereach; +//extern int itereach0; +//extern int itereach1; +//extern int itereach2; +//extern int itereach3; +//extern double sigmaDef; +//extern double sigmaDef10; +//extern double sigmaDef15; +//extern double sigmaDef20; +//extern double sigmaDef25; +//extern double sigmaDef30; +//extern double sigmaDef35; + +typedef itk::Statistics::MABMISSimulateData DataSimulatorType; +extern DataSimulatorType::Pointer datasimulator; +typedef itk::Statistics::MABMISImageOperationFilter ImageOperationFilterType; +extern ImageOperationFilterType::Pointer imgoperator; +typedef itk::Statistics::MABMISDeformationFieldFilter DeformationFieldOperationFilterType; +extern DeformationFieldOperationFilterType::Pointer dfoperator; +typedef itk::Statistics::MABMISImageRegistrationFilter ImageRegistrationFilterType; +extern ImageRegistrationFilterType::Pointer regoperator; +typedef itk::Statistics::MABMISTreeOperation TreeOperationType; +extern TreeOperationType::Pointer treeoperator; +typedef itk::Statistics::MABMISBasicOperationFilter BasicOperationFilterType; +extern BasicOperationFilterType::Pointer basicoperator; + +constexpr char ImageSuffix[] = "_cbq_000.nii.gz"; +constexpr char DeformationSuffix[] = "_deform_000.nii.gz"; +constexpr char SegmentationSuffix[] = "_seg_000.nii.gz"; +constexpr char RegisteredSuffix[] = "_cbq_reg.nii.gz"; + +#endif //__commonMABMIS_h diff --git a/itkMABMISAtlasXMLFile.cxx b/itkMABMISAtlasXMLFile.cxx index 428ab40..a1afca2 100644 --- a/itkMABMISAtlasXMLFile.cxx +++ b/itkMABMISAtlasXMLFile.cxx @@ -9,18 +9,6 @@ #include #include -#define RAISE_EXCEPTION(s) \ - { \ - ExceptionObject \ - exception( \ - __FILE__, \ - __LINE__); \ - exception. \ - SetDescription( \ - s); \ - throw exception; \ - } - #define SPACES " \t\r\n" inline std::string trim_right (const std::string & s, const std::string & t = SPACES) @@ -31,14 +19,14 @@ inline std::string trim_right (const std::string & s, const std::string & t = SP return ""; else return d.erase (d.find_last_not_of (t) + 1) ; -} +} inline std::string trim_left (const std::string & s, const std::string & t = SPACES) { std::string d (s); return d.erase (0, s.find_first_not_of (t)) ; -} - +} + inline std::string trim (const std::string & s, const std::string & t = SPACES) { std::string d (s); @@ -121,7 +109,7 @@ MABMISImageDataXMLFileReader::EndElement(const char *name) { if( m_ImageData == 0 ) { - RAISE_EXCEPTION("The xml format is not right. Please check it!"); + itkGenericExceptionMacro("The xml format is not right. Please check it!"); } } if( itksys::SystemTools::Strucmp(name, "DataPath") == 0 ) @@ -138,18 +126,6 @@ MABMISImageDataXMLFileReader::EndElement(const char *name) } } -static std::string StripLastNewline(const std::string & input) -{ - std::string output = input; - size_t last_element = input.size() - 1; - - if( input[last_element] == '\n' ) - { - output = input.substr(0, last_element); - } - return output; -} - void MABMISImageDataXMLFileReader::CharacterDataHandler(const char *inData, int inLength) { @@ -158,7 +134,6 @@ MABMISImageDataXMLFileReader::CharacterDataHandler(const char *inData, int inLen { m_CurCharacterData = m_CurCharacterData + inData[i]; } - //m_CurCharacterData = StripLastNewline(m_CurCharacterData); m_CurCharacterData = trim(m_CurCharacterData); } @@ -184,13 +159,13 @@ MABMISAtlasXMLFileReader::GenerateOutputInformation() this->parse(); // validate the results are right. - int numRealImages = this->m_OutputObject->m_NumberAllAtlases - this->m_OutputObject->m_NumberSimulatedAtlases; - int numSimImages = this->m_OutputObject->m_NumberSimulatedAtlases; + unsigned numRealImages = this->m_OutputObject->m_NumberAllAtlases - this->m_OutputObject->m_NumberSimulatedAtlases; + unsigned numSimImages = this->m_OutputObject->m_NumberSimulatedAtlases; this->m_OutputObject->m_SimulatedImageIDs.resize(numSimImages); this->m_OutputObject->m_RealImageIDs.resize(numRealImages); - int countRealImages = 0, countSimImages = 0; - for( int n = 0; n < this->m_OutputObject->m_IsSimulatedImage.size(); n++ ) + unsigned countRealImages = 0, countSimImages = 0; + for( unsigned n = 0; n < this->m_OutputObject->m_IsSimulatedImage.size(); n++ ) { if( this->m_OutputObject->m_IsSimulatedImage[n] ) { @@ -229,7 +204,7 @@ MABMISAtlasXMLFileReader::StartElement( const char *name, const char * *atts ) else if( itksys::SystemTools::Strucmp(name, "ATLAS") == 0 ) { int i = 0; - int id; + int id = -4; std::string fname = ""; std::string segFName = ""; bool isSimulated = false; @@ -255,9 +230,9 @@ MABMISAtlasXMLFileReader::StartElement( const char *name, const char * *atts ) } else if( itksys::SystemTools::Strucmp(key.c_str(), "SIMULATED") == 0 ) { - if( itksys::SystemTools::Strucmp(value.c_str(), "TRUE") || - itksys::SystemTools::Strucmp(value.c_str(), "1") || - itksys::SystemTools::Strucmp(value.c_str(), "T") ) + if (itksys::SystemTools::Strucmp(value.c_str(), "TRUE") == 0 || + itksys::SystemTools::Strucmp(value.c_str(), "1") == 0 || + itksys::SystemTools::Strucmp(value.c_str(), "T") == 0) { isSimulated = true; } @@ -279,7 +254,7 @@ MABMISAtlasXMLFileReader::StartElement( const char *name, const char * *atts ) else if( itksys::SystemTools::Strucmp(name, "Tree") == 0 ) { m_Atlas->m_Tree.resize(m_Atlas->m_NumberAllAtlases); - for( int n = 0; n < m_Atlas->m_Tree.size(); n++ ) + for( unsigned n = 0; n < m_Atlas->m_Tree.size(); n++ ) { m_Atlas->m_Tree[n] = -1; } @@ -287,7 +262,7 @@ MABMISAtlasXMLFileReader::StartElement( const char *name, const char * *atts ) else if( itksys::SystemTools::Strucmp(name, "Node") == 0 ) { int i = 0; - int id, parentID; + int id = -2, parentID = -3; ; while( atts[i] ) { @@ -321,7 +296,7 @@ MABMISAtlasXMLFileReader::EndElement(const char *name) { if( m_Atlas == 0 ) { - RAISE_EXCEPTION("The xml format is not right"); + itkGenericExceptionMacro("The xml format is not right"); } } if( itksys::SystemTools::Strucmp(name, "AtlasDirectory") == 0 ) @@ -363,10 +338,73 @@ MABMISAtlasXMLFileReader::CharacterDataHandler(const char *inData, int inLength) { m_CurCharacterData = m_CurCharacterData + inData[i]; } - //m_CurCharacterData = StripLastNewline(m_CurCharacterData); m_CurCharacterData = trim(m_CurCharacterData); } +///------------------------------------------------------------ +// ---- class MABMISImageDataXMLFileWriter ------------------ +///------------------------------------------------------------ + +int +MABMISImageDataXMLFileWriter::CanWriteFile(const char *) +{ + return true; +} + +int +MABMISImageDataXMLFileWriter::WriteFile() +{ + // sanity checks + if( m_InputObject == 0 ) + { + itkGenericExceptionMacro("No MABMISAtlas to Write"); + } + if( m_Filename.length() == 0 ) + { + itkGenericExceptionMacro("No filename given"); + } + std::ofstream output( m_Filename.c_str() ); + if( output.fail() ) + { + itkGenericExceptionMacro("Can't Open "<< m_Filename); + } + + WriteStartElement("?xml version=\"1.0\"?", output); + output << std::endl; + WriteStartElement("!DOCTYPE MABMISImageData", output); + output << std::endl; + + WriteStartElement("MABMISImageData", output); + output << std::endl; + + WriteStartElement("DataPath", output); + output << this->m_InputObject->m_DataDirectory; + WriteEndElement("DataPath", output); + output << std::endl; + + WriteStartElement("NumberOfImageData", output); + output << this->m_InputObject->m_NumberImageData; + WriteEndElement("NumberOfImageData", output); + output << std::endl; + + WriteStartElement("ImageDataSets", output); + output << std::endl; + for( unsigned n = 0; n < this->m_InputObject->m_ImageFileNames.size(); n++ ) + { + output << "m_InputObject->m_ImageFileNames[n] << "\""; + output << " SegmentationFile=\"" << this->m_InputObject->m_SegmentationFileNames[n] << "\"/>" << std::endl; + } + + WriteEndElement("ImageDataSets", output); + output << std::endl; + + WriteEndElement("MABMISImageData", output); + output << std::endl; + output.close(); + + return 0; +} + ///------------------------------------------------------------ // ---- class MABMISAtlasXMLFileWriter ---------------------- ///------------------------------------------------------------ @@ -374,7 +412,7 @@ MABMISAtlasXMLFileReader::CharacterDataHandler(const char *inData, int inLength) int MABMISAtlasXMLFileWriter::CanWriteFile( const char *itkNotUsed(name) ) { - return true; // not sure what else to say + return true; } int @@ -383,20 +421,16 @@ MABMISAtlasXMLFileWriter::WriteFile() // sanity checks if( m_InputObject == 0 ) { - std::string errmsg("No MABMISAtlas to Write"); - RAISE_EXCEPTION(errmsg); + itkGenericExceptionMacro("No MABMISAtlas to Write"); } if( m_Filename.length() == 0 ) { - std::string errmsg("No filename given"); - RAISE_EXCEPTION(errmsg); + itkGenericExceptionMacro("No filename given"); } std::ofstream output( m_Filename.c_str() ); if( output.fail() ) { - std::string errmsg("Can't Open "); - errmsg += m_Filename; - RAISE_EXCEPTION(errmsg); + itkGenericExceptionMacro("Can't Open "<< m_Filename); } WriteStartElement("?xml version=\"1.0\"?", output); @@ -430,7 +464,7 @@ MABMISAtlasXMLFileWriter::WriteFile() // atlas file list WriteStartElement("AtlasFileList", output); output << std::endl; - for( int n = 0; n < this->m_InputObject->m_AtlasFilenames.size(); n++ ) + for( unsigned n = 0; n < this->m_InputObject->m_AtlasFilenames.size(); n++ ) { output << "m_InputObject->m_Tree.size(); ++i ) + for( unsigned i = 0; i < this->m_InputObject->m_Tree.size(); ++i ) { output << " #include -#ifdef _WIN32 -static const char FILESEP = '\\'; -#else -static const char FILESEP = '/'; -#endif - namespace itk { class MABMISImageData @@ -18,7 +12,7 @@ class MABMISImageData public: std::string m_DataDirectory; std::string m_OutputDirectory; - int m_NumberImageData; + unsigned m_NumberImageData; std::vector m_ImageFileNames; std::vector m_SegmentationFileNames; @@ -32,7 +26,7 @@ class MABMISImageData } ~MABMISImageData() { - }; + } }; class MABMISAtlas @@ -40,8 +34,8 @@ class MABMISAtlas public: std::string m_AtlasDirectory; - int m_NumberAllAtlases; - int m_NumberSimulatedAtlases; + unsigned m_NumberAllAtlases; + unsigned m_NumberSimulatedAtlases; std::vector m_AtlasFilenames; std::vector m_AtlasSegmentationFilenames; @@ -51,7 +45,7 @@ class MABMISAtlas std::vector m_SimulatedImageIDs; std::vector m_Tree; - int m_TreeSize; + unsigned int m_TreeSize; int m_TreeRoot; int m_TreeHeight; @@ -65,14 +59,14 @@ class MABMISAtlas m_TreeSize = 0; m_TreeRoot = -1; m_TreeHeight = 0; - }; + } ~MABMISAtlas() { - }; + } }; class MABMISImageDataXMLFileReader : - public XMLReader + public XMLReader { public: /** Standard typedefs */ @@ -87,7 +81,7 @@ class MABMISImageDataXMLFileReader : itkNewMacro(Self); public: /** Determine if a file can be read */ - virtual int CanReadFile(const char *name); + int CanReadFile(const char *name) ITK_OVERRIDE; protected: MABMISImageDataXMLFileReader() @@ -95,32 +89,62 @@ class MABMISImageDataXMLFileReader : m_OutputObject = new itk::MABMISImageData; m_ImageCount = 0; m_ImageData = 0; - }; + } virtual ~MABMISImageDataXMLFileReader() { - }; + } - virtual void StartElement(const char *name, const char * *atts); + void StartElement(const char *name, const char * *atts) ITK_OVERRIDE; - virtual void EndElement(const char *name); + void EndElement(const char *name) ITK_OVERRIDE; - virtual void CharacterDataHandler(const char *inData, int inLength); + void CharacterDataHandler(const char *inData, int inLength) ITK_OVERRIDE; private: - MABMISImageDataXMLFileReader(const Self &); // purposely not - // implemented - void operator=(const Self &); // purposely not - - // implemented + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISImageDataXMLFileReader); MABMISImageData* m_ImageData; std::string m_CurCharacterData; - int m_ImageCount; + unsigned m_ImageCount; +}; + +class MABMISImageDataXMLFileWriter : + public XMLWriterBase +{ +public: + /** standard typedefs */ + typedef XMLWriterBase Superclass; + typedef MABMISImageDataXMLFileWriter Self; + typedef SmartPointer Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(MABMISImageDataXMLFileWriter, XMLWriterBase ); + + /** Test whether a file is writable. */ + int CanWriteFile(const char *name) ITK_OVERRIDE; + + /** Actually write out the file in question */ + int WriteFile() ITK_OVERRIDE; + +protected: + MABMISImageDataXMLFileWriter() + { + } + + virtual ~MABMISImageDataXMLFileWriter() + { + } + +private: + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISImageDataXMLFileWriter); }; class MABMISAtlasXMLFileReader : - public XMLReader + public XMLReader { public: /** Standard typedefs */ @@ -134,41 +158,37 @@ class MABMISAtlasXMLFileReader : /** Method for creation through the object factory. */ itkNewMacro(Self); - virtual void GenerateOutputInformation(); + void GenerateOutputInformation() ITK_OVERRIDE; public: /** Determine if a file can be read */ - virtual int CanReadFile(const char *name); + int CanReadFile(const char *name) ITK_OVERRIDE; protected: MABMISAtlasXMLFileReader() { m_OutputObject = new itk::MABMISAtlas; m_Atlas = 0; - }; + } virtual ~MABMISAtlasXMLFileReader() { - }; + } - virtual void StartElement(const char *name, const char * *atts); + void StartElement(const char *name, const char * *atts) ITK_OVERRIDE; - virtual void EndElement(const char *name); + void EndElement(const char *name) ITK_OVERRIDE; - virtual void CharacterDataHandler(const char *inData, int inLength); + void CharacterDataHandler(const char *inData, int inLength) ITK_OVERRIDE; private: - MABMISAtlasXMLFileReader(const Self &); // purposely not - // implemented - void operator=(const Self &); // purposely not - - // implemented + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISAtlasXMLFileReader); MABMISAtlas* m_Atlas; std::string m_CurCharacterData; }; class MABMISAtlasXMLFileWriter : - public XMLWriterBase + public XMLWriterBase { public: /** standard typedefs */ @@ -183,10 +203,10 @@ class MABMISAtlasXMLFileWriter : itkTypeMacro(MABMISAtlasXMLFileWriter, XMLWriterBase ); /** Test whether a file is writable. */ - virtual int CanWriteFile(const char *name); + int CanWriteFile(const char *name) ITK_OVERRIDE; /** Actually write out the file in question */ - virtual int WriteFile(); + int WriteFile() ITK_OVERRIDE; protected: MABMISAtlasXMLFileWriter() @@ -198,11 +218,7 @@ class MABMISAtlasXMLFileWriter : } private: - MABMISAtlasXMLFileWriter(const Self &); // purposely not - // implemented - void operator=(const Self &); // purposely not - - // implemented + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISAtlasXMLFileWriter); }; } diff --git a/itkMABMISBasicOperationFilter.h b/itkMABMISBasicOperationFilter.h index 48cd893..a2c3c37 100644 --- a/itkMABMISBasicOperationFilter.h +++ b/itkMABMISBasicOperationFilter.h @@ -1,14 +1,7 @@ #ifndef __itkMABMISBasicOperationFilter_h #define __itkMABMISBasicOperationFilter_h -#include -#include - -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" -#include "itkImageRegionIterator.h" - -#define ImageDimension 3 +#include "commonMABMIS.h" namespace itk { @@ -37,18 +30,18 @@ class MABMISBasicOperationFilter : public ImageToImageFilter matrix, int iSize, int jSize, std::string martixFileName); private: - MABMISBasicOperationFilter(const Self &); // purposely not implemented - void operator=(const Self &); // purposely not implemented + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISBasicOperationFilter); protected: MABMISBasicOperationFilter(); ~MABMISBasicOperationFilter(); - void PrintSelf(std::ostream& os, Indent indent) const; + void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE; }; } // namespace itk } // namespace Statistics diff --git a/itkMABMISBasicOperationFilter.hxx b/itkMABMISBasicOperationFilter.hxx index da3756e..7d73de6 100644 --- a/itkMABMISBasicOperationFilter.hxx +++ b/itkMABMISBasicOperationFilter.hxx @@ -2,6 +2,7 @@ #define __itkMABMISBasicOperationFilter_hxx #include "itkMABMISBasicOperationFilter.h" +#include namespace itk { @@ -32,20 +33,13 @@ void MABMISBasicOperationFilter ::RemoveFile(std::string filename) { + int retVal; if( itksys::SystemTools::FileExists(filename.c_str(), true) ) { -#if defined(_WIN32) && !defined(__CYGWIN__) - std::ostringstream oss(std::ostringstream::out); - oss << "del /F " << filename << std::endl; - // std::cerr << oss.str().c_str(); - system(oss.str().c_str() ); - -#else - std::ostringstream oss(std::ostringstream::out); - oss << "rm -f " << filename << std::endl; - // std::cerr << oss.str().c_str(); - system(oss.str().c_str() ); -#endif + if (!itksys::SystemTools::RemoveFile(filename)) + { + std::cout << " Error deleting " << filename << std::endl; + } } return; } @@ -84,26 +78,13 @@ MABMISBasicOperationFilter } template -void +std::string MABMISBasicOperationFilter -::myitoa(int num, std::string& str, int digit) +::myitoa(int number, int digits) { - str = ""; - - if( num < 10 ) - { - str.append("00"); - } - else if( num < 100 ) - { - str.append("0"); - } - - std::stringstream st; - st << num; - str.append(st.str() ); - - return; + std::stringstream ss; + ss << std::setw(digits) << std::setfill('0') << number; + return ss.str(); } template diff --git a/itkMABMISDeformationFieldFilter.h b/itkMABMISDeformationFieldFilter.h index 2fa87a5..cdde7fa 100644 --- a/itkMABMISDeformationFieldFilter.h +++ b/itkMABMISDeformationFieldFilter.h @@ -1,14 +1,8 @@ #ifndef __itkMABMISDeformationFieldFilter_h #define __itkMABMISDeformationFieldFilter_h -#include -#include +#include "commonMABMIS.h" -#include -#include -#include - -#include "itkImage.h" #include "itkWarpImageFilter.h" #include "itkAddImageFilter.h" #include "itkDivideImageFilter.h" @@ -21,8 +15,6 @@ #include "itkMABMISImageOperationFilter.h" #include "itkWarpVectorImageFilter.h" -#define ImageDimension 3 - namespace itk { namespace Statistics @@ -107,9 +99,6 @@ class MABMISDeformationFieldFilter : public ImageToImageFilter WarpVectorFilterType::Pointer vectorWarper = WarpVectorFilterType::New(); vectorWarper->SetInput( input ); - vectorWarper->SetDeformationField( deformationField ); + vectorWarper->SetDisplacementField( deformationField ); vectorWarper->SetOutputOrigin(deformationField->GetOrigin() ); vectorWarper->SetOutputSpacing(deformationField->GetSpacing() ); vectorWarper->SetOutputDirection(deformationField->GetDirection() ); @@ -119,7 +119,7 @@ MABMISDeformationFieldFilter warper->SetOutputSpacing( movingImage->GetSpacing() ); warper->SetOutputOrigin( movingImage->GetOrigin() ); warper->SetOutputDirection( movingImage->GetDirection() ); - warper->SetDeformationField( deformationField ); + warper->SetDisplacementField( deformationField ); try { warper->Update(); @@ -142,7 +142,7 @@ MABMISDeformationFieldFilter warper->SetOutputSpacing( movingImage->GetSpacing() ); warper->SetOutputOrigin( movingImage->GetOrigin() ); warper->SetOutputDirection( movingImage->GetDirection() ); - warper->SetDeformationField( deformationField ); + warper->SetDisplacementField( deformationField ); try { warper->Update(); @@ -157,28 +157,6 @@ MABMISDeformationFieldFilter } } -// apply deformation field on image and write deformed image -template -void -MABMISDeformationFieldFilter -::ApplyDeformationFieldAndWriteWithFileNames(std::string movingImageName, std::string deformationFieldFileName, - std::string deformedImageName, bool isLinearInterpolator) -{ - typename ImageOperationType::Pointer imageoperator = ImageOperationType::New(); - - DeformationFieldType::Pointer deformationField = DeformationFieldType::New(); - ReadDeformationField(deformationFieldFileName, deformationField); - - InternalImageType::Pointer movingImage = InternalImageType::New(); - imageoperator->ReadImage(movingImageName, movingImage); - - InternalImageType::Pointer deformedImage = InternalImageType::New(); - - ApplyDeformationField(movingImage, deformationField, deformedImage, isLinearInterpolator); - - imageoperator->WriteImage(deformedImageName, deformedImage); -} - template void MABMISDeformationFieldFilter @@ -186,33 +164,8 @@ MABMISDeformationFieldFilter std::string deformationFieldFileName, std::string deformedImageFileName, bool isLinear) { - itk::ImageIOBase::Pointer imageIO; + const itk::ImageIOBase::IOComponentType ioType = imgoperator->GetIOPixelType(movingImageFileName); - try - { - imageIO = itk::ImageIOFactory::CreateImageIO(movingImageFileName.c_str(), itk::ImageIOFactory::ReadMode); - if( imageIO ) - { - imageIO->SetFileName(movingImageFileName); - imageIO->ReadImageInformation(); - } - else - { - std::cout << "Could not read the image information of " << movingImageFileName << "." << std::endl; - exit( EXIT_FAILURE ); - } - } - catch( itk::ExceptionObject& err ) - { - std::cout << "Could not read the image information of " << movingImageFileName << "." << std::endl; - std::cout << err << std::endl; - exit( EXIT_FAILURE ); - } - // {UNKNOWNCOMPONENTTYPE,UCHAR,CHAR,USHORT,SHORT,UINT,INT,ULONG,LONG,FLOAT,DOUBLE} IOComponentType; - // 0 1 2 3 4 5 6 7 8 9 10 - const int input_type = imageIO->GetComponentType(); // 9:float, 10:double - - // DeformationFieldType::Pointer deformationField = DeformationFieldType::New(); ReadDeformationField(deformationFieldFileName, deformationField); @@ -223,31 +176,9 @@ MABMISDeformationFieldFilter ApplyDeformationField(movingImage, deformationField, deformedImage, isLinear); - if( input_type == 1 ) // UCHAR - { - imgoperator->WriteImageUCHAR(deformedImageFileName, deformedImage); - } - else if( input_type == 4 ) // SHORT - { - imgoperator->WriteImageSHORT(deformedImageFileName, deformedImage); - } - else if( input_type == 6 ) // INT - { - imgoperator->WriteImageINT(deformedImageFileName, deformedImage); - } - else if( input_type == 9 || input_type == 10 ) // FLOAT || DOUBLE - { - imgoperator->WriteImageFLOAT(deformedImageFileName, deformedImage); - } - else - { - std::cerr << "input_type type not supported " << input_type << std::endl; - } + imgoperator->WriteImage(deformedImageFileName, deformedImage, ioType); } -// end -// - template void MABMISDeformationFieldFilter @@ -282,7 +213,7 @@ MABMISDeformationFieldFilter ry = sampleRate; rz = sampleRate; - DeformationFieldType::Pointer originImage = 0; + DeformationFieldType::Pointer originImage = nullptr; ReadDeformationField(deformationFieldFileName, originImage); @@ -447,7 +378,7 @@ MABMISDeformationFieldFilter ry = sampleRate; rz = sampleRate; - DeformationFieldType::Pointer originImage = 0; + DeformationFieldType::Pointer originImage = nullptr; ReadDeformationField(deformationFieldFileName, originImage); diff --git a/itkMABMISImageOperationFilter.h b/itkMABMISImageOperationFilter.h index f563963..558f4e4 100644 --- a/itkMABMISImageOperationFilter.h +++ b/itkMABMISImageOperationFilter.h @@ -1,17 +1,10 @@ #ifndef __itkMABMISImageOperationFilter_h #define __itkMABMISImageOperationFilter_h -#include -#include +#include "commonMABMIS.h" -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" -#include "itkImageRegionIterator.h" #include "itkCastImageFilter.h" - -#include "itkImage.h" - -#define ImageDimension 3 +#include "itkWarpImageFilter.h" namespace itk { @@ -42,32 +35,18 @@ class MABMISImageOperationFilter : public ImageToImageFilter VectorPixelType; // basic image type - typedef itk::Image CharImageType; - typedef itk::Image IntImageType; - typedef itk::Image ShortImageType; - typedef itk::Image FloatImageType; typedef itk::Image InternalImageType; typedef itk::Image DeformationFieldType; // basic iterator type typedef itk::ImageRegionIterator DeformationFieldIteratorType; typedef itk::ImageRegionIterator InternalImageIteratorType; - typedef itk::ImageRegionIterator CharImageIteratorType; // basic image reader/writer related type - typedef itk::ImageFileReader CharImageReaderType; typedef itk::ImageFileReader InternalImageReaderType; typedef itk::ImageFileWriter InternalImageWriterType; - typedef itk::CastImageFilter Internal2CharCastFilterType; - typedef itk::CastImageFilter Internal2FloatCastFilterType; - typedef itk::CastImageFilter Internal2IntCastFilterType; - typedef itk::CastImageFilter Internal2ShortCastFilterType; typedef itk::WarpImageFilter InternalWarpFilterType; - typedef itk::ImageFileWriter CharImageWriterType; - typedef itk::ImageFileWriter IntImageWriterType; - typedef itk::ImageFileWriter FloatImageWriterType; - typedef itk::ImageFileWriter ShortImageWriterType; typedef itk::ImageFileReader DeformationFieldReaderType; typedef itk::ImageFileWriter DeformationFieldWriterType; @@ -78,34 +57,31 @@ class MABMISImageOperationFilter : public ImageToImageFilter imageFileNames, int totalNumber, + void PairwiseDistanceAmongImages(const std::vector& imageFileNames, int totalNumber, vnl_matrix& distanceMatrix); - double calculateDistanceMSD(std::string& filename1, std::string& filename2); - - int ReadImage(std::string filename, InternalImageType::Pointer& image); - - void WriteImageUCHAR(std::string filename, InternalImageType::Pointer image); + double calculateDistanceMSD(const std::string& filename1, const std::string& filename2); - void WriteImageINT(std::string filename, InternalImageType::Pointer image); + int ReadImage(const std::string filename, InternalImageType::Pointer& image); - void WriteImageSHORT(std::string filename, InternalImageType::Pointer image); + itk::ImageIOBase::IOComponentType GetIOPixelType(const std::string& filename); - void WriteImageFLOAT(std::string filename, InternalImageType::Pointer image); + /** Write the image to disk, after casting it the the provided pixel type. */ + template + int WriteImage(const std::string filename, InternalImageType::Pointer image); - void WriteImage(std::string filename, InternalImageType::Pointer image); + int WriteImage(const std::string filename, InternalImageType::Pointer image, itk::ImageIOBase::IOComponentType ioType); - void WriteImage(std::string filename, InternalImageType::Pointer image, char* outputType); + void WriteImage(const std::string filename, InternalImageType::Pointer image, char* outputType); private: - MABMISImageOperationFilter(const Self &); // purposely not implemented - void operator=(const Self &); // purposely not implemented + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISImageOperationFilter); int m_Root; protected: MABMISImageOperationFilter(); ~MABMISImageOperationFilter(); - void PrintSelf(std::ostream& os, Indent indent) const; + void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE; }; } // namespace itk } // namespace Statistics diff --git a/itkMABMISImageOperationFilter.hxx b/itkMABMISImageOperationFilter.hxx index d86e2d0..cd0c4a7 100644 --- a/itkMABMISImageOperationFilter.hxx +++ b/itkMABMISImageOperationFilter.hxx @@ -1,3 +1,4 @@ +#include "itkMABMISImageOperationFilter.h" #ifndef __itkMABMISImageOperationFilter_hxx #define __itkMABMISImageOperationFilter_hxx @@ -30,12 +31,10 @@ MABMISImageOperationFilter template int MABMISImageOperationFilter -::ReadImage(std::string filename, InternalImageType::Pointer& image) +::ReadImage(const std::string filename, InternalImageType::Pointer& image) { // std::cout << "Reading "<< filename << std:: endl; - InternalImageReaderType::Pointer internalImageReader = InternalImageReaderType::New(); - internalImageReader->SetFileName(filename); try { @@ -50,159 +49,140 @@ MABMISImageOperationFilter return 0; } -template -void +template +itk::ImageIOBase::IOComponentType MABMISImageOperationFilter -::WriteImageUCHAR(std::string filename, InternalImageType::Pointer image) +::GetIOPixelType(const std::string & filename) { - Internal2CharCastFilterType::Pointer caster = Internal2CharCastFilterType::New(); + itk::ImageIOBase::Pointer imageIO; - caster->SetInput(image); - CharImageWriterType::Pointer writer = CharImageWriterType::New(); - writer->SetFileName(filename); - writer->SetInput(caster->GetOutput() ); - writer->SetUseCompression( false ); try { - writer->Update(); + imageIO = itk::ImageIOFactory::CreateImageIO(filename.c_str(), itk::ImageIOFactory::ReadMode); + if( imageIO ) + { + imageIO->SetFileName(filename); + imageIO->ReadImageInformation(); + return imageIO->GetComponentType(); + } + else + { + std::cout << "Could not create the imageIO for file " << filename << "." << std::endl; + exit( EXIT_FAILURE ); + } } - catch( itk::ExceptionObject & err ) + catch( itk::ExceptionObject& err ) { - std::cerr << err << std::endl; - return; + std::cout << "Could not read the image information of " << filename << "." << std::endl; + std::cout << err << std::endl; + exit( EXIT_FAILURE ); } - return; + + return itk::ImageIOBase::UCHAR; //this should never be reached } template -void +template +int MABMISImageOperationFilter -::WriteImageINT(std::string filename, InternalImageType::Pointer image) +::WriteImage(const std::string filename, InternalImageType::Pointer image) { - Internal2IntCastFilterType::Pointer caster = Internal2IntCastFilterType::New(); - + typedef itk::Image GivenImageType; + typedef itk::CastImageFilter CasterType; + typename CasterType::Pointer caster = CasterType::New(); + typedef itk::ImageFileWriter WriterType; + typename WriterType::Pointer writer = WriterType::New(); caster->SetInput(image); - IntImageWriterType::Pointer writer = IntImageWriterType::New(); + writer->SetInput(caster->GetOutput()); writer->SetFileName(filename); - writer->SetInput(caster->GetOutput() ); - writer->SetUseCompression( false ); + writer->SetUseCompression(true); try { writer->Update(); } - catch( itk::ExceptionObject & err ) + catch (itk::ExceptionObject & err) { std::cerr << err << std::endl; - return; + return -1; } - return; + return 0; } -template -void +template +int MABMISImageOperationFilter -::WriteImageSHORT(std::string filename, InternalImageType::Pointer image) +::WriteImage(const std::string filename, InternalImageType::Pointer image, itk::ImageIOBase::IOComponentType ioType) { - Internal2ShortCastFilterType::Pointer caster = Internal2ShortCastFilterType::New(); - - caster->SetInput(image); - ShortImageWriterType::Pointer writer = ShortImageWriterType::New(); - writer->SetFileName(filename); - writer->SetInput(caster->GetOutput() ); - writer->SetUseCompression( false ); - try + if( ioType == itk::ImageIOBase::UCHAR ) { - writer->Update(); + return this->WriteImage(filename, image); } - catch( itk::ExceptionObject & err ) + else if( ioType == itk::ImageIOBase::CHAR ) { - std::cerr << err << std::endl; - return; + return this->WriteImage(filename, image); } - return; -} - -template -void -MABMISImageOperationFilter -::WriteImageFLOAT(std::string filename, InternalImageType::Pointer image) -{ - Internal2FloatCastFilterType::Pointer caster = Internal2FloatCastFilterType::New(); - - caster->SetInput(image); - FloatImageWriterType::Pointer writer = FloatImageWriterType::New(); - writer->SetFileName(filename); - writer->SetInput(caster->GetOutput() ); - writer->SetUseCompression( false ); - try + else if( ioType == itk::ImageIOBase::USHORT ) { - writer->Update(); + return this->WriteImage(filename, image); } - catch( itk::ExceptionObject & err ) + else if( ioType == itk::ImageIOBase::SHORT ) { - std::cerr << err << std::endl; - return; + return this->WriteImage(filename, image); } - return; -} - -template -void -MABMISImageOperationFilter -::WriteImage(std::string filename, InternalImageType::Pointer image) -{ - Internal2CharCastFilterType::Pointer caster = Internal2CharCastFilterType::New(); - - caster->SetInput(image); - CharImageWriterType::Pointer writer = CharImageWriterType::New(); - writer->SetFileName(filename); - writer->SetInput(caster->GetOutput() ); - writer->SetUseCompression( false ); - try + else if( ioType == itk::ImageIOBase::UINT ) { - writer->Update(); + return this->WriteImage(filename, image); } - catch( itk::ExceptionObject & err ) + else if( ioType == itk::ImageIOBase::INT ) { - std::cerr << err << std::endl; - return; + return this->WriteImage(filename, image); } - return; -} - -template -void -MABMISImageOperationFilter -::WriteImage(std::string filename, InternalImageType::Pointer image, char* outputType) -{ - // - if( strcmp(outputType, "uchar") == 0 ) + else if( ioType == itk::ImageIOBase::ULONG ) { - WriteImageUCHAR(filename, image); + return this->WriteImage(filename, image); } - else if( strcmp(outputType, "short") == 0 ) + else if( ioType == itk::ImageIOBase::LONG ) { - WriteImageSHORT(filename, image); + return this->WriteImage(filename, image); } - else if( strcmp(outputType, "int") == 0 ) + //else if( ioType == itk::ImageIOBase::ULONGLONG ) + // { + // return this->WriteImage(filename, image); + // } + //else if( ioType == itk::ImageIOBase::LONGLONG ) + // { + // return this->WriteImage(filename, image); + // } + else if(ioType == itk::ImageIOBase::FLOAT) { - WriteImageINT(filename, image); + return this->WriteImage(filename, image); } - else if( strcmp(outputType, "float") == 0 ) + else if(ioType == itk::ImageIOBase::DOUBLE) { - WriteImageFLOAT(filename, image); + return this->WriteImage(filename, image); } - else // default + else { - WriteImage(filename, image); + std::cerr << itk::ImageIOBase::GetComponentTypeAsString(ioType) << " not supported as I/O pixel type" << std::endl; + return -1; } + return -2; +} + +template +void +MABMISImageOperationFilter +::WriteImage(const std::string filename, InternalImageType::Pointer image, char* outputType) +{ + itk::ImageIOBase::IOComponentType ioType = itk::ImageIOBase::GetComponentTypeFromString(outputType); + return this->WriteImage(filename, image, ioType); } // calculate distance between two images by mean squared difference template double MABMISImageOperationFilter -::calculateDistanceMSD(std::string& imageName1, std::string& imageName2) +::calculateDistanceMSD(const std::string& imageName1, const std::string& imageName2) { double dist; @@ -245,7 +225,7 @@ MABMISImageOperationFilter template void MABMISImageOperationFilter -::PairwiseDistanceAmongImages(std::vector imageFileNames, int totalNumber, vnl_matrix& distanceMatrix) +::PairwiseDistanceAmongImages(const std::vector& imageFileNames, int totalNumber, vnl_matrix& distanceMatrix) { for( int i = 0; i < totalNumber; ++i ) diff --git a/itkMABMISImageRegistrationFilter.h b/itkMABMISImageRegistrationFilter.h index 0e14c97..1fe5940 100644 --- a/itkMABMISImageRegistrationFilter.h +++ b/itkMABMISImageRegistrationFilter.h @@ -1,22 +1,12 @@ #ifndef __itkMABMISImageRegistrationFilter_h #define __itkMABMISImageRegistrationFilter_h -#include -#include - -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" -#include "itkImageRegionIterator.h" +#include "commonMABMIS.h" #include "itkHistogramMatchingImageFilter.h" - -#include "itkImage.h" - #include "itkMABMISDeformationFieldFilter.h" #include "itkMABMISImageOperationFilter.h" -#define ImageDimension 3 - namespace itk { namespace Statistics @@ -104,14 +94,13 @@ class MABMISImageRegistrationFilter : public ImageToImageFilter // float sigmaDef = 1.5; if( sigmaDef > 0.1 ) { - filter->SmoothDeformationFieldOn(); + filter->SmoothDisplacementFieldOn(); filter->SetStandardDeviations( sigmaDef ); } else { - filter->SmoothDeformationFieldOff(); + filter->SmoothDisplacementFieldOff(); } // set up smoothing kernel for update field @@ -112,16 +112,15 @@ MABMISImageRegistrationFilter InternalImageType, InternalImageType, DeformationFieldType, InternalPixelType> MultiResRegistrationFilterType; MultiResRegistrationFilterType::Pointer multires = MultiResRegistrationFilterType::New(); - typedef itk::VectorLinearInterpolateNearestNeighborExtrapolateImageFunction FieldInterpolatorType; + typedef itk::LinearInterpolateImageFunction FieldInterpolatorType; + typedef itk::NearestNeighborExtrapolateImageFunction FieldExtrapolatorType; FieldInterpolatorType::Pointer VectorInterpolator = FieldInterpolatorType::New(); + FieldExtrapolatorType::Pointer VectorExtrapolator = FieldExtrapolatorType::New(); -#if ( ITK_VERSION_MAJOR > 3 ) || ( ITK_VERSION_MAJOR == 3 && ITK_VERSION_MINOR > 8 ) multires->GetFieldExpander()->SetInterpolator(VectorInterpolator); -#endif + multires->GetFieldExpander()->SetExtrapolator(VectorExtrapolator); std::vector curNumIterations; - // unsigned int curNumOfIterations[] = {15,10,5}; for( int i = 0; i < res; ++i ) { curNumIterations.push_back(iterInResolutions[i]); @@ -191,7 +190,7 @@ MABMISImageRegistrationFilter { int res = iterInResolutions.size(); // read initial deformation field file - DeformationFieldType::Pointer initDeformationField = 0; + DeformationFieldType::Pointer initDeformationField = nullptr; dfoperator->ReadDeformationField(initDeformationFieldFileName, initDeformationField); @@ -236,12 +235,12 @@ MABMISImageRegistrationFilter // float sigmaDef = 1.5; if( sigmaDef > 0.1 ) { - filter->SmoothDeformationFieldOn(); + filter->SmoothDisplacementFieldOn(); filter->SetStandardDeviations( sigmaDef ); } else { - filter->SmoothDeformationFieldOff(); + filter->SmoothDisplacementFieldOff(); } // set up smoothing kernel for update field @@ -263,14 +262,14 @@ MABMISImageRegistrationFilter InternalImageType, InternalImageType, DeformationFieldType, InternalPixelType> MultiResRegistrationFilterType; MultiResRegistrationFilterType::Pointer multires = MultiResRegistrationFilterType::New(); - typedef itk::VectorLinearInterpolateNearestNeighborExtrapolateImageFunction FieldInterpolatorType; + typedef itk::LinearInterpolateImageFunction FieldInterpolatorType; + typedef itk::NearestNeighborExtrapolateImageFunction FieldExtrapolatorType; FieldInterpolatorType::Pointer VectorInterpolator = FieldInterpolatorType::New(); + FieldExtrapolatorType::Pointer VectorExtrapolator = FieldExtrapolatorType::New(); -#if ( ITK_VERSION_MAJOR > 3 ) || ( ITK_VERSION_MAJOR == 3 && ITK_VERSION_MINOR > 8 ) multires->GetFieldExpander()->SetInterpolator(VectorInterpolator); -#endif + multires->GetFieldExpander()->SetExtrapolator(VectorExtrapolator); std::vector curNumIterations; // unsigned int curNumOfIterations[] = {15,10,5}; for( int i = 0; i < res; ++i ) @@ -292,7 +291,7 @@ MABMISImageRegistrationFilter } // set initial - multires->SetArbitraryInitialDeformationField( initDeformationField ); + multires->SetArbitraryInitialDisplacementField( initDeformationField ); // Compute the deformation field try diff --git a/itkMABMISSimulateData.h b/itkMABMISSimulateData.h index 94023be..9e03826 100644 --- a/itkMABMISSimulateData.h +++ b/itkMABMISSimulateData.h @@ -1,21 +1,12 @@ #ifndef __itkMABMISSimulateData_h #define __itkMABMISSimulateData_h -#include -#include - -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" -#include "itkImageRegionIterator.h" - -#include "itkImage.h" +#include "commonMABMIS.h" #include "itkMABMISDeformationFieldFilter.h" #include "itkMABMISImageOperationFilter.h" #include "itkMABMISBasicOperationFilter.h" -#define ImageDimension 3 - namespace itk { namespace Statistics @@ -71,12 +62,12 @@ class MABMISSimulateData : public ImageToImageFilter typename ImageOperationType::Pointer imgoperator; typename BasicOperationType::Pointer basicoperator; - int DoPCATraining(std::vector deformationFieldFileNames, int numFiles, - std::vector allImgFileName, int root); + int DoPCATraining(const std::vector& deformationFieldFileNames, int numFiles, + const std::vector& allImgFileName, int root); - void LoadIntoArray(std::string resampledDeformationFieldFileName, float* df_vector); + void LoadIntoArray(const std::string& resampledDeformationFieldFileName, float* df_vector); - void SaveFromArray(std::string deformationFieldFileName, float* df_vector, int sx, int sy, int sz); + void SaveFromArray(const std::string& deformationFieldFileName, float* df_vector, int sx, int sy, int sz); itkSetMacro(Root, int); itkSetMacro(Imx, int); @@ -85,8 +76,7 @@ class MABMISSimulateData : public ImageToImageFilter itkSetMacro(AtlasSize, int); itkSetMacro(SimulateSize, int); private: - MABMISSimulateData(const Self &); // purposely not implemented - void operator=(const Self &); // purposely not implemented + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISSimulateData); int m_Root; @@ -99,7 +89,7 @@ class MABMISSimulateData : public ImageToImageFilter protected: MABMISSimulateData(); ~MABMISSimulateData(); - void PrintSelf(std::ostream& os, Indent indent) const; + void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE; }; } // namespace itk } // namespace Statistics diff --git a/itkMABMISSimulateData.hxx b/itkMABMISSimulateData.hxx index 710852b..7c993f1 100644 --- a/itkMABMISSimulateData.hxx +++ b/itkMABMISSimulateData.hxx @@ -3,11 +3,7 @@ #include "itkMABMISSimulateData.h" -//For definition of FILESEP -#include "itkMABMISAtlasXMLFile.h" - -int numEigenVector = 4; // t -int numSampleEachDirection = 4; // n n^t total number of intermediate templates +constexpr int numSampleEachDirection = 4; // n n^numEigenVector total number of intermediate templates namespace itk { @@ -39,8 +35,8 @@ MABMISSimulateData template int MABMISSimulateData -::DoPCATraining(std::vector deformationFieldFileNames, int numFiles, - std::vector allImgFileName, int root) +::DoPCATraining(const std::vector& deformationFieldFileNames, int numFiles, + const std::vector& allImgFileName, int root) { bool doPCATraining = true; int sampleRate = 2; @@ -48,8 +44,9 @@ MABMISSimulateData int size_x = 0; int size_y = 0; int size_z = 0; int size_xn = 0; int size_yn = 0; int size_zn = 0; int size_dfn = 0; // size of sub-sampled deformation field + int numEigenVector = 4; - float c4[] = {-0.8416, -0.2533, 0.2533, 0.8416}; + float c4[] = {-0.8416f, -0.2533f, 0.2533f, 0.8416f}; float* c = NULL; @@ -73,7 +70,7 @@ MABMISSimulateData /////////////////////////////////////// // initialize - DeformationFieldType::Pointer curDeformationField = 0; + DeformationFieldType::Pointer curDeformationField = nullptr; dfoperator->ReadDeformationField(deformationFieldFileNames[0], curDeformationField); @@ -92,7 +89,7 @@ MABMISSimulateData // do PCA training //////////////////////////////////////// // std::cerr << "Resample deformation fields ..." << std::endl; - if (numEigenVector > numFiles) numEigenVector = numFiles; + if (numEigenVector > numFiles) numEigenVector = numFiles; vnl_matrix df_eigenvector(size_dfn, numEigenVector); vnl_vector df_eigenvalues(numEigenVector); //////////////////////////////////////// @@ -171,13 +168,13 @@ MABMISSimulateData // for debugging // cerr << "Pass: PCA" << std::endl; - + // make a temporary folder to store the intermediate files std::string tempFolder = "temp_PCATraining"; char numStr[10]; int num = rand() % 10000; - num = rand() % 10000; + num = rand() % 10000; sprintf(numStr, "%04d", num); tempFolder = tempFolder + numStr; @@ -191,19 +188,20 @@ MABMISSimulateData outputFolder = allImgFileName[0].substr(0, sep); } - if( 1 ) // if (0) // if (1) + if( 1 ) // if (0) { std::vector coeff(numEigenVector,0); const int numAllCombinations = (int)pow( (float)numSampleEachDirection, numEigenVector); // template image: the root of the tree - InternalImageType::Pointer templateImage = 0; + InternalImageType::Pointer templateImage = nullptr; if( imgoperator->ReadImage(allImgFileName[root], templateImage) != 0 ) { std::cerr << " Cannot read image: " << allImgFileName[root] << std::endl; std::cerr << "Please verify the file exists. " << std::endl; return -1; } + const itk::ImageIOBase::IOComponentType ioType = imgoperator->GetIOPixelType(allImgFileName[root]); std::cout << "Generate intermediate deformations from PCA results... " << std::endl; for( int i = 0; i < numAllCombinations; ++i ) @@ -232,7 +230,7 @@ MABMISSimulateData df_intermediate_sub += c[coeff[j]] * df_eigenvector.get_column(j); } - std::string index_string; basicoperator->myitoa( i, index_string, 3 ); + std::string index_string = basicoperator->myitoa(i, 3); std::string intermediateSubDeformationFieldFileName = "inter_deform_sub_000.nii.gz"; intermediateSubDeformationFieldFileName.erase( intermediateSubDeformationFieldFileName.end() - 7, intermediateSubDeformationFieldFileName.end() ); @@ -277,7 +275,7 @@ MABMISSimulateData intermediateReversedDeformationFieldFileName = tempFolder + FILESEP + intermediateReversedDeformationFieldFileName; - DeformationFieldType::Pointer deformationField = 0; + DeformationFieldType::Pointer deformationField = nullptr; if( dfoperator->ReadDeformationField(intermediateDeformationFieldFileName, deformationField) != 0 ) { std::cerr << " Cannot read deformation field: " << intermediateDeformationFieldFileName << std::endl; @@ -295,10 +293,10 @@ MABMISSimulateData dfoperator->WriteDeformationField(intermediateReversedDeformationFieldFileName, reversedDeformationField); // apply intermediate deformation field to get intermediate template - InternalImageType::Pointer intermediateTemplateImage = 0; + InternalImageType::Pointer intermediateTemplateImage = nullptr; dfoperator->ApplyDeformationField(templateImage, reversedDeformationField, intermediateTemplateImage, true); // write image - imgoperator->WriteImage(intermediateTemplateFileName, intermediateTemplateImage); + imgoperator->WriteImage(intermediateTemplateFileName, intermediateTemplateImage, ioType); // intermediateFileNamesFile << intermediateTemplateFileName << std::endl; } // end for numOfAllCombinations @@ -322,7 +320,7 @@ MABMISSimulateData { for( int j = 0; j < numAllCombinations; ++j ) { - std::string index_string; basicoperator->myitoa( j, index_string, 3 ); + std::string index_string = basicoperator->myitoa( j, 3 ); std::string curInterTempFileName = "inter_template_000.nii.gz"; curInterTempFileName.erase(curInterTempFileName.end() - 7, curInterTempFileName.end() ); curInterTempFileName.append(index_string); @@ -359,14 +357,14 @@ MABMISSimulateData { int index_inter = index[i + (numAllCombinations - m_SimulateSize)]; - std::string index_string2; basicoperator->myitoa( index_inter, index_string2, 3 ); + std::string index_string2 = basicoperator->myitoa( index_inter, 3 ); std::string curInterTempFileName = "inter_deform_000.nii.gz"; curInterTempFileName.erase(curInterTempFileName.end() - 7, curInterTempFileName.end() ); curInterTempFileName.append(index_string2); curInterTempFileName.append(".nii.gz"); curInterTempFileName = tempFolder + FILESEP + curInterTempFileName; - std::string index_string3; basicoperator->myitoa( i, index_string3, 3 ); + std::string index_string3 = basicoperator->myitoa( i, 3 ); std::string outDFName = "simulated_deform_" + index_string3 + ".nii.gz"; outDFName = outputFolder + FILESEP + outDFName; @@ -384,7 +382,7 @@ MABMISSimulateData // remove irrelevant files for( int i = 0; i < numAllCombinations; ++i ) { - std::string index_string; basicoperator->myitoa( i, index_string, 3 ); + std::string index_string = basicoperator->myitoa( i, 3 ); std::string curInterTempFileName = "inter_template_" + index_string + ".nii.gz"; std::string curInterTempDeformFileName = "inter_deform_" + index_string + ".nii.gz"; @@ -416,7 +414,7 @@ MABMISSimulateData template void MABMISSimulateData -::SaveFromArray(std::string deformationFieldFileName, float* df_vector, int sx, int sy, int sz) +::SaveFromArray(const std::string& deformationFieldFileName, float* df_vector, int sx, int sy, int sz) { std::cerr << "deformationFieldFileName: " << deformationFieldFileName << std::endl; @@ -462,15 +460,13 @@ MABMISSimulateData template void MABMISSimulateData -::LoadIntoArray(std::string resampledDeformationFieldFileName, float* df_vector) +::LoadIntoArray(const std::string& resampledDeformationFieldFileName, float* df_vector) { // read deformation field and load it into array - DeformationFieldType::Pointer dfImage = 0; + DeformationFieldType::Pointer dfImage = nullptr; dfoperator->ReadDeformationField(resampledDeformationFieldFileName, dfImage); - DeformationFieldType::SizeType im_size = dfImage->GetLargestPossibleRegion().GetSize(); - // load original image DeformationFieldIteratorType itOrigin(dfImage, dfImage->GetLargestPossibleRegion() ); int index = 0; diff --git a/itkMABMISTreeOperation.h b/itkMABMISTreeOperation.h index 4110036..118b74e 100644 --- a/itkMABMISTreeOperation.h +++ b/itkMABMISTreeOperation.h @@ -1,21 +1,14 @@ #ifndef __itkMABMISTreeOperation_h #define __itkMABMISTreeOperation_h -#include #include - -#include "itkImageFileReader.h" -#include "itkImageFileWriter.h" #include "itkImageRegionIterator.h" -#include - +#include "commonMABMIS.h" #include "itkMABMISDeformationFieldFilter.h" #include "itkMABMISImageOperationFilter.h" #include "itkMABMISBasicOperationFilter.h" -#define ImageDimension 3 - namespace itk { namespace Statistics @@ -100,8 +93,7 @@ class MABMISTreeOperation : public ImageToImageFilter itkSetMacro(SimulateSize, int); itkSetMacro(AllAtlasNumber, int); private: - MABMISTreeOperation(const Self &); // purposely not implemented - void operator=(const Self &); // purposely not implemented + ITK_DISALLOW_COPY_AND_ASSIGN(MABMISTreeOperation); int m_Root; @@ -117,7 +109,7 @@ class MABMISTreeOperation : public ImageToImageFilter protected: MABMISTreeOperation(); ~MABMISTreeOperation(); - void PrintSelf(std::ostream& os, Indent indent) const; + void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE; }; } // namespace itk } // namespace Statistics diff --git a/itkMABMISTreeOperation.hxx b/itkMABMISTreeOperation.hxx index dd36db2..26e31f0 100644 --- a/itkMABMISTreeOperation.hxx +++ b/itkMABMISTreeOperation.hxx @@ -2,6 +2,7 @@ #define __itkMABMISTreeOperation_hxx #include "itkMABMISTreeOperation.h" +#include namespace itk { @@ -222,7 +223,7 @@ MABMISTreeOperation ///////////////////////////////////////////// // rebuild the distance matrix based on connectivity on MST vnl_matrix curDistanceTemp(nnode, nnode); - static const double big_number = vcl_sqrt(itk::NumericTraits< double >::max()); + static const double big_number = std::sqrt(itk::NumericTraits< double >::max()); for( int i = 0; i < nnode; ++i ) { for( int j = 0; j < nnode; ++j ) @@ -382,7 +383,7 @@ MABMISTreeOperation // rebuild the distance matrix based on connectivity on MST double* * curDistanceTemp = new double *[nnode]; - static const double big_number = vcl_sqrt(itk::NumericTraits< double >::max()); + static const double big_number = std::sqrt(itk::NumericTraits< double >::max()); for( int i = 0; i < nnode; ++i ) { curDistanceTemp[i] = new double[nnode];