#include <stdio.h>

#include <stdlib.h>
#include <math.h>
#include <tiffio.h>

#define TABLE_SIZE 1024
#define PI 3.1415926535
#define SAMPLE_RATE 44100

/* transition time between volume levels */
#define RAMP_TIME 441

#define SAMPLE uint16

#define N_CHANNELS 2

/* number of different waveforms */
#define N_WAVE_TYPES 3

#define WAVE_SINE 0
#define WAVE_SQUARE 1
#define WAVE_SAW 2

float wav_tab[N_WAVE_TYPES][TABLE_SIZE];
/* below for temporary */
float ***wav_tab_scaled;

/**
 * Converts hz to wavelength given the
 * sampling rate
 */
#define HZ2WLFORSL(hz,sr) ((int)((float)sr / (float)hz))
#define HZ2WL(hz) ((int)(44100.0 / (float)hz))

#define BLOCK_SIZE 32768
SAMPLE buffer[BLOCK_SIZE * N_CHANNELS];
int buf_ix = 0;

float *float_buf;

#define ENERGY_THRESHOLD 0.00005

void make_tables(int wave_type) {
     int i, wt;
     float x, step = (float)(2 * PI) / (float) TABLE_SIZE;
     int half = TABLE_SIZE / 2;
     for(i = 0, x = 0; i < TABLE_SIZE; i++, x += step) {
	  wav_tab[WAVE_SINE]  [i] = sin(x);
	  wav_tab[WAVE_SQUARE][i] = (i < half) ? -1 : 1;
	  wav_tab[WAVE_SAW]   [i] =
	       (i < half) ?
	       (((float)i / (float)half) * 2.0) - 1.0 :
	       (((float)(TABLE_SIZE - i) / (float)half) * 2.0) - 1.0;
     }
}

/**
 * n = number of frequency bands
 * pitch = array with wavelength for each band (in samples)
 * energy = 2d array with energy for each band (1.0 = max, 0.0 = min) and for
 * each color r,g,b
 * energy_prev = 2d what the energy in the band was last frame
 * offset = how many samples into the frame to wait before the transition to
 *  the new energy
 * pan = panning of each freq band
 * ss = starting sample number (to align phase)
 * ns = number of samples to generate
 */
void do_timeSegment(int n, int *pitch, float **energy,
		    float **energy_prev, int *offset, float *pan,
		    long ss, int ns, SAMPLE *buf) {
     long cs; /* current sample */
     long es = ss + ns; /* end sample */
     int i;
     int s; /* sample # within result segment */
     if(!float_buf)
	  float_buf = calloc(ns * N_CHANNELS, sizeof(float));
     for(s = 0; s < ns * N_CHANNELS; s++) {
	  float_buf[s] = 0.0;
     }
     for(i = 0; i < n; i++) { /* for each freq band */
	  int rgb;
	  int wl = pitch[i];
	  float scale = (float)TABLE_SIZE / (float)wl;
	  int start_trans; /* beginning of power transition */
	  int end_trans; /* end of power transition */
	  start_trans = offset[i] - RAMP_TIME;	
	  if(start_trans < 0) start_trans = 0;
	  end_trans = start_trans + RAMP_TIME;
	  if(end_trans >= ns) end_trans = ns;
	  for(rgb = 0; rgb < 3; rgb++) {
	       if(energy_prev[rgb][i] > ENERGY_THRESHOLD ||
		  energy[rgb][i] > ENERGY_THRESHOLD) {
		    float e = energy_prev[rgb][i];
		    for(cs = ss, s = 0; cs < es; cs++, s++) {
			 float amp;
			 if(s > start_trans && s < end_trans) {
			      e = energy_prev[rgb][i] +
				   (((float)(s - start_trans) / (float)RAMP_TIME) *
				    (energy[rgb][i] - energy_prev[rgb][i]));
			 } else if(s >= end_trans) {
			      e = energy[rgb][i];
			 }
			 /* if(s > offset[i]) e = energy[rgb][i]; */
			 /* compute amplitude */
			 amp = (wav_tab_scaled[rgb][i][cs % wl] * e);
			 /* now distribute amp across channels */
			 float_buf[s*2]   += amp * pan[i];
			 float_buf[s*2+1] += amp * (1.0 - pan[i]);
		    }
	       }
	  }
     }
     /* now convert results to unsigned shorts */
     for(s = 0; s < ns * N_CHANNELS; s++) {
	  buf[s] = (SAMPLE)(float_buf[s] * 32767);
     }
}

void display(SAMPLE *buf, int n, int step) {
     int i;
     for(i = 0; i < n; i += step) {
	  int j;
	  int d = buf[i] / (65536.0 / 50.0);
	  for(j = 0; j < d; j++) {
	       putchar(' ');
	  }
	  printf("*\n");
     }
}

void output(SAMPLE *buf, int n) {
     int i;
     for(i = 0; i < n * N_CHANNELS; i++) {
	  buffer[buf_ix++] = buf[i];
	  if(buf_ix==BLOCK_SIZE) {
	       fwrite(buffer, BLOCK_SIZE, sizeof(SAMPLE), stdout);
	       buf_ix = 0;
	  }
     }
}

void final_output() {
     fwrite(buffer, buf_ix, sizeof(SAMPLE), stdout);
}

void show_progress(int seg, int nseg) {
     int i, n;
     float percent = ((float)seg / (float)nseg) * 100;
     fprintf(stderr,"\rcompleted: [");
     n = percent / 3;
     for(i = 0; i < n; i++) {
	  fputc('*',stderr);
     }
     for(; i < 33; i++) {
	  fputc(' ',stderr);
     }
     fprintf(stderr,"] %.2f%%",percent);
     fflush(stderr);
}

