/*=========================================================================
  *
  *  Copyright RTK Consortium
  *
  *  Licensed under the Apache License, Version 2.0 (the "License");
  *  you may not use this file except in compliance with the License.
  *  You may obtain a copy of the License at
  *
  *         http://www.apache.org/licenses/LICENSE-2.0.txt
  *
  *  Unless required by applicable law or agreed to in writing, software
  *  distributed under the License is distributed on an "AS IS" BASIS,
  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  *  See the License for the specific language governing permissions and
  *  limitations under the License.
  *
  *=========================================================================*/

  //RTK inclusions
  #include "rtkJosephForwardProjectionImageFilter.h"
  #include <rtkThreeDCircularProjectionGeometry.h>
  #include <rtkConstantImageSource.h>

  //ITK inclusions
  #include <itkImageFileReader.h>
  #include <itkImageFileWriter.h>
  #include <itkExtractImageFilter.h>

  //Inclusions
  #include <fstream>
  #include <iostream>
  #include <vector>

  int
  main()
  {
  // parameter list
  double sid_arg = 1000;
  double sdd_arg = 1536;
  double first_angle = 0;
  double proj_iso_x_arg = 0;
  double proj_iso_y_arg = 0;
  double out_angle = 0;
  double in_angle = 0;
  double source_x_arg = 0;
  double source_y_arg = 0;
  double rad_cyl = 0;
  const std::string input_arg = "target.mhd";
  const std::string output_arg = "projectedImage.png";
  std::vector<double> rotations;
  double input;

  using InputPixelType = unsigned short;
  using OutputPixelType = unsigned short;

  using InputImageType = itk::Image<InputPixelType, 3 >;
  using OutputImageType = itk::Image<OutputPixelType, 2 >;


  // Create RTK geometry object
  using GeometryType = rtk::ThreeDCircularProjectionGeometry;
  GeometryType::Pointer geometry = GeometryType::New();
  // Projection matrix
  geometry->AddProjection(sid_arg,
                          sdd_arg,
                          first_angle,
                          proj_iso_x_arg,
                          proj_iso_y_arg,
                          out_angle,
                          in_angle,
                          source_x_arg,
                          source_y_arg);

 // Input reader
  using ReaderType = itk::ImageFileReader<InputImageType>;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(input_arg);
  TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->Update())


  // Create a stack of empty projection images
  using ConstantImageSourceType = rtk::ConstantImageSource<InputImageType>;
  ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
  ConstantImageSourceType::PointType origin;
  ConstantImageSourceType::SpacingType spacing;
  ConstantImageSourceType::SizeType sizeOutput;

  const InputImageType::PointType & og = reader->GetOutput()->GetOrigin();
   for(int i = 0; i < 3 ; i++)
  {
    origin[i] = og[i];
  }
  constantImageSource->SetOrigin(origin);

  const InputImageType::SpacingType & sp = reader->GetOutput()->GetSpacing();
  for(int i = 0; i < 3 ; i++)
  {
    spacing[i] = sp[i];
  }
  constantImageSource->SetSpacing(spacing);

  sizeOutput[0] = 512;
  sizeOutput[1] = 512;
  sizeOutput[2] = 1;
  constantImageSource->SetSize(sizeOutput);


  // Create forward projection image filter
  rtk::ForwardProjectionImageFilter<InputImageType, InputImageType>::Pointer forwardProjection;
  forwardProjection = rtk::JosephForwardProjectionImageFilter<InputImageType, InputImageType>::New();
  forwardProjection->SetInput(constantImageSource->GetOutput());
  forwardProjection->SetInput(1, reader->GetOutput());
  forwardProjection->SetGeometry(geometry);
  TRY_AND_EXIT_ON_ITK_EXCEPTION(forwardProjection->Update())
  
  // Create slice extraction filter
  using FilterType = itk::ExtractImageFilter< InputImageType,
                                      OutputImageType >;
  FilterType::Pointer filter = FilterType::New();
  filter->InPlaceOn();
  filter->SetDirectionCollapseToSubmatrix();

  forwardProjection -> UpdateOutputInformation();
  InputImageType::RegionType inputRegion =
            forwardProjection->GetOutput()->GetLargestPossibleRegion();
  InputImageType::SizeType size = inputRegion.GetSize();
  size[2] = 0;

  InputImageType::IndexType start = inputRegion.GetIndex();
  start[2] = 0;

  InputImageType::RegionType desiredRegion;
  desiredRegion.SetSize(  size  );
  desiredRegion.SetIndex( start );

  filter->SetExtractionRegion( desiredRegion );
  filter->SetInput(forwardProjection->GetOutput());

  // Write
  using WriterType = itk::ImageFileWriter<OutputImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(output_arg);
  writer->SetInput(filter->GetOutput());
  TRY_AND_EXIT_ON_ITK_EXCEPTION(writer->Update())
  
  return EXIT_SUCCESS;
  }