#include <vtkVersion.h>
#include <vtkRenderer.h>
#include <vtkImageActor.h>
#include <vtkStructuredPoints.h>
#include <vtkSmartPointer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkPNGReader.h>
#include <vtkImageLuminance.h>
#include <vtkProperty.h>
#include <vtkImageGaussianSmooth.h>
#include <vtkImageCast.h>
#include <vtkImageGradient.h>
#include <vtkImageMagnitude.h>
#include <vtkImageNonMaximumSuppression.h>
#include <vtkImageConstantPad.h>
#include <vtkImageToStructuredPoints.h>
#include <vtkLinkEdgels.h>
#include <vtkThreshold.h>
#include <vtkGeometryFilter.h>
#include <vtkSubPixelPositionEdgels.h>
#include <vtkStripper.h>
#include <vtkPolyDataMapper.h>
#include <vtkCamera.h>
#include <vtkDICOMImageReader.h>

int main(int argc, char* argv[])
{
    //if(argc != 2)
    //{
    //  std::cerr << "Required args: filename.png" << std::endl;
    //  return EXIT_FAILURE;
    //}

    //std::string filename = argv[1];
    std::string filename = "D:\\2017_Projects\\CannyEdgeDetector\\Data\\CT_IQI_2.DCM";


    // Define viewport ranges
    // (xmin, ymin, xmax, ymax)
    double originalViewport[4] = {0.0, 0.0, 0.5, 1.0};
    double edgeViewport[4] = {0.5, 0.0, 1.0, 1.0};

    vtkSmartPointer<vtkRenderer> originalRenderer =
        vtkSmartPointer<vtkRenderer>::New();
    originalRenderer->SetViewport(originalViewport);
    vtkSmartPointer<vtkRenderer> edgeRenderer =
        vtkSmartPointer<vtkRenderer>::New();
    edgeRenderer->SetViewport(edgeViewport);

    vtkSmartPointer<vtkRenderWindow> renderWindow =
        vtkSmartPointer<vtkRenderWindow>::New();
    renderWindow->SetSize(600,300);
    renderWindow->SetMultiSamples(0);
    renderWindow->AddRenderer(originalRenderer);
    renderWindow->AddRenderer(edgeRenderer);

    vtkSmartPointer<vtkRenderWindowInteractor> interactor =
        vtkSmartPointer<vtkRenderWindowInteractor>::New();
    interactor->SetRenderWindow(renderWindow);

    std::string inputFilename = filename; //argv[1];

    // Read all the DICOM files in the specified directory.
    vtkSmartPointer<vtkDICOMImageReader> imageIn =
        vtkSmartPointer<vtkDICOMImageReader>::New();
    imageIn->SetFileName(inputFilename.c_str());
    imageIn->Update();

    //vtkSmartPointer<vtkPNGReader> imageIn =
    //  vtkSmartPointer<vtkPNGReader>::New();
    //imageIn->SetFileName(filename.c_str());
    //imageIn->Update();

    vtkSmartPointer<vtkImageActor> imageActor =
        vtkSmartPointer<vtkImageActor>::New();
#if VTK_MAJOR_VERSION <= 5
    imageActor->SetInput(imageIn->GetOutput());
#else
    imageActor->SetInputData(imageIn->GetOutput());
#endif
    originalRenderer->AddActor(imageActor);

    vtkSmartPointer<vtkImageLuminance> il =
        vtkSmartPointer<vtkImageLuminance>::New();
    il->SetInputConnection(imageIn->GetOutputPort());

    vtkSmartPointer<vtkImageCast> ic =
        vtkSmartPointer<vtkImageCast>::New();
    ic->SetOutputScalarTypeToFloat();
    ic->SetInputConnection(il->GetOutputPort());

    // Smooth the image
    vtkSmartPointer<vtkImageGaussianSmooth> gs =
        vtkSmartPointer<vtkImageGaussianSmooth>::New();
    gs->SetInputConnection(ic->GetOutputPort());
    gs->SetDimensionality(2);
    // gs->SetRadiusFactors(.1, .1, 0);
    gs->SetRadiusFactors(1.0, 1.0);
    gs->SetStandardDeviation(1.0);

    // Gradient the image
    vtkSmartPointer<vtkImageGradient> imgGradient =
        vtkSmartPointer<vtkImageGradient>::New();
    imgGradient->SetInputConnection(gs->GetOutputPort());
    imgGradient->SetDimensionality(2);

    vtkSmartPointer<vtkImageMagnitude> imgMagnitude =
        vtkSmartPointer<vtkImageMagnitude>::New();
    imgMagnitude->SetInputConnection(imgGradient->GetOutputPort());

    // Non maximum suppression
    vtkSmartPointer<vtkImageNonMaximumSuppression> nonMax =
        vtkSmartPointer<vtkImageNonMaximumSuppression>::New();
#if VTK_MAJOR_VERSION <= 5
    nonMax->SetMagnitudeInput(imgMagnitude->GetOutput());
    nonMax->SetVectorInput(imgGradient->GetOutput());
#else
    imgMagnitude->Update();
    nonMax->SetMagnitudeInputData(imgMagnitude->GetOutput());
    imgGradient->Update();
    nonMax->SetVectorInputData(imgGradient->GetOutput());
#endif
    nonMax->SetDimensionality(2);

    vtkSmartPointer<vtkImageConstantPad> pad =
        vtkSmartPointer<vtkImageConstantPad>::New();
    pad->SetInputConnection(imgGradient->GetOutputPort());
    pad->SetOutputNumberOfScalarComponents(3);
    pad->SetConstant(0);

    vtkSmartPointer<vtkImageToStructuredPoints> i2sp1 =
        vtkSmartPointer<vtkImageToStructuredPoints>::New();
    i2sp1->SetInputConnection(nonMax->GetOutputPort());
#if VTK_MAJOR_VERSION <= 5
    i2sp1->SetVectorInput(pad->GetOutput());
#else
    pad->Update();
    i2sp1->SetVectorInputData(pad->GetOutput());
#endif

    // Link edgles
    vtkSmartPointer<vtkLinkEdgels> imgLink =
        vtkSmartPointer<vtkLinkEdgels>::New();
    imgLink->SetInputConnection(i2sp1->GetOutputPort());
    double mGradient = 450;
    imgLink->SetGradientThreshold(450);

    // Threshold links
    vtkSmartPointer<vtkThreshold> thresholdEdgels =
        vtkSmartPointer<vtkThreshold>::New();
    thresholdEdgels->SetInputConnection(imgLink->GetOutputPort());
    double mThreshold = 450;
    thresholdEdgels->ThresholdByUpper(mThreshold);
    thresholdEdgels->AllScalarsOff();

    vtkSmartPointer<vtkGeometryFilter> gf =
        vtkSmartPointer<vtkGeometryFilter>::New();
    gf->SetInputConnection(thresholdEdgels->GetOutputPort());

    vtkSmartPointer<vtkImageToStructuredPoints> i2sp =
        vtkSmartPointer<vtkImageToStructuredPoints>::New();
    i2sp->SetInputConnection(imgMagnitude->GetOutputPort());
#if VTK_MAJOR_VERSION <= 5
    i2sp->SetVectorInput(pad->GetOutput());
#else
    pad->Update();
    i2sp->SetVectorInputData(pad->GetOutput());
#endif

    // Subpixel them
    vtkSmartPointer<vtkSubPixelPositionEdgels> spe =
        vtkSmartPointer<vtkSubPixelPositionEdgels>::New();
    spe->SetInputConnection(gf->GetOutputPort());
#if VTK_MAJOR_VERSION <= 5
    spe->SetGradMaps(i2sp->GetStructuredPointsOutput());
#else
    i2sp->Update();
    spe->SetGradMapsData(i2sp->GetStructuredPointsOutput());
#endif


    vtkSmartPointer<vtkStripper> strip =
        vtkSmartPointer<vtkStripper>::New();
    strip->SetInputConnection(spe->GetOutputPort());
    strip->Update();

    FILE *fp = fopen("D:\\2017_Projects\\CannyEdgeDetector\\Data\\Edges.txt","w");
    int nPoints = spe->GetOutput()->GetNumberOfPoints();
    float** edges=new float*[nPoints];
    int i;
    float pixsiz = 1;
    float xht(0), yht(0), zht(0);
    char mode = 'z';
    int file_type = 0, ord(0);;
    for (i=0;i<nPoints;i++)
    {

        edges[i]= new float[3];
        if ( mode == 'z' ) {
            if ( file_type == 0 ) {
                edges[i][0]=((spe->GetOutput()->GetPoint(i))[0])*pixsiz + xht;
                edges[i][1]=((spe->GetOutput()->GetPoint(i))[1])*pixsiz + yht;
            } else {
                edges[i][1]=((spe->GetOutput()->GetPoint(i))[0])*pixsiz + xht;
                edges[i][0]=((spe->GetOutput()->GetPoint(i))[1])*pixsiz + yht;
            }

            edges[i][2] = zht;
        } else if ( mode == 'y' ) {
            if ( file_type == 0 ) {
                edges[i][0]=((spe->GetOutput()->GetPoint(i))[0])*pixsiz + xht;
                edges[i][1]=yht;
            } else {
                edges[i][1]=((spe->GetOutput()->GetPoint(i))[0])*pixsiz + xht;
                edges[i][0]=yht;
            }

            if ( ord == 0 )
                edges[i][2]=((spe->GetOutput()->GetPoint(i))[1])*pixsiz + zht;
            else
                edges[i][2]=-((spe->GetOutput()->GetPoint(i))[1])*pixsiz + zht;

        }

        fprintf(fp, "%f %f %f\n",edges[i][0],  edges[i][1], edges[i][2]); 
    }
    fclose(fp);

    vtkSmartPointer<vtkPolyDataMapper> dsm =
        vtkSmartPointer<vtkPolyDataMapper>::New();
    dsm->SetInputConnection(strip->GetOutputPort());
    dsm->ScalarVisibilityOff();

    vtkSmartPointer<vtkActor> planeActor =
        vtkSmartPointer<vtkActor>::New();
    planeActor->SetMapper(dsm);
    planeActor->GetProperty()->SetAmbient(1.0);
    planeActor->GetProperty()->SetDiffuse(0.0);

    // Add the actors to the renderer, set the background and size
    edgeRenderer->AddActor(planeActor);

    // Render the image
    interactor->Initialize();
    renderWindow->Render();
    renderWindow->Render();

    interactor->Start();

    return EXIT_SUCCESS;
}

