// Computes median intensity image of any number of possibly resized input files
// Quicksort from the Algorithms wiki at http://www.algorithm-code.com/wiki/Quick_Sort
// Compile: gcc median_image.c -o median_image -lgd

#include <stdio.h> 
#include <gd.h> 
#include <stdlib.h>  
#include <string.h>  
#define XSIZE 640 
#define YSIZE 480 

void swap(int l[], int i, int j)
{  int dummy;
   dummy=l[j];
   l[j]=l[i];
   l[i]=dummy;  
 }

int partions(int l[],int low,int high)
{
    int prvotkey=l[low];
    while (low<high)
    {
        while (low<high && l[high]>=prvotkey)
            --high;
        swap(l,high,low);
        while (low<high && l[low]<=prvotkey) 
            ++low;
        swap(l,high,low);      
    }

    return low;
}

void qsorter(int l[],int low,int high)
{
    int prvotloc;
    if(low<high)
    {
        prvotloc=partions(l,low,high);    
        qsorter(l,low,prvotloc); 
        qsorter(l,prvotloc+1,high); 
    }
}

void quicksort(int l[],int n)
{
    qsorter(l,0,n); 
}

int main(int argc, char *argv[]) {
 gdImagePtr *images,median_image,img;
 int x,y,c,i,nb_images,median,*intensities;
 FILE *fp;
 if (argc<2) {
  printf("Syntax: median_image <image1.jpg> <image2.jpg> ... <out.png>\n");
  exit(1);
 }
 nb_images=argc-2;
 printf("Number of images: %d\n",nb_images);
 fp = fopen(argv[argc-1], "rb");
 if(fp != NULL) {
	printf("ERROR: %s already exists.\n",argv[argc-1]);
	fclose(fp);
	exit(1);
 }
 images = malloc(sizeof(gdImagePtr) * nb_images);
 for (i=1;i<argc-1;i++) {
	fp = fopen(argv[i], "rb");
	if(fp == NULL) {
		printf("ERROR: Unable to open %s\n",argv[i]);
		exit(1);
	}
	img=gdImageCreateFromJpeg(fp);
	images[i-1]=gdImageCreateTrueColor(XSIZE,YSIZE);
	gdImageCopyResized(images[i-1], img, 0, 0, 0, 0, XSIZE, YSIZE, img->sx, img->sy); 
	fclose(fp);
	if(images[i-1] == NULL) {
		printf("ERROR: %s is no JPEG file.\n",argv[i]);
		exit(1);
	}
	printf("%s loaded.\n",argv[i]);
 }
 intensities=malloc(sizeof(int) * nb_images);
 median_image=gdImageCreateTrueColor(images[0]->sx,images[0]->sy);
 for (x=0;x<images[0]->sx;x++) {
  for (y=0;y<images[0]->sy;y++) {
   for (i=0;i<nb_images;i++) {
    c=gdImageGetPixel(images[i],x,y);
    intensities[i]=(gdTrueColorGetRed(c)+gdTrueColorGetGreen(c)+gdTrueColorGetBlue(c))/3;
   }
   quicksort(intensities,nb_images);
   if (nb_images % 2 == 0) median=(intensities[(nb_images/2)-1]+intensities[nb_images/2])/2;
   else median=intensities[nb_images/2];
   c = gdImageColorAllocate(median_image, median, median, median); 
   gdImageSetPixel(median_image, x, y, c);
  } 
 }
 for (i=0;i<nb_images;i++) gdImageDestroy(images[i]);
 fp = fopen(argv[argc-1], "wb");
 gdImagePng(median_image, fp);
 fclose(fp);
 printf("Saved to %s\n",argv[argc-1]);
 gdImageDestroy(median_image);
 free(intensities);
 free(images);
 return 0;
}