How to generate MPR Reconstruction for PET / CT images?

I am trying to edit the files

"Programs / dicomtodicom.cxx"
"Source / vtkDICOMPETGenerator.h"
"Source / vtkDICOMPETGenerator.cxx" to add MPR rebuild option for PT mode

below the source codes:

/*=========================================================================

vtkDICOMPETGenerator.h

=========================================================================*/
#ifndef vtkDICOMPETGenerator_h
#define vtkDICOMPETGenerator_h

#include "vtkDICOMModule.h" // For export macro
#include "vtkDICOMGenerator.h"

//! Generate DICOM data objects for PET images.
/*!
 *  Generate a DICOM data set belonging to one of the PET SOP classes.
 *  The assumption is that images have been read into VTK as PET, and
 *  are being written out as derived images after being processed.
 *  The specific IOD classes supported are as follows:
 *  - PET Image, 1.2.840.10008.5.1.4.1.1.128
 */
class VTKDICOM_EXPORT vtkDICOMPETGenerator : public vtkDICOMGenerator
{
public:
  //! Static method for construction.
  static vtkDICOMPETGenerator *New();
  vtkTypeMacro(vtkDICOMPETGenerator, vtkDICOMGenerator);

  //! Print information about this object.
  void PrintSelf(ostream& os, vtkIndent indent) VTK_DICOM_OVERRIDE;

  //! Generate an instance of one of the supported classes.
  /*!
   *  This is the primary interface method of this class.  Given the
   *  information for a vtkImageData object, it will populate the
   *  attributes of the supplied vtkDICOMMetaData object.
   */
  bool GenerateInstance(vtkInformation *info) VTK_DICOM_OVERRIDE;

protected:
  vtkDICOMPETGenerator();
  ~vtkDICOMPETGenerator();

  //! Generate the Series Module.
  virtual bool GeneratePETSeriesModule(vtkDICOMMetaData *source);

  //! Generate the Image Module.
  virtual bool GeneratePETImageModule(vtkDICOMMetaData *source);

  //! Instantiate a DICOM Secondary Capture Image object.
  virtual bool GeneratePETInstance(vtkInformation *info);

private:
#ifdef VTK_DICOM_DELETE
  vtkDICOMPETGenerator(const vtkDICOMPETGenerator&) VTK_DICOM_DELETE;
  void operator=(const vtkDICOMPETGenerator&) VTK_DICOM_DELETE;
#else
  vtkDICOMPETGenerator(const vtkDICOMPETGenerator&) = delete;
  void operator=(const vtkDICOMPETGenerator&) = delete;
#endif
};

#endif // vtkDICOMPETGenerator_h

/*=========================================================================
vtkDICOMPETGenerator.cxx
=========================================================================*/
#include "vtkDICOMPETGenerator.h"
#include "vtkDICOMMetaData.h"
#include "vtkDICOMSequence.h"
#include "vtkDICOMItem.h"
#include "vtkDICOMTagPath.h"

#include "vtkObjectFactory.h"
#include "vtkDataSetAttributes.h"
#include "vtkInformation.h"

vtkStandardNewMacro(vtkDICOMPETGenerator);

//----------------------------------------------------------------------------
vtkDICOMPETGenerator::vtkDICOMPETGenerator()
{
}

//----------------------------------------------------------------------------
vtkDICOMPETGenerator::~vtkDICOMPETGenerator()
{
}

//----------------------------------------------------------------------------
void vtkDICOMPETGenerator::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
}

//----------------------------------------------------------------------------
bool vtkDICOMPETGenerator::GeneratePETSeriesModule(vtkDICOMMetaData *source)
{
  vtkDICOMMetaData *meta = this->MetaData;
  meta->Set(DC::Modality, "PT");

  // optional and conditional: direct copy of values with no checks
  static const DC::EnumType optional[] = {
    DC::ReferencedPerformedProcedureStepSequence, // 1C
    DC::ItemDelimitationItem
  };

  return this->CopyOptionalAttributes(optional, source);
}

