/*                                                                                                        */
/* Xtractor.c (c) Magnus Palmblad, Leiden University Medical Center 2009-                                 */
/*                                                                                                        */
/* This program is free software; you can redistribute it and/or modify it under the terms of the         */
/* Creative Commons Attribution-Share Alike 3.0 License (http://creativecommons.org/licenses/by-sa/3.0/)  */
/*                                                                                                        */
/* 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.              */
/*                                                                                                        */
/* Contact information: n.m.palmblad@lumc.nl                                                              */
/*                                                                                                        */
/* Xtractor integrates signal in a window around a given m/z and outputs these integrals to standard      */
/* output. This may be useful for binning or peak integration in mass spectrometry profiling.             */
/*                                                                                                        */
/* Usage: ./Xtractor -i <spectrum .xy file> -r <reference_list> -o <output filename> [-b <background>]    */
/*                                                                                                        */
/* All *.xy files in one directory can be processed with the following instructions on the command        */
/* line or in a shell script, like this:                                                                  */
/*                                                                                                        */
/* /*!/bin/sh                                                                                             */
/* for x in *.xy; do                                                                                      */
/* ./Xtractor -i "$x" -r reference_list.ref -o "$x.peaks" -b 5000 -nographics                             */
/* done                                                                                                   */
/*                                                                                                        */
/* compile with e.g. gcc -o Xtractor Xtractor.c -lgd                                                      */
/*                                                                                                        */

#include <stdio.h>
#include <stdlib.h>  
#include <gd.h>
#include <gdfontl.h>

#define MAX_PEAKS 1000

