#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Images/PNG_FILE.h>
#include <Core/Parsing/PARSE_ARGS.h>
#include <Core/Read_Write/FILE_UTILITIES.h>
#include <fstream>
#include <iostream>

using namespace PhysBAM;

typedef double T;
typedef VECTOR<T,3> TV;

int approx_max = 5;

int main(int argc,char *argv[])
{
    std::string f0,f1,f2;
    PARSE_ARGS parse_args(argc,argv);
    parse_args.Extra(&f0,"file","first image file");
    parse_args.Extra(&f1,"file","second image file");
    parse_args.Extra(&f2,"file","output image file");
    parse_args.Parse();

    ARRAY<TV,VECTOR<int,2>> img0,img1;
    PNG_FILE<T>::Read(f0,img0);
    PNG_FILE<T>::Read(f1,img1);

    if(img0.domain!=img1.domain)
    {
        printf("DIFF size\n");
        PNG_FILE<T>::Write(f2,img0);
        return 0;
    }

    ARRAY<int> approx(approx_max+1);
    int num_diff=0,total_approx=0;
    
    for(RANGE_ITERATOR<2> it(img0.domain);it.Valid();it.Next())
    {
        int e=rint((img0(it.index)-img1(it.index)).Max_Abs()*255);
        if(e>approx_max)
        {
            num_diff++;
            img0(it.index)=TV(1,0,0);
        }
        else if(e)
        {
            total_approx++;
            approx(e)++;
            img0(it.index)=TV(0,1,0);
        }
        else img0(it.index)/=8;
    }
    PNG_FILE<T>::Write(f2,img0);

    if(!num_diff && !total_approx) printf("SAME\n");
    else
    {
        int n=img0.domain.Size();
        if(num_diff)
            printf("DIFF %i of %i (%8.5f)%s", num_diff, n, (double)num_diff/n,total_approx?"    ":"\n");
        if(total_approx)
        {
            printf("APPROX");
            for(int i=1;i<approx.m;i++)
                printf(" %i", approx(i));
            printf(" -> %i of %i (%8.5f)\n", total_approx, n, (double)total_approx/n);
        }
    }

    return 0;
}