//----------------------------------------------------------------------------
bool vtkDICOMPETGenerator::GeneratePETImageModule(vtkDICOMMetaData *source)
{
  // ImageType is specialized from GeneralImageModule,
  // by adding a third value that is specific to CT:
  // AXIAL or LOCALIZER
  const char *it = 0;
  if (source)
  {
    it = source->Get(DC::ImageType).GetCharData();
  }
  if (it == 0 || it[0] == '\0')
  {
    it = "DERIVED\\SECONDARY\\AXIAL";
  }
  vtkDICOMMetaData *meta = this->MetaData;
  meta->Set(DC::ImageType, it);

  // These specialized from ImagePixelModule:
  // SamplesPerPixel must be 1
  // PhotometricInterpretation must be MONOCHROME1 or MONOCHROME2
  // BitsAllocated must be 16
  // BitsStored must be 12 to 16

  // these are mandatory
  meta->Set(DC::RescaleIntercept, this->RescaleIntercept);
  meta->Set(DC::RescaleSlope, this->RescaleSlope);

  // SpacingBetweenSlices is not part of the CT IOD, but it is a
  // frequently used extended attribute.  If it is present we want to
  // ensure that it is correct.  If absent, we don't set it.
  if (meta->Has(DC::SpacingBetweenSlices))
  {
    meta->Set(DC::SpacingBetweenSlices, this->Spacing[2]);
  }

  // required items: use simple read/write validation
  DC::EnumType required[] = {
    DC::KVP,
    DC::AcquisitionNumber,
    DC::ItemDelimitationItem
  };

  // optional and conditional: direct copy of values with no checks
  static const DC::EnumType optional[] = {
    DC::RescaleType, // 1C
    DC::ScanOptions,
    DC::DataCollectionDiameter,
    DC::DataCollectionCenterPatient,
    DC::ReconstructionDiameter,
    DC::ReconstructionTargetCenterPatient,
    DC::DistanceSourceToDetector,
    DC::DistanceSourceToPatient,
    DC::GantryDetectorTilt,
    DC::TableHeight,
    DC::RotationDirection,
    DC::ExposureTime,
    DC::XRayTubeCurrent,
    DC::Exposure,
    DC::ExposureInuAs,
    DC::FilterType,
    DC::GeneratorPower,
    DC::FocalSpots,
    DC::ConvolutionKernel,
    DC::RevolutionTime,
    DC::SingleCollimationWidth,
    DC::TotalCollimationWidth,
    DC::TableSpeed,
    DC::TableFeedPerRotation,
    DC::ExposureModulationType,
    DC::AnatomicRegionSequence,
    DC::PrimaryAnatomicStructureSequence,
    DC::ItemDelimitationItem
  };

  return (this->CopyRequiredAttributes(required, source) &&
          this->CopyOptionalAttributes(optional, source));
}

//----------------------------------------------------------------------------
bool vtkDICOMPETGenerator::GeneratePETInstance(vtkInformation *info)
{
  this->SetPixelRestrictions(
    RepresentationSigned | RepresentationUnsigned,
    BitsStored12 | BitsStored16,
    1);

  const char *SOPClass = "1.2.840.10008.5.1.4.1.1.128";
  this->InitializeMetaData(info);

  vtkDICOMMetaData *source = this->SourceMetaData;
  if (!this->GenerateSOPCommonModule(source, SOPClass) ||
      !this->GeneratePatientModule(source) ||
      !this->GenerateClinicalTrialSubjectModule(source) ||
      !this->GenerateGeneralStudyModule(source) ||
      !this->GeneratePatientStudyModule(source) ||
      !this->GenerateClinicalTrialStudyModule(source) ||
      !this->GenerateGeneralSeriesModule(source) ||
      !this->GenerateFrameOfReferenceModule(source) ||
      !this->GenerateClinicalTrialSeriesModule(source) ||
      !this->GenerateGeneralEquipmentModule(source) ||
      !this->GenerateGeneralImageModule(source) ||
      !this->GenerateImagePlaneModule(source) ||
      !this->GenerateImagePixelModule(source) ||
      !this->GenerateContrastBolusModule(source) ||
      !this->GenerateDeviceModule(source) ||
      !this->GenerateSpecimenModule(source) ||
      !this->GeneratePETSeriesModule(source) ||
      !this->GeneratePETImageModule(source) ||
      !this->GenerateOverlayPlaneModule(source) ||
      !this->GenerateVOILUTModule(source))
  {
    return false;
  }

  return true;
}

