/*                                                                                                 */                     
/*  mzXMLplot - plot mzXML files as "heat maps" or "virtual gels"                                  */
/*                                                                                                 */
/* (c) Magnus Palmblad, The University of Reading, 2007-                                           */ 
/*                                                                                                 */
/* 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: magnus.palmblad@gmail.com                                                  */
/*                                                                                                 */
/*                                                                                                 */
/* compile with e.g. gcc -o mzXMLplot mzXMLplot.c base64.c ramp.c -I. -lgd -lm -lz                 */
/*                                                                                                 */
/*                                                                                                 */

#include <stdio.h>
#include <stdlib.h>  
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <gd.h>
#include <gdfontl.h>
#include "ramp.h"
#include "base64.h"



/* mzXML read functions adopted from Pep3D */

typedef struct {
  int size;
  double * xval;
  double * yval;
} spectStrct;

void freeSpectStrct(spectStrct spectrum)
{
  free(spectrum.xval);
  free(spectrum.yval);
  
  return;
}

ramp_fileoffset_t *pScanIndex;
int iLastScan;
struct ScanHeaderStruct scanHeader;
struct RunHeaderStruct runHeader;
ramp_fileoffset_t indexOffset;
spectStrct *spects;
double *wghs;
int spectNum;
long i;
RAMPREAL *pPeaks;
int n, n_MS_spectra=0;
RAMPFILE *pFI;

int getMsSpect(spectStrct *msSpect, RAMPFILE *pFI, int scanNum[2])
{
  void getCmbSpect(spectStrct *cmbSpect, int spectNum, spectStrct *spects, double *wghs);

  msSpect->size = -1;
  scanNum[0]=scanNum[0]>1?scanNum[0]:1;
  scanNum[1]=scanNum[1]<iLastScan?scanNum[1]:iLastScan;
  spectNum=scanNum[1]-scanNum[0]+1;
  if(spectNum < 1){
    printf("invalid scan number: %d-%d (full scan range: 1-%d)\n",scanNum[0], scanNum[1], iLastScan);
    fflush(stdout);
    free (pScanIndex);
    return -1;    
  }

  spects=(spectStrct *) calloc(spectNum, sizeof(spectStrct));
  spectNum=0;
  for(i=scanNum[0];i<=scanNum[1];++i) 
    {
      if((scanHeader.msLevel==1)&&(scanHeader.peaksCount>0)) /* MS ? */
	{                         
	  spects[spectNum].size=scanHeader.peaksCount;
	  spects[spectNum].xval=(double *) calloc(spects[spectNum].size, sizeof(double));
	  spects[spectNum].yval=(double *) calloc(spects[spectNum].size, sizeof(double));
	  
	  pPeaks=readPeaks(pFI,pScanIndex[i]);
	  
	  spects[spectNum].size=0;
	  n = 0;
	  while(pPeaks[n]!=-1)
	    {
	      spects[spectNum].xval[spects[spectNum].size]=pPeaks[n];
	      n++;
	      spects[spectNum].yval[spects[spectNum].size]=pPeaks[n];
	      n++;
	      ++(spects[spectNum].size);
	    }
	  free (pPeaks);
	  n_MS_spectra++;
	  if(spects[spectNum].size>0) ++spectNum; 
	  else freeSpectStrct(spects[spectNum]);
	}
      
    } 
  
  if(spectNum>0) 
    {
    wghs= (double *) calloc(spectNum, sizeof(double));
    for (i=0;i<spectNum;++i)
      wghs[i]=1.;
    getCmbSpect(msSpect, spectNum, spects, wghs);
    free(wghs);
  }
  else 
    {
      printf("cannot find an MS spectrum..."); fflush(stdout);
      for(i=0;i<spectNum;++i) freeSpectStrct(spects[i]);
      free(spects);
      return -1;
    }
  
  for(i=0;i<spectNum;++i) freeSpectStrct(spects[i]);
  free(spects);
  
  return 0;
}

