/* PCA_true2false v1.0.0 Ver.:OK (c) 4/08 - Ch Iossif */
/*
   Name:         PCA_true2false
   Version:      1.2.0
   Author:       Ch Iossif <chiossif@yahoo.com>
   Date:         12/04/08 23:03
   Modified:     21/04/08 17:40
   Description:  Principal components <http://en.wikipedia.org/wiki/Principal_components_analysis>
                 of a 16bit TIFF color image
   Usage:        $PCA_true2false tiffimage scale triplette ex. $PCA_true2false mytiff.tif 65536 RGB
   Uses: alpha.  libtiff, see <http://directory.fsf.org/project/libtiff>
         beta.   GSL, see <http://www.gnu.org/software/gsl/>
   Compile line: $gcc PCA_true2false.c -ltiff -lgsl -lgslcblas

   Copyright (C) April 2008 Ch Iossif <chiossif@yahoo.com>

   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 3 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, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <tiffio.h> /* libtiff */
#include <gsl/gsl_math.h> /* GSL */
#include <gsl/gsl_eigen.h> /* GSL */

void error(const char *);

int main(int argc, char **argv) {
    FILE *fp;
    TIFF* tif; /* libtiff */
    tdata_t buf; /* libtiff */
    tstrip_t strip; /* libtiff */
    uint32* bc; /* libtiff */
    uint32 stripsize; /* libtiff */
    char buff[BUFSIZ], *pca;
    unsigned char bv_c;
    unsigned short int bv_s, *p, bv_sh, bv_sl;
    unsigned int rows, columns, bands, sizeofsample, tifftag, greys, bv;
    double *s, **sxy, *eval, **evec, *d, seval;
    double **one, mino, maxo, **two, mint, maxt, **three, minr, maxr;
    double **red, redmin, redmax, **green, greenmin, greenmax, **blue, bluemin, bluemax;
    register unsigned int r, c, b, bb;

    if (argc!=4) {
        error("\nUsage:\n\t$PCA_true2false tiffimage scale triplette\nex.\t$PCA_true2false mttiff.tif 65536 RGB");
    }
    greys=atoi(argv[2]);
    if (greys>65536 || greys<2) {
        error("Scale must be in space [ 2, 65536 ]");
    }
    pca=argv[3];
    if (strlen(pca)!=3) {
        error("Choose one of the following triplettes: [ RGB RBG GRB GBR BRG BGR ]");
    }
    if ((tif = TIFFOpen(argv[1], "r"))==NULL) { /* libtiff */
        sprintf(buff,"Bad file name <%s>", argv[1]);
        error(buff);
    }
    TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &tifftag); /* libtiff */
    rows=tifftag;
    TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &tifftag); /* libtiff */
    columns=tifftag;
    TIFFGetField(tif, TIFFTAG_BITSPERSAMPLE, &tifftag); /* libtiff */
    sizeofsample=tifftag>>3;
    TIFFGetField(tif, TIFFTAG_SAMPLESPERPIXEL, &tifftag); /* libtiff */
    bands=tifftag;
    if (bands!=3) {
        error("There must be three samples per pixel [RGB]");
    }
    if ((s=(double *)malloc(bands*sizeof(double)))==NULL) {
        error("Bad memory allocation");
    }
    if ((sxy=(double **)malloc(bands*sizeof(double *)))==NULL) {
        error("Bad memory allocation");
    }
    for (b=0;b<bands;b++) {
        if ((sxy[b]=(double *)malloc(bands*sizeof(double)))==NULL) {
            error("Bad memory allocation");
        }
    }
    for (b=0;b<bands;b++) {
        for (bb=0;bb<bands;bb++) {
            sxy[b][bb]=0;
        }
        s[b]=0;
    }
    TIFFGetField(tif, TIFFTAG_STRIPBYTECOUNTS, &bc); /* libtiff */
    stripsize = bc[0];
    if ((buf = _TIFFmalloc(stripsize))==NULL) { /* libtiff */
        error("Bad memory allocation");
    }
    if (TIFFNumberOfStrips(tif)==rows) { /* libtiff */
        for (strip = 0; strip < rows; strip++) {
            if (bc[strip] > stripsize) {
                buf = _TIFFrealloc(buf, bc[strip]); /* libtiff */
                stripsize = bc[strip];
            }
            p=buf;
            TIFFReadRawStrip(tif, strip, buf, bc[strip]); /* libtiff */
            for (c=0;c<columns;c++)
                for (b=0;b<bands;b++) {
                    s[b]+=(double)p[b+c*bands];
                    for (bb=0;bb<bands;bb++) {
                        sxy[b][bb]+=(double)p[b+c*bands]*p[bb+c*bands];
                    }
                }
        }
    } else {
        error("Sorry: possibly not a uncompressed 16-bit per color tiff image");
    }
    if ((d=(double *)malloc(bands*bands*sizeof(double )))==NULL) {
        error("Bad memory allocation");
    }
    if ((evec=(double **)malloc(bands*sizeof(double *)))==NULL) {
        error("Bad memory allocation");
    }
    if ((eval=(double *)malloc(bands*sizeof(double)))==NULL) {
        error("Bad memory allocation");
    }
    for (b=0;b<bands;b++) {
        if ((evec[b]=(double *)malloc(bands*sizeof(double)))==NULL) {
            error("Bad memory allocation");
        }
    }
    for (b=0;b<bands;b++)
        for (bb=0;bb<bands;bb++)
            d[b*bands+bb]=(sxy[b][bb]-(s[b]*s[bb])/rows/columns)/(rows*columns-1);
    gsl_matrix_view m=gsl_matrix_view_array(d,bands,bands); /* GSL */
    gsl_vector *gsl_eval=gsl_vector_alloc(bands); /* GSL */
    gsl_matrix *gsl_evec=gsl_matrix_alloc(bands,bands); /* GSL */
    gsl_eigen_symmv_workspace *w=gsl_eigen_symmv_alloc(bands); /* GSL */
    gsl_eigen_symmv(&m.matrix,gsl_eval,gsl_evec,w); /* GSL */
    gsl_eigen_symmv_free(w); /* GSL */
    gsl_eigen_symmv_sort(gsl_eval,gsl_evec,GSL_EIGEN_SORT_ABS_DESC); /* GSL */
    for (b=0;b<bands;b++)
        eval[b]=gsl_vector_get(gsl_eval,b); /* GSL */
    for (b=0;b<bands;b++)
        for (bb=0;bb<bands;bb++)
            evec[b][bb]=gsl_matrix_get(gsl_evec,b,bb); /* GSL */
    gsl_vector_free(gsl_eval); /* GSL */
    gsl_matrix_free(gsl_evec); /* GSL */
    mino=1e99;
    maxo=-mino;
    if ((one=(double **)malloc(rows*sizeof(double *)))==NULL) {
        error("Bad memory allocation");
    }
    for (r = strip = 0; r < rows; strip++, r++ ) {
        if ((one[r]=(double *)malloc(columns*sizeof(double)))==NULL) {
            error("Bad memory allocation");
        }
        if (bc[strip] > stripsize) {
            buf = _TIFFrealloc(buf, bc[strip]); /* libtiff */
            stripsize = bc[strip];
        }
        p=buf;
        TIFFReadRawStrip(tif, strip, buf, bc[strip]); /* libtiff */
        for (c=0;c<columns;c++) {
            seval=0.0;
            for (b=0;b<bands;b++)
                seval+=(double)p[b+c*bands]*evec[b][0];
            one[r][c]=seval;
            if (seval>maxo) maxo=seval;
            if (seval<mino) mino=seval;
        }
    }
    mint=1e99;
    maxt=-mint;
    if ((two=(double **)malloc(rows*sizeof(double *)))==NULL) {
        error("Bad memory allocation");
    }
    for (r = strip = 0; r < rows; strip++, r++ ) {
        if ((two[r]=(double *)malloc(columns*sizeof(double)))==NULL) {
            error("Bad memory allocation");
        }
        if (bc[strip] > stripsize) {
            buf = _TIFFrealloc(buf, bc[strip]); /* libtiff */
            stripsize = bc[strip];
        }
        p=buf;
        TIFFReadRawStrip(tif, strip, buf, bc[strip]); /* libtiff */
        for (c=0;c<columns;c++) {
            seval=0.0;
            for (b=0;b<bands;b++)
                seval+=(double)p[b+c*bands]*evec[b][1];
            two[r][c]=seval;
            if (seval>maxt) maxt=seval;
            if (seval<mint) mint=seval;
        }
    }
    minr=1e99;
    maxr=-minr;
    if ((three=(double **)malloc(rows*sizeof(double *)))==NULL) {
        error("Bad memory allocation");
    }
    for (r = strip = 0; r < rows; strip++, r++ ) {
        if ((three[r]=(double *)malloc(columns*sizeof(double)))==NULL) {
            error("Bad memory allocation");
        }
        if (bc[strip] > stripsize) {
            buf = _TIFFrealloc(buf, bc[strip]); /* libtiff */
            stripsize = bc[strip];
        }
        p=buf;
        TIFFReadRawStrip(tif, strip, buf, bc[strip]); /* libtiff */
        for (c=0;c<columns;c++) {
            seval=0.0;
            for (b=0;b<bands;b++)
                seval+=(double)p[b+c*bands]*evec[b][2];
            three[r][c]=seval;
            if (seval>maxr) maxr=seval;
            if (seval<minr) minr=seval;
        }
    }
    _TIFFfree(buf); /* libtiff */
    TIFFClose(tif); /* libtiff */
    sprintf(buff,"%s.%d.%s.ppm", argv[1], greys, pca);
    if ((fp=fopen(buff,"wb"))==NULL) {
        error("Bad output filename");
    }
    sprintf(buff,"P6\n# Creator: PCA_true2false v1.0.0 - 4/08 (c) Ch Iossif\n# Source file: %s\n# Principal Component triplette: %s\n%d %d\n%d\n", argv[1], pca, columns, rows, greys-1);
    if (fputs(buff, fp)==EOF) {
        error("Disk full");
    }
    red=green=blue=NULL;
    if (pca[0]=='r' || pca[0]=='R') red=one, redmax=maxo, redmin=mino;
    else if (pca[0]=='g' || pca[0]=='G') red=two, redmax=maxt, redmin=mint;
    else if (pca[0]=='b' || pca[0]=='B') red=three, redmax=maxr, redmin=minr;
    else error("Bad value at 1st color triplette");
    if (pca[1]=='r' || pca[1]=='R') green=one, greenmax=maxo, greenmin=mino;
    else if (pca[1]=='g' || pca[1]=='G') green=two, greenmax=maxt, greenmin=mint;
    else if (pca[1]=='b' || pca[1]=='B') green=three, greenmax=maxr, greenmin=minr;
    else error("Bad value at 2nd color triplette");
    if (pca[2]=='r' || pca[2]=='R') blue=one, bluemax=maxo, bluemin=mino;
    else if (pca[2]=='g' || pca[2]=='G') blue=two, bluemax=maxt, bluemin=mint;
    else if (pca[2]=='b' || pca[2]=='B') blue=three, bluemax=maxr, bluemin=minr;
    else error("Bad value at 3rd color triplette");
    if (red==NULL||green==NULL||blue==NULL)
        error("Bad values at color triplette. Choose one of the following triplettes: [ RGB RBG GRB GBR BRG BGR ]");
    for (r=0;r<rows;r++)
        for (c=0;c<columns;c++) {
            for (b=0;b<3;b++) {
                if (b==0) {
                    bv=(red[r][c]-redmin)/(redmax-redmin)*greys;
                } else if (b==1) {
                    bv=(green[r][c]-greenmin)/(greenmax-greenmin)*greys;
                } else if (b==2) {
                    bv=(blue[r][c]-bluemin)/(bluemax-bluemin)*greys;
                } else {
                    bv=0;
                }
                if (greys>256) {
                    bv_s=(unsigned short int)(bv>=greys)?greys-1:((bv<0)?0:bv);
                    bv_sh=bv_s&0xff00;
                    bv_sl=bv_s&0x00ff;
                    bv_s=(bv_sh>>8)|(bv_sl<<8);
                    if (fwrite(&bv_s, sizeof(short),1,fp)!=1) {
                        error("Disk full");
                    }
                } else {
                    bv_c=(unsigned char)(bv>=greys)?greys-1:((bv<0)?0:bv);
                    if (fwrite(&bv_c, sizeof(char),1,fp)!=1) {
                        error("Disk full");
                    }
                }
            }
        }
    fclose(fp);
    for (b=0;b<bands;b++) {
        free(sxy[b]);
        free(evec[b]);
    }
    for (r=0;r<rows;r++) {
        free(one[r]);
        free(two[r]);
        free(three[r]);
    }
    free(sxy);
    free(one);
    free(two);
    free(three);
    free(d);
    free(evec);
    free(eval);
    return 0;
}

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