//----------------------------------------------------------------------------
bool vtkDICOMPETGenerator::GenerateInstance(vtkInformation *info)
{
  if (this->MultiFrame)
  {
    vtkErrorMacro("Enhanced Multi-Frame PET is not yet supported.");
    return false;
  }

  return this->GeneratePETInstance(info);
}

/*=========================================================================
 dicomtodicom.cxx
=========================================================================*/

#include "vtkDICOMConfig.h"
#include "vtkDICOMBuild.h"
#include "vtkDICOMMetaData.h"
#include "vtkDICOMParser.h"
#include "vtkDICOMReader.h"
#include "vtkDICOMWriter.h"
#include "vtkDICOMFileSorter.h"
#include "vtkDICOMMRGenerator.h"
#include "vtkDICOMPETGenerator.h"
#include "vtkDICOMCTGenerator.h"
#include "vtkDICOMToRAS.h"
#include "vtkDICOMPETRectifier.h"
#include "vtkDICOMUtilities.h"
#include "vtkDICOMFile.h"
#include "vtkDICOMFileDirectory.h"

#include "vtkImageData.h"
#include "vtkMatrix4x4.h"
#include "vtkImageReslice.h"
#include "vtkStringArray.h"
#include "vtkIntArray.h"
#include "vtkErrorCode.h"
#include "vtkSortFileNames.h"
#include "vtkSmartPointer.h"
#include "vtkVersion.h"

#if VTK_MAJOR_VERSION >= 6 || VTK_MINOR_VERSION >= 10
#include "vtkImageResize.h"
#include "vtkImageSincInterpolator.h"
#endif

#include <algorithm>
#include <string>
#include <vector>
#include <set>

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>

// from dicomcli
#include "vtkConsoleOutputWindow.h"
#include "mainmacro.h"

// Kinds of reformats
enum MPREnum
{
  MPRAxial = 1,
  MPRCoronal = 2,
  MPRSagittal = 3
};

static const char *dicomtodicom_description[] = {
  " NEW", " MPR Ax", " MPR Cor", " MPR Sag"
};

// Simple structure for command-line options
struct dicomtodicom_options
{
  const char *modality;
  const char *series_description;
  const char *series_number;
  const char *uid_prefix;
  int mpr;
  bool resample;
  bool silent;
  bool verbose;
  const char *output;
};

// Print the version
void dicomtodicom_version(FILE *file, const char *command_name, bool verbose)
{
  const char *cp = command_name + strlen(command_name);
  while (cp != command_name && cp[-1] != '\\' && cp[-1] != '/') { --cp; }

  if (!verbose)
  {
    fprintf(file, "%s %s\n", cp, DICOM_VERSION);
    fprintf(file, "\n"
      "Copyright (c) 2012-2019, David Gobbi.\n\n"
      "This software is distributed under an open-source license.  See the\n"
      "Copyright.txt file that comes with the vtk-dicom source distribution.\n");
  }
  else
  {
    fprintf(file,
      "Head %8.8s, Built %s, %s\n",
      DICOM_SOURCE_VERSION, DICOM_BUILD_DATE, DICOM_BUILD_TIME);
  }
}

// Print the options
void dicomtodicom_usage(FILE *file, const char *command_name)
{
  const char *cp = command_name + strlen(command_name);
  while (cp != command_name && cp[-1] != '\\' && cp[-1] != '/') { --cp; }

  fprintf(file,
    "usage:\n"
    "  %s -o directory file.dcm ...\n\n", cp);
  fprintf(file,
    "options:\n"
    "  -o directory            The output directory.\n"
    "  -s --silent             Do not print anything while executing.\n"
    "  -v --verbose            Verbose error reporting.\n"
    "  --resample              Resample to produce square pixels.\n"
    "  --axial                 Produce axial slices.\n"
    "  --coronal               Produce coronal slices.\n"
    "  --sagittal              Produce sagittal slices.\n"
    "  --series-description    Textual description of the series.\n"
    "  --series-number         The series number to use.\n"
    "  --modality              The modality: MR or CT or SC or PT.\n"
    "  --uid-prefix            A DICOM uid prefix (optional).\n"
    "  --version               Print the version and exit.\n"
    "  --build-version         Print source and build version.\n"
    "  --help                  Documentation for dicomtodicom.\n"
  );
}

