diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp index ef672c393f..7880f8c5bf 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp @@ -1,217 +1,246 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include +#include using namespace std; typedef itksys::SystemTools ist; typedef itk::Image ItkUcharImgType; typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType; void CreateFolderStructure(const std::string& path) { if (ist::PathExists(path)) ist::RemoveADirectory(path); ist::MakeDirectory(path); ist::MakeDirectory(path + "/known_tracts/"); ist::MakeDirectory(path + "/candidate_tracts/"); ist::MakeDirectory(path + "/implausible_tracts/"); ist::MakeDirectory(path + "/skipped_masks/"); } ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); mitk::CastToItkImage(img, itkMask); return itkMask; } std::vector< MaskType > get_file_list(const std::string& path, int dropout, const std::string& skipped_path) { - std::srand(std::time(0)); + std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds >(std::chrono::system_clock::now().time_since_epoch()); + std::srand(ms.count()); std::vector< MaskType > mask_list; itk::Directory::Pointer dir = itk::Directory::New(); + int skipped = 0; if (dir->Load(path.c_str())) { int n = dir->GetNumberOfFiles(); + int num_images = 0; + + std::vector< int > im_indices; for (int r = 0; r < n; r++) { const char *filename = dir->GetFile(r); + std::string ext = ist::GetFilenameExtension(filename); + if (ext==".nii" || ext==".nii.gz" || ext==".nrrd") + { + ++num_images; + im_indices.push_back(r); + } + } + + int skipping_num = 0; + if (dropout>1) + skipping_num = (float)num_images/dropout; + if (dropout>=num_images) + skipping_num = (float)num_images/2.0; + + std::random_shuffle(im_indices.begin(), im_indices.end()); + MITK_INFO << num_images; + MITK_INFO << dropout; + MITK_INFO << "Skipping " << skipping_num << " images"; + + int c = -1; + for (int r : im_indices) + { + c++; + const char *filename = dir->GetFile(r); - if (dropout>1 && std::rand()%dropout==0 && ist::GetFilenameWithoutExtension(filename)!="CC") + if (c matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); matrix->Identity(); matrix->SetElement(0,0,-matrix->GetElement(0,0)); matrix->SetElement(1,1,-matrix->GetElement(1,1)); geometry->SetIndexToWorldTransformByVtkMatrix(matrix); vtkSmartPointer transformFilter = vtkSmartPointer::New(); transformFilter->SetInputData(fib->GetFiberPolyData()); transformFilter->SetTransform(geometry->GetVtkTransform()); transformFilter->Update(); mitk::FiberBundle::Pointer transformed_fib = mitk::FiberBundle::New(transformFilter->GetOutput()); return transformed_fib; } /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Tract Plausibility"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output folder", us::Any(), false); parser.addArgument("reference_mask_folder", "m", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", false); parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "", ""); parser.addArgument("dropout", "", mitkCommandLineParser::Int, "", "", false); + parser.addArgument("overlap", "", mitkCommandLineParser::Float, "", "", false, 0.7); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string fibFile = us::any_cast(parsedArgs["input"]); string reference_mask_folder = us::any_cast(parsedArgs["reference_mask_folder"]); string out_folder = us::any_cast(parsedArgs["out"]); string gray_matter_mask = ""; if (parsedArgs.count("gray_matter_mask")) gray_matter_mask = us::any_cast(parsedArgs["gray_matter_mask"]); int dropout = 1; if (parsedArgs.count("dropout")) dropout = us::any_cast(parsedArgs["dropout"]); + float overlap = 0.7; + if (parsedArgs.count("overlap")) + overlap = us::any_cast(parsedArgs["overlap"]); + try { CreateFolderStructure(out_folder); std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder, dropout, out_folder + "/skipped_masks/"); if (known_tract_masks.empty()) return EXIT_FAILURE; mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); if (gray_matter_mask.compare("")!=0) { ItkUcharImgType::Pointer gm_image = LoadItkMaskImage(gray_matter_mask); // filter gray matter fibers mitk::FiberBundle::Pointer not_gm_fibers = inputTractogram->ExtractFiberSubset(gm_image, false, true, true); mitk::IOUtil::Save(not_gm_fibers, out_folder + "/implausible_tracts/no_gm_endings.trk"); inputTractogram = inputTractogram->ExtractFiberSubset(gm_image, false, false, true); } // resample fibers float minSpacing = 1; if(std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[1] && std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[2]) minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]; else if (std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] < std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]) minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[1]; else minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]; inputTractogram->ResampleLinear(minSpacing/5); // find known tracts via overlap match mitk::FiberBundle::Pointer all_known_tracts = nullptr; for ( MaskType mask : known_tract_masks ) { ItkUcharImgType::Pointer mask_image = std::get<0>(mask); std::string mask_name = std::get<1>(mask); - mitk::FiberBundle::Pointer known_tract = inputTractogram->ExtractFiberSubset(mask_image, true, false, false, 0.7, false); + mitk::FiberBundle::Pointer known_tract = inputTractogram->ExtractFiberSubset(mask_image, true, false, false, overlap, false); mitk::IOUtil::Save(known_tract, out_folder + "/known_tracts/" + mask_name + ".trk"); if (all_known_tracts.IsNull()) all_known_tracts = mitk::FiberBundle::New(known_tract->GetFiberPolyData()); else all_known_tracts = all_known_tracts->AddBundle(known_tract); } mitk::IOUtil::Save(all_known_tracts, out_folder + "/known_tracts/all_known_tracts.trk"); mitk::IOUtil::Save(TransformToMRtrixSpace(all_known_tracts), out_folder + "/known_tracts/all_known_tracts_mrtrixspace.fib"); mitk::FiberBundle::Pointer remaining_tracts = inputTractogram->SubtractBundle(all_known_tracts); mitk::IOUtil::Save(remaining_tracts, out_folder + "/candidate_tracts/remaining_tracts.trk"); mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), out_folder + "/candidate_tracts/remaining_tracts_mrtrixspace.fib"); - MITK_INFO << "step_size: " << minSpacing/5; } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; }