int main(int argc, char *argv[]) 
{
  int nographics=0,included,n_peaks;
  long i,j,spectrum_size;
  char *p,s[100],infile[250],outfile[250],line[250],spectrum_filename[250],reference_filename[250],output_filename[250],peak_name[MAX_PEAKS][100];
  double background=0,*mz,*intensity,peak_mz[MAX_PEAKS],peak_window[MAX_PEAKS],peak_area[MAX_PEAKS];
  FILE *inp,*outp;

  int black,white,red,green,blue,XWIDTH,YWIDTH;
  double intensity_sum=0.0,intensity_max=0.0,xstepsize,ystepsize,ystepsize2,x0,x1,y0,y1,rel_intensity_max=0.0,rel_out_max=0.0;
  gdImagePtr im;
  
  /* parsing command line parameters */
  
  if( (argc==2) && ( (strcmp(argv[1],"--help")==0) || (strcmp(argv[1],"-help")==0) || (strcmp(argv[1],"-h")==0)) ) /* want help? */
    {
      printf("Xtractor - (c) Magnus Palmblad 2009-\n\nusage: Xtractor -i <spectrum .xy file> -r <reference_list> -o <output filename> [-b <background> -nographics]\n");
      return 0;
    }
  
  if (argc<4 || argc>12) /* incorrect number of parameters */
    {
      printf("usage: Xtractor -i <spectrum .xy file> -r <reference_list> -o <output filename> [-b <background> -nographics]\n");
      return -1;
    }
  
  
  /* set default values and parse flags */
  
  for(i=1;i<argc;i++) {
    if( (argv[i][0]=='-') && (argv[i][1]=='i') ) strcpy(spectrum_filename,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( (argv[i][0]=='-') && (argv[i][1]=='r') ) strcpy(reference_filename,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( (argv[i][0]=='-') && (argv[i][1]=='o') ) strcpy(output_filename,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( (argv[i][0]=='-') && (argv[i][1]=='b') ) background=atof(&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( strcmp(argv[i],"-nographics")==0 ) {nographics=1; printf("graphical output off\n");}
  }


  /* read spectrum and count number of entries*/
  
  strcpy(infile,spectrum_filename);
  if ((inp = fopen(infile, "r"))==NULL) {printf("error opening spectrum *.xy file %s",infile); return -1;}
  i=0;
  while (fgets(line, 100, inp) != NULL)
    {
      if (strcmp(line,"\n")==0) continue;
      i++;
    }
  close(inp);  
  spectrum_size=i;

  mz=(double*)malloc(spectrum_size*sizeof(double));
  intensity=(double*)malloc(spectrum_size*sizeof(double));

  
  /* read spectrum and store data */
  
  strcpy(infile,spectrum_filename);
  if ((inp = fopen(infile, "r"))==NULL) {printf("error opening spectrum *.xy file %s",infile); return -1;}
  printf("reading spectrum (%i m/z channels)...",spectrum_size); fflush(stdout); 
  i=0;
  while (fgets(line, 100, inp) != NULL)
    {
      if (strcmp(line,"\n")==0) continue;
      p=strtok(line," ");   
      mz[i]=atof(p); // printf("%f ",mz[i]); fflush(stdout);
      p=strtok('\0'," ");
      intensity[i]=atof(p); // printf("%f\n",intensity[i]); fflush(stdout);
      intensity_sum+=intensity[i]; if (intensity[i]>intensity_max) intensity_max=intensity[i];
      i++;
    }
  close(inp);  
  printf("done (read %i m/z and intensity values)\n",i); fflush(stdout);
  

  /* reading reference file */
  
  strcpy(infile,reference_filename);
  if ((inp = fopen(infile, "r"))==NULL) {printf("error opening reference file %s",reference_filename);return -1;}
  printf("reading reference file %s...",reference_filename); fflush(stdout); 
  i=0;
  fgets(line, 250, inp); /* read reference file header */
  while (fgets(line, 250, inp) != NULL)
    {
      if (strcmp(line,"\n")==0) continue;
      p=strtok(line,"\t");
      strcpy(peak_name[i],p); // printf("%s ",peak_name[i]); fflush(stdout);
      p=strtok('\0',"\t");
      peak_mz[i]=atof(p); // printf("%f ",peak_mz[i]); fflush(stdout);
      p=strtok('\0',"\t");
      p=strtok('\0',"\t");
      peak_window[i]=atof(p); // printf("%f\n",peak_window[i]); fflush(stdout);
      peak_area[i]=0;
      i++; 
    }
  close(inp);
  n_peaks=i;
  printf("done (read reference file %s)\n",reference_filename); fflush(stdout); 


  /* opening output file */
  
  strcpy(outfile,output_filename);
  if ((outp = fopen(outfile, "w"))==NULL) {printf("error opening output file %s",output_filename);return -1;}
  printf("writing peak areas to %s...",output_filename); fflush(stdout); 
  
  for(i=0;i<spectrum_size;i++)
    {
      for(j=0;j<n_peaks;j++)
	{
	  if(fabs(mz[i]-peak_mz[j])<=peak_window[j]) peak_area[j]+=intensity[i]-background;
	}
    }

  for(j=0;j<n_peaks;j++) fprintf(outp,"%s %f %f %f\n",peak_name[j],peak_mz[j],peak_window[j],peak_area[j]);

  close(outp);
  printf("done (wrote %i peak areas)\n",n_peaks); fflush(stdout); 
  
  if(!nographics)
    {
      strcpy(outfile,output_filename);
      strcat(outfile,".png");
      if ((outp = fopen(outfile, "w"))==NULL) {printf("error opening output image file %s",outfile);return -1;}
      printf("making spectrum image..."); fflush(stdout); 
      
      XWIDTH=(int)10*ceil(mz[spectrum_size-1]-mz[0])+20; // printf("%i\n",XWIDTH); /* 10 points per m/z unit */
      YWIDTH=500;

      im = gdImageCreate(XWIDTH, YWIDTH); 
      white = gdImageColorAllocate(im, 255, 255, 255);  
      black = gdImageColorAllocate(im, 0, 0, 0);  
      red = gdImageColorAllocate(im, 255, 0, 0);  
      green = gdImageColorAllocate(im, 0, 255, 0);  
      blue = gdImageColorAllocate(im, 0, 0, 255);  
	
      gdImageLine(im, 10, YWIDTH-20, XWIDTH-10, YWIDTH-20, black); /* x-axis */
      sprintf(s,"m/z");
      gdImageString(im,gdFontGetSmall(),XWIDTH-20-2.75*strlen(s),YWIDTH-13,s,black);
      sprintf(s,"spectrum %s",spectrum_filename);
      gdImageString(im,gdFontGetSmall(),20,10,s,black);
	
      x0=mz[0]; xstepsize=(XWIDTH-20)/(mz[spectrum_size-1]-mz[0]);
      while (x0<mz[spectrum_size-1])
	{
	  x1=10+(x0-mz[0])*xstepsize;
	  if((x0-floor(x0))<((x0-0.01)-floor(x0-0.01))) 
	    {
	      gdImageLine(im,x1,YWIDTH-20,x1,YWIDTH-15,black); /* make a tick mark for every integer m/z */
	      if(((int)floor(x0)%10)==0) /* print m/z labels under axis */
		{
		  sprintf(s,"%i",(int)floor(x0));
		  if( (x1-2.75*strlen(s))<(XWIDTH-50) ) gdImageString(im,gdFontGetSmall(),x1-2.75*strlen(s),YWIDTH-13,s,black);
		}
	    }
	  x0+=0.01;
	}
      
      ystepsize=(YWIDTH-50)/intensity_max; /* scale to highest peak */
      
      for(i=0;i<spectrum_size-1;i++) 
	{
	  x0=10+(mz[i]-mz[0])*xstepsize; x1=10+(mz[i+1]-mz[0])*xstepsize; 
	  y0=YWIDTH-20-intensity[i]*ystepsize; y1=YWIDTH-20-intensity[i+1]*ystepsize;
	  if((y0>0)&&(y1>0)) gdImageLine(im,x0,y0,x1,y1,green);
	}
      
      for(i=0;i<spectrum_size-1;i++) 
	{
	  x0=10+(mz[i]-mz[0])*xstepsize; x1=10+(mz[i+1]-mz[0])*xstepsize; 
	  y0=YWIDTH-20-intensity[i]*ystepsize; y1=YWIDTH-20-intensity[i+1]*ystepsize;
	  if((y0>0)&&(y1>0)) 
	    {
	      included=0;
	      for(j=0;j<n_peaks;j++)
		{
		  if(fabs(mz[i]-peak_mz[j])<=peak_window[j]) included=1;
		}
	      if (included) gdImageLine(im,x0,y0,x1,y1,red);
	    }
	}

      if((outp=fopen(outfile, "wb"))==NULL) printf("could not write PNG image %s\n",outfile);  
      gdImagePng(im, outp);
      fclose(outp);
      gdImageDestroy(im);
      printf("done (wrote image to %s)\n",outfile); fflush(stdout); 
    }
  
  return 0;
}