// Print the help
void dicomtodicom_help(FILE *file, const char *command_name)
{
  dicomtodicom_usage(file, command_name);

  fprintf(file,
    "\n");

  fprintf(file,
    "This program allows adjustment of a DICOM series.\n"
    "\n");
  fprintf(file,
    "In its most basic functionality, this command reads a DICOM series and\n"
    "then writes it out with a new UID.  It strips out any meta data that is\n"
    "not recognized as part of the SOP class that it is writing.  In other\n"
    "words, its main purpose is to coerce the data to fit the requested SOP\n"
    "class.\n"
    "\n");
  fprintf(file,
    "Currently, only three output SOP classes can be written: Secondary\n"
    "Capture, CT, PT and MR.  Enhanced CT, PT and Enhanced MR cannot be written,\n"
    "but they can be read, therefore this program can be used to convert\n"
    "enhanced single-file CT, PT and MR DICOM into a series of DICOM files.\n"
    "\n");
  fprintf(file,
    "The written data has the ImageType set to DERIVED\\SECONDARY\\OTHER\n"
    "and has a new series number and name.  It isn't meant to replace the\n"
    "original data, it is simply meant to coerce the data into a format\n"
    "that might be more compatible with other software.\n"
    "\n");
  fprintf(file,
    "Reformatting of the data (MPR) is permitted during the conversion.\n"
    "This is an experimental feature and causes much of the per-instance\n"
    "meta data to be discarded.  Reformatting be combined with resampling\n"
    "to produce an output with square pixels via Lanczos interpolation.\n"
    "\n");
}

// Print error
void dicomtodicom_check_error(vtkObject *o)
{
  vtkDICOMReader *reader = vtkDICOMReader::SafeDownCast(o);
  vtkDICOMFileSorter *sorter = vtkDICOMFileSorter::SafeDownCast(o);
  vtkDICOMWriter *writer = vtkDICOMWriter::SafeDownCast(o);
  vtkDICOMParser *parser = vtkDICOMParser::SafeDownCast(o);
  const char *filename = 0;
  unsigned long errorcode = 0;
  if (writer)
  {
    filename = writer->GetFileName();
    errorcode = writer->GetErrorCode();
  }
  else if (reader)
  {
    filename = reader->GetInternalFileName();
    errorcode = reader->GetErrorCode();
  }
  else if (sorter)
  {
    filename = sorter->GetInternalFileName();
    errorcode = sorter->GetErrorCode();
  }
  else if (parser)
  {
    filename = parser->GetFileName();
    errorcode = parser->GetErrorCode();
  }
  if (!filename)
  {
    filename = "";
  }

  switch(errorcode)
  {
    case vtkErrorCode::NoError:
      return;
    case vtkErrorCode::FileNotFoundError:
      fprintf(stderr, "File not found: %s\n", filename);
      break;
    case vtkErrorCode::CannotOpenFileError:
      fprintf(stderr, "Cannot open file: %s\n", filename);
      break;
    case vtkErrorCode::UnrecognizedFileTypeError:
      fprintf(stderr, "Unrecognized file type: %s\n", filename);
      break;
    case vtkErrorCode::PrematureEndOfFileError:
      fprintf(stderr, "File is truncated: %s\n", filename);
      break;
    case vtkErrorCode::FileFormatError:
      fprintf(stderr, "Bad DICOM file: %s\n", filename);
      break;
    case vtkErrorCode::NoFileNameError:
      fprintf(stderr, "Output filename could not be used: %s\n", filename);
      break;
    case vtkErrorCode::OutOfDiskSpaceError:
      fprintf(stderr, "Out of disk space while writing file: %s\n", filename);
      break;
    default:
      fprintf(stderr, "An unknown error occurred.\n");
      break;
  }

  exit(1);
}

