/* MODIS NDVI Calculator v1.0.0 Ver.: OK (c) 7/07 - Christos Iossifides */
/*
  Name:			MNDVI_Calc
  Version:		1.0.0
  Copyright:	Ch Iossif @ 2007
  Author:		Christos Iosifidis
  Date cr.:		02/07/07 17:39
  Date md.:		02/07/07 18:54
  Description:	MODIS NDVI Calculator

  This code calculates the Normalized Difference Vegetation
  Index (NDVI) on MODIS 15-bit bands 1 (red) and 2 (nir).
  NDVI file is a float in [0, 1]. Value -1 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 *fpred, *fpnir, *fpndvi;
    int row, col, g;
    unsigned short **redb, **nirb, dr, dn;
    float **ndvi, f;
    char buff[BUFSIZ];
    register int r, c;

    if (argc!=6)
        error("Usage: MNDVI_Calc RED NIR NDVI rows cols");
    if ((fpred=fopen(argv[1],"rb"))==NULL)
        error("Bad red band input file name");
    if ((fpnir=fopen(argv[2],"rb"))==NULL)
        error("Bad nir band input file name");
    if ((fpndvi=fopen(argv[3],"wb"))==NULL)
        error("Bad NDVI 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 ((redb=malloc(row*sizeof(short*)))==NULL)
        error("Not enough memory avaliable for red band");
    if ((nirb=malloc(row*sizeof(short*)))==NULL)
        error("Not enough memory avaliable for nir band");
    if ((ndvi=malloc(row*sizeof(float*)))==NULL)
        error("Not enough memory avaliable for ndvi band");
    for (r=0;r<row;r++){
        if ((redb[r]=malloc(col*sizeof(short)))==NULL)
            error("Not enough memory avaliable for red band");
        if ((nirb[r]=malloc(col*sizeof(short)))==NULL)
            error("Not enough memory avaliable for nir band");
        if ((ndvi[r]=malloc(col*sizeof(float)))==NULL)
            error("Not enough memory avaliable for ndvi band");
    }

    for (r=0;r<row;r++)
        if (fread(redb[r], sizeof(short), col, fpred)!=col)
            error("End Of File encountered at red band input file");
    fclose(fpred);
    for (r=0;r<row;r++)
        if (fread(nirb[r], sizeof(short), col, fpnir)!=col)
            error("End Of File encountered at nir band input file");
    fclose(fpnir);

    for (r=0;r<row;r++)
        for (c=0;c<col;c++){
            if ((dr=redb[r][c])<32768 && (dn=nirb[r][c])<32768 )
                ndvi[r][c]=(float)(dn-dr)/((float)dn+(float)dr+1);
            else
                ndvi[r][c]=-1.0;
        }

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

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

    return 0;
}

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

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