Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
TumorProgressionMapping.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 // MITK
13 #include <mitkIOUtil.h>
14 #include <mitkImageAccessByItk.h>
15 #include <mitkImageCast.h>
18 #include <mitkImageToItk.h>
19 #include <mitkPixelType.h>
20 #include <mitkPlaneGeometry.h>
21 
22 // ITK
23 #include "itkImageDuplicator.h"
24 #include "itkPermuteAxesImageFilter.h"
25 #include "itkTileImageFilter.h"
26 #include <itkExtractImageFilter.h>
27 #include <itkImageFileWriter.h>
28 #include <itkRGBPixel.h>
29 #include <itkResampleImageFilter.h>
30 #include <itkVTKImageImport.h>
31 
32 #include <itkDirectory.h>
33 #include <itksys/SystemTools.hxx>
34 // VTK
35 #include <vtkCellArray.h>
36 #include <vtkFreeTypeStringToImage.h>
37 #include <vtkImageData.h>
38 #include <vtkTextProperty.h>
39 #include <vtkUnicodeString.h>
40 
41 // MITK
42 #include <mitkCommandLineParser.h>
43 
44 #include "itkRescaleIntensityImageFilter.h"
45 
46 #include <typeExtension.h>
47 
48 using namespace std;
49 
50 typedef unsigned short PixelType;
51 typedef double InputPixelType;
52 typedef itk::Image<itk::RGBPixel<PixelType>, 2> RGB2DImage;
53 typedef itk::RGBPixel<PixelType> RGBPixelType;
54 
55 typedef std::vector<std::string> FileList;
56 typedef std::vector<mitk::Image::Pointer> ImageList;
57 
59 static FileList CreateFileList(std::string folder, std::string postfix)
60 {
61  itk::Directory::Pointer dir = itk::Directory::New();
62  FileList fileList;
63 
64  if (dir->Load(folder.c_str()))
65  {
66  int n = dir->GetNumberOfFiles();
67  for (int r = 0; r < n; r++)
68  {
69  std::string filename = dir->GetFile(r);
70  if (filename == "." || filename == "..")
71  continue;
72  filename = folder + filename;
73  if (!itksys::SystemTools::FileExists(filename.c_str()))
74  continue;
75 
76  if (postfix.compare(filename.substr(filename.length() - postfix.length())) == 0)
77  fileList.push_back(filename);
78  }
79  }
80  std::sort(fileList.begin(), fileList.end());
81  return fileList;
82 }
83 
84 class mitkProgressionVisualization
85 {
86 public:
87  mitkProgressionVisualization() {}
93  mitk::Image::Pointer ConvertToRGBImage(mitk::Image::Pointer grayImage)
94  {
96  unsigned int *dim = grayImage->GetDimensions();
97  rgbImage->Initialize(mitk::MakePixelType<PixelType, RGBPixelType, 3>(), 3, dim);
98  rgbImage->SetGeometry(grayImage->GetGeometry());
99 
100  itk::Image<InputPixelType, 3>::Pointer itkGrayImage = itk::Image<InputPixelType, 3>::New();
101  mitk::CastToItkImage(grayImage, itkGrayImage);
102 
104 
105  typedef itk::RescaleIntensityImageFilter<itk::Image<InputPixelType, 3>, itk::Image<PixelType, 3>> RescaleFilterType;
106  RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
107  rescaleFilter->SetInput(itkGrayImage);
108  rescaleFilter->SetOutputMinimum(0);
109  rescaleFilter->SetOutputMaximum(255 * 255);
110  rescaleFilter->Update();
111 
112  itk::Index<3> idx;
113  RGBPixelType value;
114 
115  // Fill rgb image with gray values
116  for (idx[2] = 0; (unsigned int)idx[2] < dim[2]; idx[2]++)
117  {
118  for (idx[1] = 0; (unsigned int)idx[1] < dim[1]; idx[1]++)
119  {
120  for (idx[0] = 0; (unsigned int)idx[0] < dim[0]; idx[0]++)
121  {
122  value.Fill(rescaleFilter->GetOutput()->GetPixel(idx));
123  writeAcc.SetPixelByIndex(idx, value);
124  }
125  }
126  }
127  return rgbImage;
128  }
129 
130  RGB2DImage::Pointer TextAsImage(std::string text)
131  {
132  vtkSmartPointer<vtkImageData> textImage = vtkSmartPointer<vtkImageData>::New();
133  vtkSmartPointer<vtkFreeTypeStringToImage> freetype = vtkSmartPointer<vtkFreeTypeStringToImage>::New();
134  vtkSmartPointer<vtkTextProperty> prop = vtkSmartPointer<vtkTextProperty>::New();
135  float color[3] = {1, 1, 1};
136  float opacity = 1.0;
137  prop->SetColor(color[0], color[1], color[2]);
138  prop->SetFontSize(40);
139  prop->SetOpacity(opacity);
140  textImage->AllocateScalars(VTK_UNSIGNED_SHORT, 3);
141  freetype->RenderString(prop, vtkUnicodeString::from_utf8(text.c_str()), 72, textImage);
142  textImage->Modified();
143 
144  int *extent = textImage->GetExtent();
145 
146  RGB2DImage::Pointer itkImage = RGB2DImage::New();
147 
148  RGB2DImage::IndexType start;
149  start.Fill(0);
150  RGB2DImage::SizeType size;
151  size[0] = extent[1];
152  size[1] = extent[3];
153  size[2] = 0;
154 
155  RGB2DImage::RegionType region(start, size);
156  itkImage->SetRegions(region);
157  itkImage->Allocate();
158  RGB2DImage::IndexType idx;
159 
160  for (unsigned int y = 0; y < size[1]; y++)
161  {
162  for (unsigned int x = 0; x < size[0]; x++)
163  {
164  PixelType *pixel = static_cast<PixelType *>(textImage->GetScalarPointer(x, y, 0));
165  RGBPixelType values;
166  values.Fill(0);
167  values[0] = pixel[1];
168  values[1] = pixel[1];
169  values[2] = pixel[1];
170  idx.Fill(0);
171  idx[0] = x;
172  idx[1] = y;
173 
174  itkImage->SetPixel(idx, values);
175  }
176  }
177 
178  typedef itk::PermuteAxesImageFilter<RGB2DImage> PermuteAxesImageFilterType;
179 
180  itk::FixedArray<unsigned int, 2> order;
181  order[0] = 1;
182  order[1] = 0;
183 
184  PermuteAxesImageFilterType::Pointer permuteAxesFilter = PermuteAxesImageFilterType::New();
185  permuteAxesFilter->SetInput(itkImage);
186  permuteAxesFilter->SetOrder(order);
187  permuteAxesFilter->Update();
188 
189  RGB2DImage::Pointer itkResultImage = RGB2DImage::New();
190  itkResultImage->SetRegions(region);
191  itkResultImage->Allocate();
192 
193  itkResultImage->Graft(permuteAxesFilter->GetOutput());
194 
195  return itkResultImage;
196  }
197 
208  void AddColouredOverlay(mitk::Image::Pointer rgbImage, mitk::Image::Pointer overlayImage, mitk::Color color)
209  {
210  unsigned int *dim = rgbImage->GetDimensions();
211  itk::Image<PixelType, 3>::Pointer itkOverlayImage = itk::Image<PixelType, 3>::New();
212  mitk::CastToItkImage(overlayImage.GetPointer(), itkOverlayImage);
214 
215  itk::Index<3> idx;
216  itk::RGBPixel<PixelType> value;
217  unsigned short overlayVal = 0;
218  // Fill rgb image with gray values
219  for (idx[2] = 0; (unsigned int)idx[2] < dim[2]; idx[2]++)
220  {
221  for (idx[1] = 0; (unsigned int)idx[1] < dim[1]; idx[1]++)
222  {
223  for (idx[0] = 0; (unsigned int)idx[0] < dim[0]; idx[0]++)
224  {
225  overlayVal = 255 * itkOverlayImage->GetPixel(idx);
226  value = writeAcc.GetPixelByIndex(idx);
227  value[0] = std::min((int)(value[0] + overlayVal * color[0]), 254 * 255);
228  value[1] = std::min((int)(value[1] + overlayVal * color[1]), 254 * 255);
229  value[2] = std::min((int)(value[2] + overlayVal * color[2]), 254 * 255);
230  writeAcc.SetPixelByIndex(idx, value);
231  }
232  }
233  }
234  }
235 
236  itk::Image<itk::RGBPixel<PixelType>, 2>::Pointer ExtractSlice(
237  itk::Image<itk::RGBPixel<PixelType>, 3>::Pointer itkImage, unsigned int sliceNo)
238  {
239  typedef itk::Image<RGBPixelType, 3> InputImageType;
240  typedef itk::Image<RGBPixelType, 2> OutputImageType;
241 
242  int dim[3];
243  dim[0] = itkImage->GetLargestPossibleRegion().GetSize()[0];
244  dim[1] = itkImage->GetLargestPossibleRegion().GetSize()[1];
245  dim[2] = itkImage->GetLargestPossibleRegion().GetSize()[2];
246 
247  itk::Index<3> desiredStart;
248  itk::Size<3> desiredSize;
249 
250  // case AXIAL: // 3rd dimension fixed
251  desiredStart.Fill(0);
252  desiredStart[2] = sliceNo;
253  desiredSize.Fill(0);
254  desiredSize[0] = dim[0];
255  desiredSize[1] = dim[1];
256  desiredSize[2] = 0;
257 
258  itk::ImageRegion<3> desiredRegion(desiredStart, desiredSize);
259 
260  // Extract slice
261  itk::ExtractImageFilter<InputImageType, OutputImageType>::Pointer extractSlice =
262  itk::ExtractImageFilter<InputImageType, OutputImageType>::New();
263  extractSlice->SetInput(itkImage);
264  extractSlice->SetExtractionRegion(desiredRegion);
265  extractSlice->SetDirectionCollapseToIdentity();
266  extractSlice->Update();
267 
268  return extractSlice->GetOutput();
269  }
270 };
271 
272 static std::string GetName(std::string fileName, std::string suffix = "_T2.nrrd")
273 {
274  fileName = itksys::SystemTools::GetFilenameName(fileName);
275  return fileName.substr(0, fileName.length() - suffix.length() - 11); // 10 = date length
276 }
277 
278 static std::string GetDate(std::string fileName, std::string suffix = "_T2.nrrd")
279 {
280  fileName = itksys::SystemTools::GetFilenameName(fileName);
281  fileName = fileName.substr(fileName.length() - suffix.length() - 10, 10); // 10 = date length
282  return fileName;
283 }
284 
286 {
288  for (unsigned int i = 0; i < files.size(); ++i)
289  {
290  images.push_back(mitk::IOUtil::Load<mitk::Image>(files.at(i)));
291  }
292  return images;
293 }
294 
295 int main(int argc, char *argv[])
296 {
297  // Parse Command-Line Arguments
298  mitkCommandLineParser parser;
299  parser.setArgumentPrefix("--", "-");
300 
301  parser.setTitle("Tumor Progression Mapping");
302  parser.setCategory("Preprocessing Tools");
303  parser.setContributor("German Cancer Research Center (DKFZ)");
304  parser.setDescription("Convert a set of co-registered and resampled follow-up images into a 2D png overview (and "
305  "optionally in a 4D NRRD Volume).\nNaming convecntion of files is "
306  "IDENTIFIER_YYYY-MM-DD_MODALITY.nrrd");
307 
308  parser.setArgumentPrefix("--", "-");
309  parser.addArgument("input", "i", mitkCommandLineParser::Directory, "Input folder containing all follow ups", "", us::Any(), false, false, false, mitkCommandLineParser::Input);
310  parser.addArgument("output", "o", mitkCommandLineParser::File, "Output file (PNG)", "", us::Any(), false, false, false, mitkCommandLineParser::Output);
311  parser.addArgument("blanked", "b", mitkCommandLineParser::Bool, "Only Display Morphology");
312  parser.addArgument("morphology", "m", mitkCommandLineParser::String, "Morphology postfix.", "_T2.nrrd");
313  parser.addArgument(
314  "segmentation", "s", mitkCommandLineParser::String, "Segmentation postfix. Default: _GTV.nrrd", "_GTV.nrrd");
315  parser.addArgument("4dVolume", "v", mitkCommandLineParser::File, "Filepath to merged 4d NRRD Volume.", "", us::Any(), false, false, false, mitkCommandLineParser::Output);
316  parser.addArgument(
317  "skip", "k", mitkCommandLineParser::Int, "Number of slices to be skipped from top and from button (Default 0)");
318  parser.addArgument(
319  "interval", "n", mitkCommandLineParser::Int, "1 - for all slices, 2 - every second, 3 - every third ...");
320  parser.addArgument("opacity", "c", mitkCommandLineParser::Float, "Opacity of overlay [0,1] invisible -> visible");
321 
322  map<string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
323 
324  if (parsedArgs.size() == 0)
325  return EXIT_SUCCESS;
326 
327  // Show a help message
328  if (parsedArgs.count("help") || parsedArgs.count("h"))
329  {
330  std::cout << parser.helpText();
331  return EXIT_SUCCESS;
332  }
333 
334  std::string outputFile;
335  std::string inputFolder;
336 
337  if (parsedArgs.count("input") || parsedArgs.count("i"))
338  {
339  inputFolder = us::any_cast<string>(parsedArgs["input"]) + "/";
340  }
341 
342  if (parsedArgs.count("output") || parsedArgs.count("o"))
343  {
344  outputFile = us::any_cast<string>(parsedArgs["output"]);
345  }
346 
347  int skip = 0;
348  int interval = 1;
349  float opacity = .3;
350 
351  if (parsedArgs.count("skip") || parsedArgs.count("k"))
352  {
353  skip = us::any_cast<int>(parsedArgs["skip"]);
354  }
355 
356  if (parsedArgs.count("interval") || parsedArgs.count("n"))
357  {
358  interval = us::any_cast<int>(parsedArgs["interval"]);
359  }
360 
361  if (parsedArgs.count("opacity") || parsedArgs.count("c"))
362  {
363  opacity = us::any_cast<float>(parsedArgs["opacity"]);
364  }
365 
366  FileList morphFiles;
367  FileList segFiles;
368 
369  std::string refPattern;
370  std::string segPattern;
371 
372  if (parsedArgs.count("morphology") || parsedArgs.count("m"))
373  {
374  refPattern = us::any_cast<std::string>(parsedArgs["morphology"]);
375  }
376  else
377  return EXIT_FAILURE;
378 
379  if (parsedArgs.count("segmentation") || parsedArgs.count("s"))
380  {
381  segPattern = us::any_cast<std::string>(parsedArgs["segmentation"]);
382  }
383 
384  bool blank = false;
385  if (parsedArgs.count("blanked") || parsedArgs.count("b"))
386  {
387  blank = true;
388  }
390  typedef itk::Image<RGBPixelType, 2> OutputImageType;
391  typedef itk::Image<RGBPixelType, 3> InputImageType;
392 
393  mitkProgressionVisualization progressVis;
394 
395  morphFiles = CreateFileList(inputFolder, refPattern);
396  segFiles = CreateFileList(inputFolder, segPattern);
397 
398  ImageList morphImages = LoadPreprocessedFiles(morphFiles);
399  ImageList segImages;
400  if (segFiles.size() > 0 && blank == false)
401  segImages = LoadPreprocessedFiles(segFiles);
402 
403  mitk::Image::Pointer rgbImage;
404 
405  // define color for overlay image
406  mitk::Color color;
407  color.Fill(0);
408  color[0] = 200 * opacity;
409  color[1] = 0;
410 
411  // Set up itk time image filter to contain 0..N-1 images
412  itk::TileImageFilter<OutputImageType, OutputImageType>::Pointer tileFilter =
413  itk::TileImageFilter<OutputImageType, OutputImageType>::New();
414  itk::FixedArray<unsigned int, 2> layout;
415  unsigned int noOfTimeSteps = morphImages.size();
416  layout[0] = noOfTimeSteps;
417  layout[1] = 0; // automatic number of neccessary rows
418  tileFilter->SetLayout(layout);
419 
420  // Get Reference image (all images are expected to have exact same dimensions!)
421 
422  std::string fileName = morphFiles.at(0);
423  mitk::Image *referenceImg = morphImages.at(0);
424 
425  mitk::Image::Pointer merged4D;
426  std::string volumeFile;
427 
428  if (parsedArgs.count("4dVolume") || parsedArgs.count("v"))
429  {
430  volumeFile = us::any_cast<string>(parsedArgs["4dVolume"]);
431  if (volumeFile != "")
432  {
433  unsigned int *dims = new unsigned int[4];
434  dims[0] = referenceImg->GetDimensions()[0];
435  dims[1] = referenceImg->GetDimensions()[1];
436  dims[2] = referenceImg->GetDimensions()[2];
437  dims[3] = morphImages.size();
438  merged4D = mitk::Image::New();
439  merged4D->Initialize(referenceImg->GetPixelType(), 4, dims);
440  merged4D->GetTimeGeometry()->Expand(noOfTimeSteps);
441  }
442  }
443 
444  unsigned int *dim = referenceImg->GetDimensions();
445  unsigned int sliceMaxAxial = dim[2];
446 
447  // Now iterate over all data sets, extract overlay and add it to reference image
448  mitk::Image *morphImage;
449  mitk::Image *segmentationImage = nullptr;
450 
451  for (unsigned int i = 0; i < noOfTimeSteps; i++)
452  {
453  MITK_INFO << "File : " << i << " of /" << noOfTimeSteps;
454  int currentSliceNo = 0;
455 
456  // Retrieve images of current time step
457  fileName = morphFiles.at(i);
458  morphImage = morphImages.at(i);
459 
460  // Create 4D Volume image
461  if (volumeFile != "")
462  {
464 
465  merged4D->GetGeometry(i)->SetSpacing(referenceImg->GetGeometry()->GetSpacing());
466  merged4D->GetGeometry(i)->SetOrigin(referenceImg->GetGeometry()->GetOrigin());
467  merged4D->GetGeometry(i)->SetIndexToWorldTransform(referenceImg->GetGeometry()->GetIndexToWorldTransform());
468  merged4D->SetVolume(readAc.GetData(), i);
469  }
470 
471  MITK_INFO << "-- Convert to RGB";
472  rgbImage = progressVis.ConvertToRGBImage(morphImage);
473 
474  // Add current seg in red
475  color.Fill(0);
476  color[0] = 200 * opacity;
477 
478  if (!blank)
479  {
480  segmentationImage = segImages.at(i);
481  if (segmentationImage != nullptr)
482  {
483  MITK_INFO << "-- Add Overlay";
484  progressVis.AddColouredOverlay(rgbImage, segmentationImage, color);
485  }
486  }
487  // Add Segmentation of next time step in red
488  if (i == 0)
489  {
490  MITK_INFO << "Skipping retro view in first time step";
491  }
492  else
493  {
494  color.Fill(0);
495  // Add previous in green
496  color[1] = 200 * opacity;
497  if (!blank)
498  {
499  mitk::Image *nextSeg = segImages.at(i - 1);
500  MITK_INFO << "-- Add Overlay of next Step";
501  progressVis.AddColouredOverlay(rgbImage, nextSeg, color);
502  }
503  }
504 
505  // Now save all slices from overlay in output folder
506  MITK_INFO << "-- Extract Slices";
507  for (int slice = sliceMaxAxial - skip - 1; slice > skip; slice -= interval) // sliceMaxAxial/40.0f))
508  {
509  InputImageType::Pointer itkImage = InputImageType::New();
510  InputImageType::Pointer itkImage2;
511 
512  mitk::CastToItkImage(rgbImage, itkImage);
513  typedef itk::ImageDuplicator<InputImageType> DuplicatorType;
514  DuplicatorType::Pointer duplicator = DuplicatorType::New();
515  duplicator->SetInputImage(itkImage);
516  duplicator->Update();
517  itkImage2 = duplicator->GetOutput();
518 
519  itk::Image<RGBPixelType, 2>::Pointer extractedSlice = itk::Image<RGBPixelType, 2>::New();
520  extractedSlice->Graft(progressVis.ExtractSlice(itkImage2, slice));
521  tileFilter->SetInput(((currentSliceNo + 1) * (noOfTimeSteps) + i), extractedSlice);
522 
523  tileFilter->SetInput(i, progressVis.TextAsImage(GetDate(fileName)));
524 
525  currentSliceNo++;
526  }
527  }
528 
529  MITK_INFO << "Tile Filter Update";
530  tileFilter->Update();
531 
532  // Write the output image
533  typedef itk::ImageFileWriter<OutputImageType> WriterType;
534  WriterType::Pointer writer = WriterType::New();
535  writer->SetInput(tileFilter->GetOutput());
536  std::string patientName;
537 
538  patientName = GetName(fileName);
539 
540  if (blank)
541  writer->SetFileName(outputFile);
542  else
543  writer->SetFileName(outputFile);
544 
545  writer->Update();
546  if (volumeFile != "")
547  mitk::IOUtil::Save(merged4D, volumeFile);
548 
549  return EXIT_SUCCESS;
550 }
Gives locked and index-based write access for a particular image part. The class provides several set...
Gives locked and index-based read access for a particular image part. The class provides several set-...
ImageVectorType images
#define MITK_INFO
Definition: mitkLogMacros.h:18
static std::string GetName(std::string fileName, std::string suffix="_T2.nrrd")
const mitk::PixelType GetPixelType(int n=0) const
Returns the PixelType of channel n.
Definition: mitkImage.cpp:101
unsigned int * GetDimensions() const
Get the sizes of all dimensions as an integer-array.
Definition: mitkImage.cpp:1309
void setContributor(std::string contributor)
std::vector< mitk::Image::Pointer > ImageList
int main(int argc, char *argv[])
STL namespace.
ValueType * any_cast(Any *operand)
Definition: usAny.h:377
void addArgument(const std::string &longarg, const std::string &shortarg, Type type, const std::string &argLabel, const std::string &argHelp=std::string(), const us::Any &defaultValue=us::Any(), bool optional=true, bool ignoreRest=false, bool deprecated=false, mitkCommandLineParser::Channel channel=mitkCommandLineParser::Channel::None)
std::vector< std::string > FileList
std::map< std::string, us::Any > parseArguments(const StringContainerType &arguments, bool *ok=nullptr)
static std::string GetDate(std::string fileName, std::string suffix="_T2.nrrd")
class ITK_EXPORT Image
Image class for storing images.
Definition: mitkImage.h:72
Definition: usAny.h:163
const TPixel * GetData() const
Gives const access to the data.
std::string helpText() const
static ImageList LoadPreprocessedFiles(FileList files)
void setCategory(std::string category)
static FileList CreateFileList(std::string folder, std::string postfix)
Create list of all files in provided folder ending with same postfix.
const Point3D GetOrigin() const
Get the origin, e.g. the upper-left corner of the plane.
static Pointer New()
itk::RGBPixel< float > Color
Color Standard RGB color typedef (float)
void SetPixelByIndex(const itk::Index< VDimension > &idx, const TPixel &value)
Sets a pixel value at given index.
void setArgumentPrefix(const std::string &longPrefix, const std::string &shortPrefix)
static T min(T x, T y)
Definition: svm.cpp:53
itk::Image< itk::RGBPixel< PixelType >, 2 > RGB2DImage
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
itk::RGBPixel< PixelType > RGBPixelType
static void Save(const mitk::BaseData *data, const std::string &path, bool setPathProperty=false)
Save a mitk::BaseData instance.
Definition: mitkIOUtil.cpp:774
const RegionType & GetLargestPossibleRegion() const
unsigned short PixelType
void setTitle(std::string title)
double InputPixelType
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:138
void setDescription(std::string description)
mitk::AffineTransform3D * GetIndexToWorldTransform()
Get the transformation used to convert from index to world coordinates.