// Read the options
void dicomtodicom_read_options(
  int argc, char *argv[],
  dicomtodicom_options *options, vtkStringArray *files)
{
  options->mpr = 0;
  options->modality = 0;
  options->series_description = 0;
  options->series_number = 0;
  options->uid_prefix = "2.25";
  options->resample = false;
  options->silent = false;
  options->verbose = false;
  options->output = 0;

  // read the options from the command line
  int argi = 1;
  while (argi < argc)
  {
    const char *arg = argv[argi++];
    if (arg[0] == '-')
    {
      if (strcmp(arg, "--") == 0)
      {
        // stop processing switches
        break;
      }
      else if (strcmp(arg, "--modality") == 0 ||
               strcmp(arg, "--series-description") == 0 ||
               strcmp(arg, "--series-number") == 0 ||
               strcmp(arg, "--uid-prefix") == 0)
      {
        if (argi >= argc ||
            argv[argi][0] == '-')
        {
          fprintf(stderr, "\nA value must follow the \'%s\' flag\n\n", arg);
          exit(1);
        }
        if (strcmp(arg, "--modality") == 0)
        {
          options->modality = argv[argi];
        }
        else if (strcmp(arg, "--series-description") == 0)
        {
          options->series_description = argv[argi];
        }
        else if (strcmp(arg, "--series-number") == 0)
        {
          options->series_number= argv[argi];
        }
        else if (strcmp(arg, "--uid-prefix") == 0)
        {
          options->uid_prefix = argv[argi];
        }
        argi++;
      }
      else if (strcmp(arg, "--axial") == 0)
      {
        options->mpr = MPRAxial;
      }
      else if (strcmp(arg, "--coronal") == 0)
      {
        options->mpr = MPRCoronal;
      }
      else if (strcmp(arg, "--sagittal") == 0)
      {
        options->mpr = MPRSagittal;
      }
      else if (strcmp(arg, "--resample") == 0)
      {
        options->resample = true;
      }
      else if (strcmp(arg, "--silent") == 0)
      {
        options->silent = true;
      }
      else if (strcmp(arg, "--verbose") == 0)
      {
        options->verbose = true;
      }
      else if (strcmp(arg, "--version") == 0)
      {
        dicomtodicom_version(stdout, argv[0], false);
        exit(0);
      }
      else if (strcmp(arg, "--build-version") == 0)
      {
        dicomtodicom_version(stdout, argv[0], true);
        exit(0);
      }
      else if (strcmp(arg, "--help") == 0)
      {
        dicomtodicom_help(stdout, argv[0]);
        exit(0);
      }
      else if (arg[0] == '-' && arg[1] == '-')
      {
        fprintf(stderr, "\nUnrecognized option %s\n\n", arg);
        dicomtodicom_usage(stderr, argv[0]);
        exit(1);
      }
      else if (arg[0] == '-' && arg[1] != '-')
      {
        for (int argj = 1; arg[argj] != '\0'; argj++)
        {
          if (arg[argj] == 's')
          {
            options->silent = true;
          }
          else if (arg[argj] == 'v')
          {
            options->verbose = true;
          }
          else if (arg[argj] == 'o')
          {
            if (arg[argj+1] != '\0')
            {
              arg += argj+1;
            }
            else
            {
              if (argi >= argc)
              {
                fprintf(stderr, "\nA file must follow the \'-o\' flag\n\n");
                dicomtodicom_usage(stderr, argv[0]);
                exit(1);
              }
              arg = argv[argi++];
            }
            options->output = arg;
            break;
          }
          else
          {
            fprintf(stderr, "\nUnrecognized \'%c\' in option %s\n\n", arg[argj], arg);
            dicomtodicom_usage(stderr, argv[0]);
            exit(1);
          }
        }
      }
    }
    else
    {
      files->InsertNextValue(arg);
    }
  }

  while (argi < argc)
  {
    files->InsertNextValue(argv[argi++]);
  }
}

