#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include "stat.h"
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/mman.h>
#include <fcntl.h>
#include <vtkObject.h>
#include "vtkImageData.h"
#include "vtkPointData.h"
#include <vtkDataArray.h>
#include <vtkUnsignedCharArray.h>
#include <vtkDoubleArray.h>
#include <vtkImageData.h>
#include <vtkImageFFT.h>
#include <vtkArrayCalculator.h>
#include <vtkImageExtractComponents.h>

int main (int argc, char *argv[])
{
    double *data;
    int xdim, ydim, zdim;
    double kx;
    int x,y,z;
    long index;
    xdim = 128;
    ydim = 128;
    zdim = 128;

    if ((data = (double*)malloc(sizeof(double)*xdim*ydim*zdim*3)) == NULL) {
        perror("Unable to allocate memory!\n\t");
        return -1;
    }

    printf("Constructing sin(kx * x) dataset...\n");
    kx = 2.0 * 3.1415926 / (double)xdim;
    for(x=0; x<xdim; x++) {
        for(y=0; y<ydim; y++) {
            for(z=0; z<zdim; z++) {
                index = z + (y*zdim) + (x * ydim * zdim);
                data[(index*3)+0] = sin(kx * (double)x);           
                data[(index*3)+1] = cos(kx * (double)x);           
                data[(index*3)+2] = tan(kx * (double)x);           
            }
        }
    }

    vtkDoubleArray *ucPointer = vtkDoubleArray::New();
    ucPointer->SetNumberOfComponents(3);
    ucPointer->SetArray(data, xdim*ydim*zdim*3, 1);
    ucPointer->SetName("image");

    vtkImageData *image = vtkImageData::New();
    image->Initialize();
    image->SetDimensions(xdim, ydim, zdim);
    image->SetExtent(1, xdim, 1, ydim, 1, zdim);
    image->SetScalarTypeToDouble();
    image->SetNumberOfScalarComponents(3);
    image->GetPointData()->SetScalars(ucPointer);
    PrintStatistics(image);
    printf("Extracting Velocity X-Component...\n");

#if 0
    vtkArrayCalculator *extract = vtkArrayCalculator::New();
    extract->SetInput(image);
    extract->SetResultArrayName("result");
    extract->AddScalarVariable("x", "image", 0);
    extract->SetFunction("(x*1.0)+0.0");
    extract->ReplaceInvalidValuesOff();
    extract->Update();
    extract->GetOutput()->GetPointData()->RemoveArray("image");
    extract->GetOutput()->GetPointData()->SetActiveScalars("result");
#else
    vtkImageExtractComponents *extract = vtkImageExtractComponents::New();
    extract->SetInput(image);
    extract->SetComponents(0);
    extract->Update();
#endif
    PrintStatistics(extract->GetOutput());

    printf("Computing FFT...\n");
    vtkImageFFT *fft = vtkImageFFT::New();
    fft->SetInput(extract->GetOutput());
    fft->Update();
    PrintStatistics(fft->GetOutput());

}