void getCmbSpect(spectStrct *cmbSpect, int spectNum, spectStrct *spects, double *wghs)
{
  void copySpectStrct(spectStrct * tgtSpect, spectStrct srcSpect);
  spectStrct tmpSpect[2];
  int indx, indx1, indx2;
  double tmpWghs[2]={1.,1.};
  int i;
  
  if(spectNum<1) return;
  if(spectNum==1) 
    {
      copySpectStrct(cmbSpect,spects[0]);
      if(wghs[0]!=1.) for(i=0;i<cmbSpect->size;++i) cmbSpect->yval[i]*=wghs[0];
      return;
    } 
  
  return;
}

void copySpectStrct(spectStrct * tgtSpect, spectStrct srcSpect)
{
  int i;

  tgtSpect->size=srcSpect.size;
  tgtSpect->xval=(double *) calloc(tgtSpect->size, sizeof(double));
  tgtSpect->yval=(double *) calloc(tgtSpect->size, sizeof(double));

  for(i=0;i<tgtSpect->size;++i) 
    {
    tgtSpect->xval[i]=srcSpect.xval[i];
    tgtSpect->yval[i]=srcSpect.yval[i];
    }

  return;
}



/* main program starts here */

int main(int argc, char *argv[]) 
{
  FILE *inp, *outp;
  spectStrct mzXML_spectrum, temp_spectrum;
  RAMPFILE *mzXML_file;
  char *p, mzXML_filename[100], image_filename[100], scan_range[100], mass_range[100], image_type_string[100], pixels[100], s[200];
  int range[2], image_type;
  ramp_fileoffset_t offset, *scan_index;
  struct ScanHeaderStruct scan_header;
  struct RunHeaderStruct run_header;

  int **gd_data;
  int gray[256], colormap[256][3], *scan_number;
  int black,white,red,green,blue,XWIDTH,YWIDTH, color_scheme;
  gdImagePtr im;
  
  int mzXML_MS_scan, scan, first_scan, last_scan, n_scans, x, y, c, lowest_mass, highest_mass, mz_range, horizontal_pixels, vertical_pixels, n_horizontal_labels, n_vertical_labels, horizontal_label_interval, vertical_label_interval, first_horizontal_label, first_vertical_label, MS_scans, noaxes;
  long i,j,k;
  float min_yval, max_yval, log_min_yval, log_max_yval, power, mz_per_pixel, scans_per_pixel, background, saturation;

  float jet[64][3]={{0,0,0.5625}, /* OK - this is ugly, I know... */
		    {0,0,0.625},
		    {0,0,0.6875},
		    {0,0,0.75},
		    {0,0,0.8125},
		    {0,0,0.875},
		    {0,0,0.9375},
		    {0,0,1},
		    {0,0.0625,1},
		    {0,0.125,1},
		    {0,0.1875,1},
		    {0,0.25,1},
		    {0,0.3125,1},
		    {0,0.375,1},
		    {0,0.4375,1},
		    {0,0.5,1},
		    {0,0.5625,1},
		    {0,0.625,1},
		    {0,0.6875,1},
		    {0,0.75,1},
		    {0,0.8125,1},
		    {0,0.875,1},
		    {0,0.9375,1},
		    {0,1,1},
		    {0.0625,1,1},
		    {0.125,1,0.9375},
		    {0.1875,1,0.875},
		    {0.25,1,0.8125},
		    {0.3125,1,0.75},
		    {0.375,1,0.6875},
		    {0.4375,1,0.625},
		    {0.5,1,0.5625},
		    {0.5625,1,0.5},
		    {0.625,1,0.4375},
		    {0.6875,1,0.375},
		    {0.75,1,0.3125},
		    {0.8125,1,0.25},
		    {0.875,1,0.1875},
		    {0.9375,1,0.125},
		    {1,1,0.0625},
		    {1,1,0},
		    {1,0.9375,0},
		    {1,0.875,0},
		    {1,0.8125,0},
		    {1,0.75,0},
		    {1,0.6875,0},
		    {1,0.625,0},
		    {1,0.5625,0},
		    {1,0.5,0},
		    {1,0.4375,0},
		    {1,0.375,0},
		    {1,0.3125,0},
		    {1,0.25,0},
		    {1,0.1875,0},
		    {1,0.125,0},
		    {1,0.0625,0},
		    {1,0,0},
		    {0.9375,0,0},
		    {0.875,0,0},
		    {0.8125,0,0},
		    {0.75,0,0},
		    {0.6875,0,0},
		    {0.625,0,0},
		    {0.5625,0,0}};

 
  /* 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("mzXMLplot - (c) Magnus Palmblad 2007-\n\nusage: mzXMLplot -i<mzXML file> -o<image filename> -s<first scan>,<last scan> -m<lowest mass>,<highest mass> [-t<image type> -c<color scheme> -p<horizontal pixels>,<vertical pixels> -B<background (absolute value)> -S<saturation (as fraction of maximum)>].\nThe <image type> is l or L for log scale, 1 for linear scale, 2 for square root, 3 for 3rd root etc. (up to 30000). Use the negation of these values for an inverted map!\n\nfor more information, see http://www.ms-utils.org/mzXMLplot or e-mail magnus.palmblad@gmail.com\n");
      return 0;
    }
  
  if (argc<5 || argc>22) /* incorrect number of parameters */
    {
      printf("usage: mzXMLplot -i<mzXML file> -o<image filename> -s<first scan>,<last scan> -m<lowest mass>,<highest mass> [-t<image type> -p<horizontal pixels>,<vertical pixels> -B<background (absolute value)> -S<saturation (as fraction of maximum)>] (type mzXMLplot --help for more information)\n");
      return -1;
    }
  
  
  /* set default values and parse flags */

  image_type=1; color_scheme=0; background=0.0; saturation=1.0; noaxes=0; /* linear, grayscale, m/z on horizontal axis */
  
  for(i=1;i<argc;i++) {
    if( (argv[i][0]=='-') && (argv[i][1]=='i') ) strcpy(mzXML_filename,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( (argv[i][0]=='-') && (argv[i][1]=='o') ) strcpy(image_filename,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( (argv[i][0]=='-') && (argv[i][1]=='s') ) 
      {
	strcpy(scan_range,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]); p=strtok(scan_range,","); 
	first_scan=atoi(p); p=strtok('\0',","); last_scan=atoi(p);
      }
    if( (argv[i][0]=='-') && (argv[i][1]=='m') ) 
      {
	strcpy(mass_range,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]); p=strtok(mass_range,","); 
	lowest_mass=atoi(p); p=strtok('\0',","); highest_mass=atoi(p);
      }
    if( (argv[i][0]=='-') && (argv[i][1]=='t') ) 
      {
	strcpy(image_type_string,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
	if( (strchr(image_type_string,'l')!=NULL) || (strchr(image_type_string,'L')!=NULL) )
	  if( (strchr(image_type_string,'-')!=NULL) ) image_type=-31000; else image_type=31000;
	else
	  image_type=atoi(image_type_string);
      }
    if( (argv[i][0]=='-') && (argv[i][1]=='c') ) color_scheme=atoi(&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( (argv[i][0]=='-') && (argv[i][1]=='p') ) 
      {
	strcpy(pixels,&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]); p=strtok(pixels,","); 
	horizontal_pixels=atoi(p); p=strtok('\0',","); vertical_pixels=atoi(p);
      }
    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( (argv[i][0]=='-') && (argv[i][1]=='S') ) saturation=atof(&argv[strlen(argv[i])>2?i:i+1][strlen(argv[i])>2?2:0]);
    if( strcmp(argv[i],"-noaxes")==0 ) noaxes=1;
  }

  XWIDTH=horizontal_pixels+(1-noaxes)*40; YWIDTH=vertical_pixels+(1-noaxes)*20; /* add pixels for axes */

  if( (abs(image_type)!=31000) && (abs(image_type)>30000) ) {printf("maximum root is 30000\n"); fflush(stdout); return -1;}

  printf("checking mzXML dataset %s...",mzXML_filename); fflush(stdout);
  mzXML_file=rampOpenFile(mzXML_filename);
  indexOffset = getIndexOffset (mzXML_file); /* read the index offset */
  pScanIndex = readIndex(mzXML_file, indexOffset, &iLastScan); /* read the scan index into a vector and get LastScan */
  readRunHeader(mzXML_file, pScanIndex, &runHeader, iLastScan);
  printf("done\n"); fflush(stdout);
  
  printf("allocating memory..."); fflush(stdout);
  
  n_scans=last_scan-first_scan+1;
  gd_data=malloc((n_scans>vertical_pixels?n_scans:vertical_pixels)*sizeof(int *));
  if(gd_data==NULL) {fprintf(stderr,"out of memory (terminating)\n"); return -1;}
  scan_number=malloc(n_scans*sizeof(int));
  for(y=0;y<(n_scans>vertical_pixels?n_scans:vertical_pixels);y++) 
    {
      gd_data[y]=malloc(horizontal_pixels*sizeof(int));
      if(gd_data[y]==NULL) {fprintf(stderr,"out of memory (terminating)\n"); return -1;}
    }
  printf("done\nscanning mzXML data for yval range..."); fflush(stdout);

  /* get min_yval and max_yval */

  min_yval=1E10; max_yval=0; y=-1;
  for(scan=first_scan;scan<last_scan;scan++) 
    {
      range[0]=scan;
      readHeader(mzXML_file, pScanIndex[range[0]], &scanHeader);
      if(scanHeader.msLevel==1) /* MS spectrum */
	{
	  y++;
	  scan_number[y]=scan; /* store scan numbers for axis */
	  range[1]=range[0];
	  if(getMsSpect(&mzXML_spectrum,mzXML_file,range)<0) break; /* MS spectrum not found in mzXML */
	  for(i=0;i<mzXML_spectrum.size;i++) /* # of peaks in peaklist in line spectra */
	    {
	      if(mzXML_spectrum.yval[i]>min_yval) min_yval=mzXML_spectrum.yval[i];
	      if(mzXML_spectrum.yval[i]>max_yval) max_yval=mzXML_spectrum.yval[i];
	    }
	}
    }

  printf("done\nparsing mzXML data..."); fflush(stdout);

  if(min_yval>0) log_min_yval=log(min_yval); log_max_yval=log(max_yval);

  /* read in spectrum by spectrum, store in temporary array */
      
  mz_per_pixel=(float)(highest_mass-lowest_mass)/(float)horizontal_pixels; 

  for(y=0;y<vertical_pixels;y++) for(x=0;x<horizontal_pixels;x++) gd_data[y][x]=0; /* initialize gd_data */
    
  y=-1; 

  if(abs(image_type)==1)  /* linear */
    {
      for(scan=first_scan;scan<last_scan;scan++) 
	{
	  range[0]=scan;
	  readHeader(mzXML_file, pScanIndex[range[0]], &scanHeader);
	  if(scanHeader.msLevel==1) /* MS spectrum */
	    {
	      y++;
	      range[1]=range[0];
	      if(getMsSpect(&mzXML_spectrum,mzXML_file,range)<0) break; /* MS spectrum not found in mzXML */
	      for(i=0;i<mzXML_spectrum.size;i++) /* # of peaks in peaklist in line spectra */
		{
		  if(mzXML_spectrum.xval[i]<lowest_mass) continue;
		  if(mzXML_spectrum.xval[i]>highest_mass) break;
		  if(mzXML_spectrum.yval[i]>background) c=round(255*(mzXML_spectrum.yval[i]-background)/(max_yval*saturation));
		  else c=0;
		  if(c>255) c=255; /* saturated */
		  x=(int)round((mzXML_spectrum.xval[i]-lowest_mass)/mz_per_pixel);
		  if(c>gd_data[y][x]) gd_data[y][x]=c; /* store maximum in each pixel interval */
		}
	    }
	}
    } 
  
  if(abs(image_type)==31000)  /* log */
    {
      for(scan=first_scan;scan<last_scan;scan++) 
	{
	  range[0]=scan;
	  readHeader(mzXML_file, pScanIndex[range[0]], &scanHeader);
	  if(scanHeader.msLevel==1) /* MS spectrum */
	    {
	      y++;
	      range[1]=range[0];
	      if(getMsSpect(&mzXML_spectrum,mzXML_file,range)<0) break; /* MS spectrum not found in mzXML */
	      for(i=0;i<mzXML_spectrum.size;i++) /* # of peaks in peaklist in line spectra */
		{
		  if(mzXML_spectrum.xval[i]<lowest_mass) continue;
		  if(mzXML_spectrum.xval[i]>highest_mass) break;
		  if(mzXML_spectrum.yval[i]>background) c=round(255*log((mzXML_spectrum.yval[i]-background))/(log_max_yval*saturation));
		  else c=0;
		  if(c>255) c=255; /* saturated */
		  x=(int)round((mzXML_spectrum.xval[i]-lowest_mass)/mz_per_pixel);
		  if(c>gd_data[y][x]) gd_data[y][x]=c; /* store maximum in each pixel interval */
		}
	    }
	}
    }
  else 
    if(abs(image_type)>1)  /* nth root */
      {
	power=(float)1/(float)abs(image_type);
	for(scan=first_scan;scan<last_scan;scan++) 
	  {
	    range[0]=scan;
	    readHeader(mzXML_file, pScanIndex[range[0]], &scanHeader);
	    if(scanHeader.msLevel==1) /* MS spectrum */
	      {
		y++; 
		range[1]=range[0];
	      if(getMsSpect(&mzXML_spectrum,mzXML_file,range)<0) break; /* MS spectrum not found in mzXML */
	      for(i=0;i<mzXML_spectrum.size;i++) /* # of peaks in peaklist in line spectra */
		{
		  if(mzXML_spectrum.xval[i]<lowest_mass) continue;
		  if(mzXML_spectrum.xval[i]>highest_mass) break;
		  if(mzXML_spectrum.yval[i]>background) c=round(255*pow((mzXML_spectrum.yval[i]-background),power)/pow(max_yval*saturation,power));
		  else c=0;
		  if(c>255) c=255; /* saturated */
		  x=(int)round((mzXML_spectrum.xval[i]-lowest_mass)/mz_per_pixel); 
		  if(c>gd_data[y][x]) gd_data[y][x]=c; /* store maximum in each pixel interval */
		}
	      }
	  }
      }
  
 
  printf("done\n"); fflush(stdout);

  MS_scans=y; scans_per_pixel=(float)MS_scans/(float)vertical_pixels;
   
  if(abs(image_type)!=31000) printf("generating %ix%i image of type %i color scheme %i from mzXML...",horizontal_pixels,vertical_pixels,image_type,color_scheme);
  else 
    {
      if(image_type<0) printf("generating image of type -log from mzXML..."); else printf("generating image of type log from mzXML...");
    }
  fflush(stdout);
 
  if ((outp=fopen(image_filename, "wb"))==NULL) {printf("\nerror opening output file %s\n",image_filename); return -1;}

  im = gdImageCreate(XWIDTH, YWIDTH);  
  white = gdImageColorAllocate(im, 255, 255, 255);
  black = gdImageColorAllocate(im, 0, 0, 0);  

  /* grayscale color scheme by default */
  
  if(color_scheme==0)
    {
      if(image_type>0) for(i=0;i<=255;i++) gray[i]=gdImageColorAllocate(im, i, i, i); 
      if(image_type<0) for(i=0;i<=255;i++) gray[i]=gdImageColorAllocate(im, 255-i, 255-i, 255-i); 
    }
  
  if(color_scheme==1) /* MATLAB "hot" map */
    {
      for(i=0;i<=255;i++)
	{
	  colormap[i][0]=i>85?255:i*3; 
	  colormap[i][1]=i<86?0:(i-85)*3; if(colormap[i][1]>255) colormap[i][1]=255;
	  colormap[i][2]=i<171?0:(i-170)*3;
	}
      if(image_type>0) for(i=0;i<=255;i++) gray[i]=gdImageColorAllocate(im,colormap[i][0],colormap[i][1],colormap[i][2]); 
      if(image_type<0) for(i=0;i<=255;i++) gray[i]=gdImageColorAllocate(im,255-colormap[i][0],255-colormap[i][1],255-colormap[i][2]);
    }
  if(color_scheme==2) /* MATLAB "jet" map */
    {
      for(i=0;i<=255;i++)
	{
	  colormap[i][0]=(int)round(jet[(int)round(((float)i)/4.0)][0]*255); 
	  colormap[i][1]=(int)round(jet[(int)round(((float)i)/4.0)][1]*255);
	  colormap[i][2]=(int)round(jet[(int)round(((float)i)/4.0)][2]*255);
	}
      if(image_type>0) for(i=0;i<=255;i++) gray[i]=gdImageColorAllocate(im,colormap[i][0],colormap[i][1],colormap[i][2]); 
      if(image_type<0) for(i=0;i<=255;i++) gray[i]=gdImageColorAllocate(im,255-colormap[i][0],255-colormap[i][1],255-colormap[i][2]);
    }

  

  



  /* plot data */

  if(scans_per_pixel>=1)
    {
      for(y=0;y<MS_scans;y++)
	{
	  for(x=0;x<horizontal_pixels;x++) gdImageSetPixel(im,x+(1-noaxes)*40,YWIDTH-(1-noaxes)*20-(int)round((float)y/(float)scans_per_pixel),gray[gd_data[y][x]]);
	}
    }
  else
    {
      for(y=0;y<MS_scans;y++)
	{
	  for(i=0;i<=ceil(1/scans_per_pixel);i++)
	    {
	      for(x=0;x<horizontal_pixels;x++) gdImageSetPixel(im,x+(1-noaxes)*40,YWIDTH-(1-noaxes)*20-i-(int)round((float)y/(float)scans_per_pixel),gray[gd_data[y][x]]);
	    }
	}
    }

  if(noaxes==1) printf("done "); 
  else
    {
      printf("done (%i MS scans)\ndrawing axes...",MS_scans); fflush(stdout);
     
      
      gdImageLine(im, 39, YWIDTH-19, XWIDTH, YWIDTH-19, black); /* x-axis */
      gdImageLine(im, 39, YWIDTH-19, 39, 0, black); /* y-axis */
      sprintf(s,"m/z");
      gdImageString(im,gdFontGetSmall(),XWIDTH-17-2.75*strlen(s),YWIDTH-13,s,black);
      
 
      /* horizontal tick marks and labels */
      
      mz_range=highest_mass-lowest_mass;
      n_horizontal_labels=horizontal_pixels/50; /* first approximation */
      n_vertical_labels=vertical_pixels/25;
      
      horizontal_label_interval=(int)round((float)mz_range/(float)n_horizontal_labels);
      if(horizontal_label_interval<1) horizontal_label_interval=1;
      if((horizontal_label_interval>2)&&(horizontal_label_interval<8)) horizontal_label_interval=5;
      if((horizontal_label_interval>7)&&(horizontal_label_interval<15)) horizontal_label_interval=10;
      if((horizontal_label_interval>14)&&(horizontal_label_interval<35)) horizontal_label_interval=20;
      if((horizontal_label_interval>34)&&(horizontal_label_interval<75)) horizontal_label_interval=50;
      if((horizontal_label_interval>74)&&(horizontal_label_interval<150)) horizontal_label_interval=100;
      if((horizontal_label_interval>149)&&(horizontal_label_interval<250)) horizontal_label_interval=200;
      if((horizontal_label_interval>249)&&(horizontal_label_interval<350)) horizontal_label_interval=300;
      
      n_horizontal_labels=(int)floor((float)mz_range/(float)horizontal_label_interval);
      
      first_horizontal_label=(int)ceil((float)lowest_mass/(float)horizontal_label_interval)*horizontal_label_interval;
      
      for(i=0;i<n_horizontal_labels;i++)
	{
	  x=(int)round((float)(first_horizontal_label+i*horizontal_label_interval-lowest_mass)/mz_per_pixel)+40;
	  sprintf(s,"%i",first_horizontal_label+i*horizontal_label_interval);
	  if( (x-2.75*strlen(s))<(XWIDTH-50) ) 
	    {
	      gdImageString(im,gdFontGetSmall(),x-2.75*strlen(s),YWIDTH-13,s,black);
	      gdImageLine(im,x,YWIDTH-20,x,YWIDTH-15,black); /* horizontal tick marks */
	    }
	}
      
      
      /* vertical tick marks and labels */
      
      vertical_label_interval=(int)round((float)MS_scans/(float)n_vertical_labels);
      if(vertical_label_interval<1) vertical_label_interval=1;
      if((vertical_label_interval>2)&&(vertical_label_interval<8)) vertical_label_interval=5;
      if((vertical_label_interval>7)&&(vertical_label_interval<15)) vertical_label_interval=10;
      if((vertical_label_interval>14)&&(vertical_label_interval<35)) vertical_label_interval=20;
      if((vertical_label_interval>34)&&(vertical_label_interval<75)) vertical_label_interval=50;
      if((vertical_label_interval>74)&&(vertical_label_interval<150)) vertical_label_interval=100;
      if((vertical_label_interval>149)&&(vertical_label_interval<250)) vertical_label_interval=200;
      if((vertical_label_interval>249)&&(vertical_label_interval<350)) vertical_label_interval=300;
      if((vertical_label_interval>349)&&(vertical_label_interval<450)) vertical_label_interval=400;
      if((vertical_label_interval>449)&&(vertical_label_interval<750)) vertical_label_interval=500;
      if((vertical_label_interval>749)&&(vertical_label_interval<1250)) vertical_label_interval=1000;
      if((vertical_label_interval>1249)&&(vertical_label_interval<1750)) vertical_label_interval=1500;
      if((vertical_label_interval>1449)&&(vertical_label_interval<2500)) vertical_label_interval=2000;
      
      n_vertical_labels=(int)floor((float)MS_scans/(float)vertical_label_interval)-1;
      
      first_vertical_label=(int)ceil((float)first_scan/(float)vertical_label_interval)*vertical_label_interval;
      
      for(i=0;i<n_vertical_labels;i++)
	{
      y=YWIDTH-(int)round((float)(first_vertical_label+i*vertical_label_interval-first_scan)/scans_per_pixel)-20;
      sprintf(s,"%i",scan_number[first_vertical_label+i*vertical_label_interval-first_scan]);
      // sprintf(s,"%i",(int)round((float)(first_vertical_label+i*vertical_label_interval))); /* plot MS scans - future option */
      gdImageString(im,gdFontGetSmall(),30-5.5*strlen(s),y-7,s,black);
      gdImageLine(im,35,y,40,y,black); /* vertical tick marks */
	}
      gdImageString(im,gdFontGetSmall(),9,4,"scan",black);
      // gdImageString(im,gdFontGetSmall(),9,14,"(MS)",black); /* MS scans - future option */
      printf("done ");
    }
  
  gdImagePng(im, outp);
  fclose(outp);
  gdImageDestroy(im);
  printf("(image saved to '%s')\n",image_filename); fflush(stdout);
  
  return 0;
}