/* Median v1.0.0 Ver.:OK (c) 12/07 - Christos Iosifides */
/*
   Name:   Median_8UI
   Version:  1.0.0
   Author:  Christos Iosifidis
   Date:   11/12/07 11:03
   Description: Median filtering on 8bit unsigned int BIL datasets
   Usage: >Median_8ui filename_of_the_8bit_unsigned_int_BIL_dataset filename_of_the_8bit_BIL_Medianed_dataset rows columns bands structured_element

   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>

/* http://ndevilla.free.fr/median/median/src/wirth.c - START
 * Algorithm from N. Wirth's book, implementation by N. Devillard.
 * This code in public domain.
 */

typedef unsigned char elem_type;

#define ELEM_SWAP(a,b) { register elem_type t=(a);(a)=(b);(b)=t; }

/*---------------------------------------------------------------------------
   Function :   kth_smallest()
   In       :   array of elements, # of elements in the array, rank k
   Out      :   one element
   Job      :   find the kth smallest element in the array
   Notice   :   use the median() macro defined below to get the median.

                Reference:

                  Author: Wirth, Niklaus
                   Title: Algorithms + data structures = programs
               Publisher: Englewood Cliffs: Prentice-Hall, 1976
    Physical description: 366 p.
                  Series: Prentice-Hall Series in Automatic Computation

 ---------------------------------------------------------------------------*/

elem_type kth_smallest(elem_type a[], int n, int k) {
    register i,j,l,m ;
    register elem_type x ;

    l=0 ;
    m=n-1 ;
    while (l<m) {
        x=a[k] ;
        i=l ;
        j=m ;
        do {
            while (a[i]<x) i++ ;
            while (x<a[j]) j-- ;
            if (i<=j) {
                ELEM_SWAP(a[i],a[j]) ;
                i++ ;
                j-- ;
            }
        } while (i<=j) ;
        if (j<k) l=i ;
        if (k<i) m=j ;
    }
    return a[k] ;
}


#define median(a,n) kth_smallest(a,n,(((n)&1)?((n)/2):(((n)/2)-1)))

/* http://ndevilla.free.fr/median/median/src/wirth.c - END
 * Algorithm from N. Wirth's book, implementation by N. Devillard.
 * This code in public domain.
 */

void error(const char *);

int main(int argc, char *argv[]) {
    FILE *fpin, *fpout;
    int rows, cols, bands, win, win_1, win_2, rows_1, cols_1, cols2, rows_win_2, cols_win_2, cols__win_1, X;
    unsigned char *in_buffer, *out_buffer, *s_element, *a;
    register int r, c, i, j, k;

    /* Check and read argument list items */
    if (argc!=7)
        error("Usage parameters: Input_image Output_image rows cols bands structure_element");

    rows=atoi(argv[3]);
    cols=atoi(argv[4]);
    bands=atoi(argv[5]);
    if (rows<3)
        error("rows must be larger than 3");
    if (cols<3)
        error("columns must be larger than 3");
    if (bands<1)
        error("bands must be more than 1");

    /* Compute helpfull variables */
    cols*=bands;
    rows_1=rows-1;
    cols_1=cols-1;
    cols2=cols<<1;

    /* Open input element file */
    if ((fpin=fopen(argv[6],"r"))==NULL)
        error("Bad input element file name");
    /* Read structuring element size */
    if (fscanf(fpin,"%d",&win)!=1)
        error("Bad element file read 1");
    /* Check structuring element size */
    if (win<3 || win%2==0)
        error("Element size must be odd and greater or equal to 3");
    /* Allocate structuring element memory */
    if ((s_element=(unsigned char *)malloc(win*win))==NULL)
        error("Not enough memory avaliable for structuring element");
    /* Allocate median buffer */
    if ((a=(unsigned char *)malloc(win*win))==NULL)
        error("Not enough memory avaliable for median buffer");
    /* Read structuring element */
    for (i=0;i<win;i++)
        for (j=0;j<win;j++) {
            if (fscanf(fpin,"%d",&X)!=1)
                error("Bad element file read 2");
            s_element[i*win+j]=X?1:0;
        }
    /* Close input element file */
    fclose(fpin);

    win_2=win>>1;
    win_1=win-1;
    rows_win_2=rows-win_2;
    cols_win_2=cols-win_2;
    cols__win_1=cols*win_1;

    /* Allocate input buffer */
    if ((in_buffer=(unsigned char *)malloc(win*cols))==NULL)
        error("Not enough memory avaliable for input buffer");
    /* Open input image */
    if ((fpin=fopen(argv[1],"rb"))==NULL)
        error("Bad input file name");
    /* Open output image */
    if ((fpout=fopen(argv[2],"wb"))==NULL)
        error("Bad output file open");
    /* Allocate output file */
    if ((out_buffer=(unsigned char *)malloc(cols))==NULL)
        error("Not enough memory avaliable for output buffer");

    /* Store 0 in first empty line */
    memset(out_buffer, 0, cols);
    /* Write first win/2 empty lines M_Dilation*/
    for (i=0;i<win_2;i++)
        if (fwrite(out_buffer, sizeof(unsigned char), cols, fpout)!=cols)
            error("Bad write to output file 1");
    /* Read first win-1 lines in buffer */
    if (fread(in_buffer, 1, cols__win_1, fpin)!=cols__win_1)
        error("Bad read from input file 1");

    /* For every row in [ win_2, rows-win_2 ] */
    for (r=win_2;r<rows_win_2;r++) {
        /* Store first and last pixels */
        for (i=0;i<win_2;i++)
            out_buffer[i]=out_buffer[cols-i-1]=0;
        /* Read line in buffer */
        if (fread(in_buffer+cols__win_1, 1, cols, fpin)!=cols)
            error("Bad read from input file 2");
        /* For every column in [ win_2, cols-win_2 ] */
        for (c=win_2;c<cols_win_2;c++) {
            k=0;
            /* For every surrounding pixel in structure element */
            for (i=0;i<win;i++)
                for (j=0;j<win;j++)
                    if (s_element[i*win+j])
                        a[k++]=in_buffer[c-win_2+i*cols+j];
            /* Store result */
            out_buffer[c]=median(a,k);
        }
        /* Scroll cols__win_1 lines up in buffer */
        memcpy(in_buffer, in_buffer+cols, cols__win_1 );
        /* Write results */
        if (fwrite(out_buffer, sizeof(unsigned char), cols, fpout)!=cols)
            error("Bad write to output file 2");
    }
    /* Store 0 in last empty line */
    memset(out_buffer, 0, cols);
    /* Write last win/2 empty lines */
    for (i=0;i<win_2;i++)
        if (fwrite(out_buffer, sizeof(unsigned char), cols, fpout)!=cols)
            error("Bad write to output file 3");

    /* Close files and ... */
    fclose(fpin);
    fclose(fpout);
    free(in_buffer);
    free(out_buffer);
    free(s_element);
    return 0;
}

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