/* geo_pgm v1.0.0 Ver.: OK (c) 03/10 - Christos Iosifides */
/*
   This program resamples a pgm 16bit image to EGSA87 georeferenced bil image

   Name:			geo_pgm
   Version:		    1.0.0
   Copyright:	    Ch Iossif @ 2010
   Author:		    Christos Iosifidis
   Date:			13/03/10 09:24
   Modified:		14/03/10 18:58
   Description:     Resamples a pgm 16bit image to EGSA87 georeferenced bil image
   Compile:
            Linux:     gcc geo_pgm.c -o geo_pgm -lm
            Windows:   gcc geo_pgm.c -o geo_pgm.exe -lm or similar with any ANSI C compiler

    INPUT:
    At command line:
        PGM image file
    At geo_pgm source code file:
      * Check for typical pgm file - MSG Georeference TYPICAL PARAMETERS A
      * Write hdr file - MSG TYPICAL PARAMETERS B
      * Affine polynomial 2nd order - MSG TYPICAL PARAMETERS C

    OUTPUT:
        BIL Image file georeferenced to EGSA'87 and its HDR header

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

void gpl(FILE *stream) {
    fputs("\nCopyright (C) March 2009 Ch Iossif <chiossif@yahoo.com>\n",stream);

    fputs("\nThis program is free software: you can redistribute it and/or modify\n",stream);
    fputs("it under the terms of the GNU General Public License as published by\n",stream);
    fputs("the Free Software Foundation, either version 3 of the License, or\n",stream);
    fputs("(at your option) any later version.\n",stream);

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

    fputs("\nYou should have received a copy of the GNU General Public License\n",stream);
    fputs("along with this program.  If not, see <http://www.gnu.org/licenses/>.\n\n",stream);
}

#define sqr(x) ((x)*(x))
#define abs(x) ((x)>0?(x):(-(x)))

void error(const char *);
double functionforrows(double,double);
double functionforcols(double,double);

int main(int argc, char *argv[]) {
    FILE *fp;
    int row, col, PGM_levels, nrows, ncols, ir, ic;
    double xydim, ulxmap, ulymap, xmap, ymap;
    char buff[BUFSIZ];
    unsigned short int **pgm, sbv;
    register int r, c;

    if (argc!=2) {
        error("Usage: >geo_pgm PGMfilename");
    }
    if ((fp=fopen(argv[1],"rb"))==NULL) {
        sprintf(buff,"Bad filename <%s>.",argv[1]);
        error(buff);
    }
    /* Read first line */
    if (fgets(buff,BUFSIZ,fp)!=buff) {
        sprintf(buff,"Bad read at filename <%s>.",argv[1]);
        error(buff);
    }
    if (buff[0]!='P' && buff[1]!='5') {
        sprintf(buff,"This is not a valid xrit2pic generated PGM file.",argv[1]);
        error(buff);
    }
    /* Skip comment lines */
    do {
        if (fgets(buff,BUFSIZ,fp)!=buff) {
            sprintf(buff,"Bad read at filename <%s>.",argv[1]);
            error(buff);
        }
    } while (buff[0]=='#');
    /* Get row column information */
    if (sscanf(buff,"%d %d", &col, &row)!=2) {
        sprintf(buff,"These are not a valid PGM row column values.",argv[1]);
        error(buff);
    }
    /* Read last line */
    if (fgets(buff,BUFSIZ,fp)!=buff) {
        sprintf(buff,"Bad read at filename <%s>.",argv[1]);
        error(buff);
    }
    if (sscanf(buff,"%d",&PGM_levels)!=1) {
        sprintf(buff,"This is not a valid PGM levels value.",argv[1]);
        error(buff);
    }
    if (row<1 || col<1 || PGM_levels<0 || PGM_levels>32767) {
        sprintf(buff,"Bad data <%d %d %d>.\n", row, col, PGM_levels);
        error(buff);
    }
    /* Check for typical pgm file - MSG Georeference TYPICAL PARAMETERS A */
    if (row!=200)
        error("rows must be 200");
    if (col!=300)
        error("columns must be 300");
    if (PGM_levels<256)
        error("grey levels must be more than 256");
    /* Allocate memory */
    if ((pgm=malloc(row*sizeof(short*)))==NULL)
        error("bad memory allocation");
    for (r=0;r<row;r++)
        if ((pgm[r]=malloc(col*sizeof(short)))==NULL)
            error("bad memory allocation");
    /* Read image and close pgm file*/
    for (r=0;r<row;r++) {
        if (fread(pgm[r],sizeof(short),col,fp)!=col) {
            error("EOF encountered");
        }
    }
    fclose(fp);
    /* Swap'em all */
    for (r=0;r<row;r++)
        for (c=0;c<col;c++) {
            sbv=pgm[r][c];
            pgm[r][c]=( sbv >> 8 ) | ( sbv << 8 );
        }

    /* Write hdr file - MSG TYPICAL PARAMETERS B */
    xydim=3958.591633525;
    ncols=344;
    nrows=264;
    ulxmap=-160516.888595;
    ulymap=4855973.447241;
    sprintf(buff,"%s_egsa87.hdr\0",argv[1]);
    if ((fp=fopen(buff,"w"))==NULL)
        error("bad hdr file open");
    fprintf(fp,"# File created from %s file with geo_pgm (C) March 2009 Ch Iossif <chiossif@yahoo.com>\n",argv[1]);
    fprintf(fp,"byteorder I\n");
    fprintf(fp,"nbits 16\n");
    fprintf(fp,"xdim %lf\n",xydim);
    fprintf(fp,"ydim %lf\n",xydim);
    fprintf(fp,"ncols %d\n",ncols);
    fprintf(fp,"nrows %d\n",nrows);
    fprintf(fp,"nbands 1\n");
    fprintf(fp,"ulxmap %lf\n",ulxmap);
    fprintf(fp,"ulymap %lf\n",ulymap);
    fclose(fp);

    /* Write bil file */
    sprintf(buff,"%s_egsa87.bil\0",argv[1]);
    if ((fp=fopen(buff,"wb"))==NULL)
        error("bad bil file open");
    for (r=0;r<nrows;r++) {
        for (c=0;c<ncols;c++) {
            xmap=ulxmap+c*xydim;
            ymap=ulymap-r*xydim;
            ic=mydoubletoint(functionforrows(xmap,ymap));
            ir=mydoubletoint(functionforcols(xmap,ymap));
            if (ir>=0 && ir<row && ic>=0 && ic<col)
                sbv=pgm[ir][ic];
            else
                sbv=0;
            if (fwrite(&sbv,sizeof(short),1,fp)!=1)
                error("bad bil zero write");
        }
    }
    fclose(fp);

    fputs("\ngeo_pgm v1.0.0 Ver.: OK (c) 03/10 - Christos Iosifides\n",stdout);
    gpl(stdout);
    return 0;
}

void error(const char *s) {
    fprintf(stderr,"geo_pgm Error: %s",s);
    gpl(stderr);
    exit(1);
}

/* Affine polynomial 2nd order - MSG TYPICAL PARAMETERS C */
double functionforrows(double x, double y) {
    return    -2.32003319967266747330551950169798e-12L*sqr(x)+
              3.78083306378852819790610876289173e-12L*sqr(y)+
              3.37743930736601608920838742968450e-04L*x+
              -1.16259996713311052678208118779787e-04L*y+
              -1.16745866967427868818750375420221e-11L*x*y+
              4.33064519309107074401810066888174e+02L;
}

double functionforcols(double x, double y) {
    return     2.68353385332939726747929590319954e-11L*sqr(x)+
               3.07189548167543020857386866282526e-11L*sqr(y)+
               1.12312876171285400134187980606055e-05L*x+
               -4.86947187904294509276991673998522e-04L*y+
               -4.42032505390929251524465798212168e-12L*x*y+
               1.60952746246821588738118435912838e+03L;
}

int mydoubletoint(double x) {
    int k=x;
    if (abs(x-k)>0.5)
        if (k>0)
            k++;
        else
            k--;
    return k;
}