// Convert one DICOM series into another DICOM series
void dicomtodicom_convert_one(
  const dicomtodicom_options *options,
  vtkStringArray *a,
  const char *outfile)
{
  // read the files
  vtkSmartPointer<vtkDICOMReader> reader =
    vtkSmartPointer<vtkDICOMReader>::New();
  reader->SetMemoryRowOrderToFileNative();
  reader->TimeAsVectorOn();
  reader->SetFileNames(a);
  reader->Update();
  dicomtodicom_check_error(reader);

  // get a handle for the reader's output
  vtkAlgorithmOutput *lastOutput = reader->GetOutputPort();

  // The meta data object
  vtkSmartPointer<vtkDICOMMetaData> meta =
    vtkSmartPointer<vtkDICOMMetaData>::New();
  meta->DeepCopy(reader->GetMetaData());
  meta->Set(DC::SeriesNumber, meta->Get(DC::SeriesNumber).AsUnsignedInt() +
    1000*(1 + options->mpr*100));
  std::string seriesDescription =
    (meta->Get(DC::SeriesDescription).AsString() +
     dicomtodicom_description[options->mpr]);
  if (seriesDescription.size() < 64)
  {
    meta->Set(DC::SeriesDescription, seriesDescription);
  }

  // set the metadata supplied on the command line
  if (options->series_description)
  {
    meta->Set(DC::SeriesDescription, options->series_description);
  }
  if (options->series_number)
  {
    meta->Set(DC::SeriesNumber, options->series_number);
  }

  // get the matrix from the DICOM series
  vtkMatrix4x4 *inputMatrix = reader->GetPatientMatrix();
  vtkSmartPointer<vtkMatrix4x4> matrix =
    vtkSmartPointer<vtkMatrix4x4>::New();
  if (inputMatrix)
  {
    matrix->DeepCopy(inputMatrix);
  }

  // mpr reformat if requested
  vtkSmartPointer<vtkImageReslice> reformat =
    vtkSmartPointer<vtkImageReslice>::New();
  vtkSmartPointer<vtkDICOMPETRectifier> rectify =
    vtkSmartPointer<vtkDICOMPETRectifier>::New();
#if VTK_MAJOR_VERSION >= 6 || VTK_MINOR_VERSION >= 10
  vtkSmartPointer<vtkImageResize> resample =
    vtkSmartPointer<vtkImageResize>::New();
#endif
  vtkSmartPointer<vtkMatrix4x4> axes =
    vtkSmartPointer<vtkMatrix4x4>::New();
  int permutation[3] = { 0, 1, 2 };

  if (options->mpr)
  {
    // check for CT acquired with a tilted gantry
    if (fabs(vtkDICOMPETRectifier::GetGantryDetectorTilt(matrix)) > 1e-2)
    {
      // tilt is significant, so regrid as a rectangular volume
      rectify->SetInputConnection(lastOutput);
      rectify->SetVolumeMatrix(matrix);
      rectify->Update();
      lastOutput = rectify->GetOutputPort();
      matrix = rectify->GetRectifiedMatrix();
    }

    // check if resampling was requested
    if (options->resample)
    {
#if VTK_MAJOR_VERSION >= 6 || VTK_MINOR_VERSION >= 10
      // generate cube voxels
      double spacing[3] = { 1.0, 1.0, 1.0 };
      const vtkDICOMValue& v = meta->Get(DC::PixelSpacing);
      if (v.GetNumberOfValues() == 2)
      {
        v.GetValues(spacing, 2);
        double s = std::min(spacing[0], spacing[1]);
        spacing[0] = s;
        spacing[1] = s;
        spacing[2] = s;
      }
      resample->SetInputConnection(lastOutput);
      vtkSmartPointer<vtkImageSincInterpolator> interpolator =
        vtkSmartPointer<vtkImageSincInterpolator>::New();
      interpolator->SetWindowFunctionToLanczos();
      resample->SetInterpolator(interpolator);
      resample->SetResizeMethodToOutputSpacing();
      resample->SetOutputSpacing(spacing);
      resample->InterpolateOn();
      resample->BorderOn();
      resample->Update();
      lastOutput = resample->GetOutputPort();
#else
      fprintf(stderr,
              "\nTo use --resample, recompile with VTK 5.10 or later.\n\n");
#endif
    }

    // create a permutation matrix to make the slices axial
    axes->DeepCopy(matrix);
    axes->Invert();
    int maxidx[3] = { -1, -1, -1 };
    double value[3] = { 1.0, 1.0, 1.0 };
    int prevmaxj = -1;
    int prevmaxi = -1;
    for (int kdim = 0; kdim < 2; kdim++)
    {
      int maxj = 0;
      int maxi = 0;
      double maxv = -0.0;
      for (int jdim = 0; jdim < 3; jdim++)
      {
        if (jdim == prevmaxj) { continue; }
        for (int idim = 0; idim < 3; idim++)
        {
          if (idim == prevmaxi) { continue; }
          double v = axes->GetElement(jdim, idim);
          if (v*v >= maxv)
          {
            maxi = idim;
            maxj = jdim;
            maxv = v*v;
          }
        }
      }
      maxidx[maxj] = maxi;
      value[maxj] = (axes->GetElement(maxj, maxi) < 0 ? -1.0 : 1.0);
      prevmaxj = maxj;
      prevmaxi = maxi;
    }

    axes->Zero();
    axes->SetElement(3, 3, 1.0);
    for (int jdim = 0; jdim < 3; jdim++)
    {
      int idim = maxidx[jdim];
      if (idim < 0)
      {
        idim = 3 - maxidx[(jdim+1)%3] - maxidx[(jdim+2)%3];
        maxidx[jdim] = idim;
        double perm = (((3 + maxidx[2] - maxidx[0])%3) == 2 ? 1.0 : -1.0);
        value[jdim] = value[(jdim+1)%3]*value[(jdim+2)%3]*perm;
      }
      permutation[jdim] = idim;
      axes->SetElement(jdim, idim, value[jdim]);
    }

    // change the permutation to the desired mpr
    if (options->mpr == MPRCoronal)
    {
      double cmatrix[16] = {
        1.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 1.0, 0.0,
        0.0,-1.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 1.0 };
      vtkMatrix4x4::Multiply4x4(*axes->Element, cmatrix, *axes->Element);
      int tperm[3] = { permutation[0], permutation[1], permutation[2] };
      permutation[0] = tperm[0];
      permutation[1] = tperm[2];
      permutation[2] = tperm[1];
    }
    else if (options->mpr == MPRSagittal)
    {
      double smatrix[16] = {
        0.0, 0.0,-1.0, 0.0,
        1.0, 0.0, 0.0, 0.0,
        0.0,-1.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 1.0 };
      vtkMatrix4x4::Multiply4x4(*axes->Element, smatrix, *axes->Element);
      int tperm[3] = { permutation[0], permutation[1], permutation[2] };
      permutation[0] = tperm[1];
      permutation[1] = tperm[2];
      permutation[2] = tperm[0];
    }

    // reformat with the permuted axes
    reformat->SetResliceAxes(axes);
    reformat->SetInputConnection(lastOutput);
    lastOutput = reformat->GetOutputPort();

    // factor out the permuted axes
    vtkMatrix4x4::Multiply4x4(matrix, axes, matrix);
  }

  // make the generator
  vtkSmartPointer<vtkDICOMMRGenerator> mrgenerator =
    vtkSmartPointer<vtkDICOMMRGenerator>::New();
  vtkSmartPointer<vtkDICOMCTGenerator> ctgenerator =
    vtkSmartPointer<vtkDICOMCTGenerator>::New();
  vtkDICOMGenerator *generator = 0;

  vtkSmartPointer<vtkDICOMPETGenerator> petgenerator =
    vtkSmartPointer<vtkDICOMPETGenerator>::New();
  vtkDICOMGenerator *generator = 0;

  // get the generator from the supplied DICOM data
  std::string SOPClass = meta->Get(DC::SOPClassUID).AsString();
  if (SOPClass == "1.2.840.10008.5.1.4.1.1.2" ||
      SOPClass == "1.2.840.10008.5.1.4.1.1.2.1" ||
      SOPClass == "1.2.840.10008.5.1.4.1.1.2.2")
  {
    generator = ctgenerator;
  }
  else if (SOPClass == "1.2.840.10008.5.1.4.1.1.4" ||
           SOPClass == "1.2.840.10008.5.1.4.1.1.4.1" ||
           SOPClass == "1.2.840.10008.5.1.4.1.1.4.4")
  {
    generator = mrgenerator;
  }
  else if (SOPClass == "1.2.840.10008.5.1.4.1.1.128" ||
           SOPClass == "1.2.840.10008.5.1.4.1.1.130" ||
           SOPClass == "1.2.840.10008.5.1.4.1.1.128.1")
  {
    generator = petgenerator;
  }

  // allow the user to override the generator
  if (options->modality)
  {
    if (strcmp(options->modality, "CT") == 0)
    {
      generator = ctgenerator;
    }
    else if (strcmp(options->modality, "MR") == 0 ||
             strcmp(options->modality, "MRI") == 0)
    {
      generator = mrgenerator;
    }
    else if (strcmp(options->modality, "PT") == 0 ||
             strcmp(options->modality, "NM") == 0)
    {
      generator = petgenerator;
    }

  }

  // prepare the writer to write the image
  vtkSmartPointer<vtkDICOMWriter> writer =
    vtkSmartPointer<vtkDICOMWriter>::New();
  if (generator)
  {
    writer->SetGenerator(generator);
  }
  writer->SetMetaData(meta);
  writer->SetFilePrefix(outfile);
  writer->SetFilePattern("%s/IM-0001-%04.4d.dcm");
  writer->TimeAsVectorOn();
  if (reader->GetTimeDimension() > 1)
  {
    writer->SetTimeDimension(reader->GetTimeDimension());
    writer->SetTimeSpacing(reader->GetTimeSpacing());
  }
  writer->SetPatientMatrix(matrix);
  if (reader->GetRescaleSlope() > 0)
  {
    writer->SetRescaleSlope(reader->GetRescaleSlope());
    writer->SetRescaleIntercept(reader->GetRescaleIntercept());
  }
  writer->SetInputConnection(lastOutput);
  writer->SetMemoryRowOrderToFileNative();
  writer->Write();
  dicomtodicom_check_error(writer);
}