int main(int argc, char **argv) {
     int i, j;

     /** output buffering */
     SAMPLE *result;
     long cs;

     /* how long each column lasts */
     int block;

     /** usage: bosch duration tif_file > output */

     /** control params */
     float minFreq = 40.0; /* freq of lowest band */
     float maxFreq = 3500.0; /* freq of highest band */
     float duration = atof(argv[1]); /* how long the image lasts in seconds */
     float freqExp = 2.5; /* log scaling of freq band spacing */
     float rolloff = 2.0; /* log scaling of amplitude of freq bands */
     float minEnergy = 0.1; /* the amplitude scale of the highest pitch band */
     float bandwidth; /* width of each band (assume linear) */
     float band;
     float ampScale;

     /** tiff file params */
     TIFF *tif;
     char *tif_fn = argv[2];
     uint32 w, h;
     uint32 *image;
     int c;

     /* pitches of the freq bands */
     int *pitch;
     /* offset of amp transition for each band */
     int *offset;
     /* energy of the freq bands */
     float *energy[3];
     float *energy_prev[3];
     /* panning of each freq band */
     float *pan;
     /* amp scale of each band (higher bands are quieter) */
     float *freqResponse;
     float x, y; /* temp vars */
     int rgb, wt;

     make_tables(WAVE_SAW);

     /* first arg is TIFF file. read TIFF file into raster */
     tif = TIFFOpen(tif_fn,"r");
     TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &w);
     TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &h);

     image = (uint32*) _TIFFmalloc(w * h * sizeof (uint32));
     if (!image) {
	  fprintf(stderr,"error opening tiff file!");
	  exit(-1);
     }

     if(!TIFFReadRGBAImage(tif, w, h, image, 0)) {
	  fprintf(stderr,"error reading tiff image!");
	  exit(-1);
     }

     fprintf(stderr,"opened tiff image width=%d height=%d\n", w, h);

     /* we know how wide the image is. map that to time */
     block = (duration / (float)w) * SAMPLE_RATE;
     fprintf(stderr,"blocksize: %d\n", block);
     result = calloc(block * N_CHANNELS, sizeof(SAMPLE));

     /* ok, now we know how high the image is. map that to frequencies */
     pitch = calloc(h, sizeof(int));
     /* log scale so low bands are narrower */
     for(i = 0; i < h; i++) {
	  band = minFreq +
	       (pow((float)i / (float)h, freqExp) *
		(maxFreq - minFreq));
	  pitch[i] = HZ2WL(band);
     }

     for(c = 0; c < 3; c++) {
	  energy[c] = calloc(h * N_WAVE_TYPES, sizeof(float));
	  energy_prev[c] = calloc(h * N_WAVE_TYPES, sizeof(float));
     }

     cs = 0;

     /* set random offsets to avoid "buzzing" */
     /* also set panning to random */
     /* also set up frequency response */
     offset = calloc(h, sizeof(int));
     pan = calloc(h, sizeof(float));
     freqResponse = calloc(h, sizeof(float));
     for(i = 0; i < h; i++) {
	  offset[i] = random() % block;
	  pan[i] = ((float) (random() % 1000)) / 1000.0;
	  freqResponse[i] =
	       (pow(1.0 - ((float)i / (float)h), rolloff) *
		(1.0 - minEnergy)) +
	       minEnergy;
     }

     /* now precompute all the wavetables */
     wav_tab_scaled = (float ***)calloc(N_WAVE_TYPES,sizeof(float **));
     for(rgb = 0; rgb < 3; rgb++) {
	  wav_tab_scaled[rgb] = (float **)calloc(h, sizeof(float *));
	  for(i = 0; i < h; i++) { /* for each freq band */
	       int wl = pitch[i];
	       float scale = (float)TABLE_SIZE / (float)wl;
	       int s;
	       wav_tab_scaled[rgb][i] = (float *)calloc(wl,sizeof(float));
	       for(s = 0; s < wl; s++) {
		    int samp = (int)((float)s * scale);
		    wav_tab_scaled[rgb][i][s] = wav_tab[rgb][samp];
	       }
	  }
     }

     /*
       ampScale = 0.0075 / (float)h;
     */
     ampScale = 0.05 / (float)h;
     /* now produce the sound */
     for(i = 0; i < w-1; i++) { /* for each time segment, */
	  for(j = 0; j < h; j++) { /* for each frequency band */
	       /* compute the energy based on color */
	       int c;
	       uint8 rgb[3];
	       uint32 pixel;

	       pixel = image[(j * w) + i];

	       rgb[0] = (uint8)((pixel & 0x00FF0000) >> 16);
	       rgb[1] = (uint8)((pixel & 0x0000FF00) >> 8);
	       rgb[2] = (uint8)((pixel & 0x000000FF) >> 0);

	       for(c = 0; c < 3; c++) {
		    energy_prev[c][j] = energy[c][j];
		    /* energy for each wave type is inverse color intensity */
		    energy[c][j] = ampScale * freqResponse[j] *
			 (255 - rgb[c]);
	       }
	  }
	  do_timeSegment(h, pitch, energy, energy_prev, offset, pan, cs, block, result);
	  output(result, block);
	  cs += block;
	  if(!(i % 10)) {
	       show_progress(i,w-1);
	  }
     }
     final_output();
     show_progress(100,100);
     fprintf(stderr,"\n");

/*
     result = do_timeSegment(3, pitch, energy, 0, block);
     output(result,block);
     result = do_timeSegment(3, pitch, energy, block, block);
     output(result,block);
*/

     return 1;
}