/* NRMedian_RGB2I v1.0.3 Ver.:OK (c) 11/08 - Ch Iossif */
/*
   Name:         NRMedian_RGB2I
   Version:      1.0.3
   Author:       Ch Iossif <chiossif@yahoo.com>
   Date:         08/11/08 20:31
   Modified:     09/11/08 09:54
   Description:  RGB to Intensity and Median Noise Reduction on a nxn window basis
   Usage:        $NRMedian_RGB2I tiffimage greys N ex. $NRMedian_RGB2I mytiff.tif 65536 3
   Uses:         libtiff, see <http://directory.fsf.org/project/libtiff>
   Compile line: $gcc NRMedian_RGB2I.c -ltiff

   Copyright (C) November 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 <tiffio.h> /* libtiff */

void error(const char *);
unsigned short int median(unsigned short int *, int);
/*
int compare (const void *, const void *);
*/

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 */
    unsigned int tifftag; /* libtiff */
    char buff[BUFSIZ];
    unsigned char bv_c;
    unsigned short int bv_s, *p, bv_sh, bv_sl, *image, *new_image, *array;
    unsigned int rows, columns, bands, win, win_2, sizeofsample, greys, bv, maxi, mini, red, green, blue;
    register int r, c, b, win_r, win_c;

    if (argc!=4) {
        error("\nUsage:\n\t$NRMedian_RGB2I tiffimage greys N ex. $NRMedian_RGB2I mytiff.tif 65536 3");
    }
    greys=atoi(argv[2]);
    if (greys>65536 || greys<2) {
        error("Greys must be in space [ 2, 65536 ]");
    }
    win=atoi(argv[3]);
    if (win<3 && win%2) {
        error("Window size must be an odd integer greater or equal than 3");
    }
    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 ((image=(unsigned short int *)malloc(rows*columns*sizeof(unsigned short int)))==NULL) {
        error("Bad memory allocation");
    }
    if ((new_image=(unsigned short int *)malloc(rows*columns*sizeof(unsigned short int)))==NULL) {
        error("Bad memory allocation");
    }
    if ((array=(unsigned short int *)malloc(win*win*sizeof(unsigned short int)))==NULL) {
        error("Bad memory allocation");
    }
    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++) {
                red=(unsigned int)p[c*bands];
                green=(unsigned int)p[c*bands+1];
                blue=(unsigned int)p[c*bands+2];
                image[strip*columns+c]=(unsigned short int)((red+green+blue)/3);
            }
        }
    } else {
        error("Sorry: possibly not a uncompressed 16-bit per color tiff image");
    }

    maxi=0;
    mini=65535;
    for (r=0;r<rows;r++) {
        for (c=0;c<columns;c++) {
            bv=image[r*columns+c];
            if (maxi<bv)
                maxi=bv;
            if (mini>bv)
                mini=bv;
        }
    }

    sprintf(buff,"%s.%d.pgm", argv[1], greys);
    if ((fp=fopen(buff,"wb"))==NULL) {
        error("Bad output filename");
    }
    sprintf(buff,"P5\n# Creator: NRMedian_RGB2I v1.0.3 - 11/08 (c) Ch Iossif\n# Source file: %s\n%d %d\n%d\n", argv[1], columns, rows, greys-1);
    if (fputs(buff, fp)==EOF) {
        error("Disk full");
    }
    for (r=0;r<rows;r++) {
        for (c=0;c<columns;c++) {
            bv=image[r*columns+c];
            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=greys*(bv-mini)/(maxi-mini);
                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);

    win_2=win/2;
    for (r = win_2; r < rows-win_2; r++) {
        for (c = win_2; c < columns-win_2; c++) {
            b=0;
            for (win_r = -win_2; win_r <= (int) win_2; win_r++) {
                for (win_c = -win_2; win_c <= (int) win_2; win_c++) {
                    array[b++]=image[((int)r+win_r)*columns+(int)c+win_c];
                }
            }
            new_image[r*columns+c]=median(array, win*win);
        }
    }

    maxi=0;
    mini=65535;
    for (r=0;r<rows;r++) {
        for (c=0;c<columns;c++) {
            bv=new_image[r*columns+c];
            if (maxi<bv)
                maxi=bv;
            if (mini>bv)
                mini=bv;
        }
    }

    sprintf(buff,"%s.%d.NR%d.pgm", argv[1], greys, win);
    if ((fp=fopen(buff,"wb"))==NULL) {
        error("Bad output filename");
    }
    sprintf(buff,"P5\n# Creator: NRMedian_RGB2I v1.0.3 - 11/08 (c) Ch Iossif\n# Source file: %s\n# Window size: %d\n%d %d\n%d\n", argv[1], win, columns, rows, greys-1);
    if (fputs(buff, fp)==EOF) {
        error("Disk full");
    }
    for (r=0;r<rows;r++) {
        for (c=0;c<columns;c++) {
            bv=new_image[r*columns+c];
            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=greys*(bv-mini)/(maxi-mini);
                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);
    free(image);
    free(new_image);
    free(array);
    return 0;
}

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

unsigned short int median(unsigned short int *a, int n) {
    unsigned short int k;
    register int i,j;

    /*  Paradox: qsort compare compination is slower than inline bubble sort */
    /*  qsort (a, n, sizeof(unsigned short int), compare); */

    /* Bubble sort follows */
    for (i=0;i<n-1;i++)
        for (j=i+1;j<n;j++)
            if (a[i]>a[j]) {
                k=a[i];
                a[i]=a[j];
                a[j]=k;
            }

    return a[n/2];
}

/*
int compare (const void * a, const void * b) {
    return ( *(unsigned short int*)a - *(unsigned short int*)b );
}
*/