/* MODIS NDVI Change Detector v1.0.0 Ver.: OK (c) 7/07 - Christos Iossifides */
/*
  Name:			MNDVI_CD
  Version:		1.0.0
  Copyright:	Ch Iossif @ 2007
  Author:		Christos Iosifidis
  Date cr.:		03/07/07 09:39
  Date md.:		03/07/07 09:52
  Description:	MODIS NDVI Change Detector

  This code is a normalized Change Detector on two MODIS NDV
  Indexes. The result is float in [-1, 1]. Value -1 also
  indicates No Value.

  Copyright (C) Ch Iossif Jul 2007

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void error(const char *);
void warning(const char *);

int main(int argc,char **argv) {
    FILE *fpold, *fpnew, *fpcd;
    int row, col, g;
    float **oldb, **newb,**cdb, f, f1, f2;
    char buff[BUFSIZ];
    register int r, c;

    if (argc!=6)
        error("Usage: MNDVI_CD OLD NEW DIFF rows cols");
    if ((fpold=fopen(argv[1],"rb"))==NULL)
        error("Bad old band input file name");
    if ((fpnew=fopen(argv[2],"rb"))==NULL)
        error("Bad new band input file name");
    if ((fpcd=fopen(argv[3],"wb"))==NULL)
        error("Bad DIFF band output file name");
    row=atoi(argv[4]);
    if (row<0)
        error("Bad rows value");
    col=atoi(argv[5]);
    if (col<0)
        error("Bad columns value");

    if ((oldb=malloc(row*sizeof(float*)))==NULL)
        error("Not enough memory avaliable for old band");
    if ((newb=malloc(row*sizeof(float*)))==NULL)
        error("Not enough memory avaliable for new band");
    if ((cdb=malloc(row*sizeof(float*)))==NULL)
        error("Not enough memory avaliable for DIFF band");
    for (r=0;r<row;r++){
        if ((oldb[r]=malloc(col*sizeof(float)))==NULL)
            error("Not enough memory avaliable for old band");
        if ((newb[r]=malloc(col*sizeof(float)))==NULL)
            error("Not enough memory avaliable for new band");
        if ((cdb[r]=malloc(col*sizeof(float)))==NULL)
            error("Not enough memory avaliable for DIFF band");
    }

    for (r=0;r<row;r++)
        if (fread(oldb[r], sizeof(float), col, fpold)!=col)
            error("End Of File encountered at old band input file");
    fclose(fpold);
    for (r=0;r<row;r++)
        if (fread(newb[r], sizeof(float), col, fpnew)!=col)
            error("End Of File encountered at new band input file");
    fclose(fpnew);

    for (r=0;r<row;r++)
        for (c=0;c<col;c++){
            if ( (f1=oldb[r][c])>0 && (f2=newb[r][c])>0 )
                cdb[r][c]=(float)((double)f1-f2)/((double)f1+f2);
            else
                cdb[r][c]=(float)-1.0;
        }

    for (r=0;r<row;r++)
        if (fwrite(cdb[r], sizeof(float), col, fpcd)!=col)
            error("Write error at DIFF output file");
    fclose(fpcd);

    sprintf(buff,"%s.pgm", argv[3]);
    if ((fpcd=fopen(buff,"wb"))==NULL)
        error("Bad output filename");
    fprintf(fpcd,"P5 %d %d 255\n", col, row);
    for (r=0;r<row;r++)
        for (c=0;c<col;c++){
            f=cdb[r][c]*255;
            if (f>255)
                g=255;
            else if (f<0)
                g=0;
            else
                g=(int)f;
            if (fputc(g,fpcd)==EOF)
                error("Disk full?");
        }
    fclose(fpcd);
    printf("\n\nMNDVI_CD by Ch Iossif (c) Jul 2007.\n\n\n");

    return 0;
}

void error(const char *s) {
    printf("MNDVI_CD reports error: %s.\n",s);
    exit(1);
}

void warning(const char *s) {
    printf("MNDVI_CD reports: %s.\n",s);
}