// Process a list of files
void dicomtodicom_convert_files(
  dicomtodicom_options *options, vtkStringArray *files,
  const char *outpath)
{
  // sort the files by filename first, as a fallback
  vtkSmartPointer<vtkSortFileNames> presorter =
    vtkSmartPointer<vtkSortFileNames>::New();
  presorter->NumericSortOn();
  presorter->IgnoreCaseOn();
  presorter->SetInputFileNames(files);
  presorter->Update();

  // sort the files by study and series
  vtkSmartPointer<vtkDICOMFileSorter> sorter =
    vtkSmartPointer<vtkDICOMFileSorter>::New();
  sorter->SetInputFileNames(presorter->GetFileNames());
  sorter->Update();
  dicomtodicom_check_error(sorter);

  dicomtodicom_convert_one(
    options, sorter->GetOutputFileNames(), outpath);
}

// This program will convert DICOM to DICOM
int MAINMACRO(int argc, char *argv[])
{
  // redirect all VTK errors to stderr
  vtkConsoleOutputWindow::Install();

  // for the list of input DICOM files
  vtkSmartPointer<vtkStringArray> files =
    vtkSmartPointer<vtkStringArray>::New();

  dicomtodicom_options options;
  dicomtodicom_read_options(argc, argv, &options, files);

  // whether to silence VTK warnings and errors
  vtkObject::SetGlobalWarningDisplay(options.verbose);

  // set the UID prefix
  if (options.uid_prefix)
  {
    vtkDICOMUtilities::SetUIDPrefix(options.uid_prefix);
  }

  // make sure that input files were provided
  if (files->GetNumberOfValues() == 0)
  {
    fprintf(stderr, "\nNo input files were specified.\n\n");
    dicomtodicom_usage(stderr, argv[0]);
    exit(1);
  }

  // the output
  const char *outpath = options.output;
  if (!outpath)
  {
    fprintf(stderr,
      "\nNo output directory was specified (\'-o\' <directory>).\n\n");
    dicomtodicom_usage(stderr, argv[0]);
    exit(1);
  }

  int code = vtkDICOMFile::Access(outpath, vtkDICOMFile::In);
  if (code != vtkDICOMFile::FileIsDirectory)
  {
    fprintf(stderr, "option -o must give a directory, not a file.\n");
    exit(1);
  }
  code = vtkDICOMFileDirectory::Create(outpath);
  if (code != vtkDICOMFileDirectory::Good)
  {
    fprintf(stderr, "Cannot create directory: %s\n", outpath);
    exit(1);
  }

  dicomtodicom_convert_files(&options, files, outpath);

  return 0;
}

Is my reasoning correct?

dgobbi / vtk-dicom

How to generate MPR Reconstruction for PET / CT images? #180