//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Mon Apr  9 14:01:16 PDT 2001
// Last Modified: Mon Apr  9 14:01:20 PDT 2001
// Last Modified: Sun May 13 16:15:27 PDT 2001 (added KS algorithm)
// Last Modified: Mon Dec  1 19:07:40 PST 2003 (misc bug fixes)
// Filename:      ...sig/examples/all/keyscape2.cpp
// Web Address:   http://sig.sapp.org/examples/museinfo/humdrum/keyscape2.cpp
// Syntax:        C++; museinfo
//
// Description:   generate a continuous key landscape picture or numeric data.
// 

#include "humdrum.h"
#include <string.h>
#include <stdio.h>
#include <math.h>

typedef Array<float> ArrayFloat;
typedef Array<double> ArrayDouble;

// function declarations
void   checkOptions        (Options& opts, int argc, char* argv[]);
void   example             (void);
void   usage               (const char* command);
void   makeGrayscale       (void);
void   printMatrix         (Array<ArrayInt>& matrix);
void   printTriangle       (Array<ArrayInt>& matrix);
void   printRoots          (Array<ArrayInt>& matrix);
void   printPPM            (Array<ArrayInt>& matrix);
void   printCAPPM          (Array<ArrayInt>& matrix, 
                              Array<ArrayFloat>& claritymatrix, 
                              Array<ArrayFloat>& ambigmatrix);
void   printCPPM           (Array<ArrayInt>& matrix,
                              Array<ArrayFloat>& claritymatrix);
void   printAPPM           (Array<ArrayInt>& matrix,
                              Array<ArrayFloat>& ambigmatrix);
void   printIAPPM          (Array<ArrayInt>& matrix,
                              Array<ArrayFloat>& ambigmatrix,
                              Array<ArrayInt>& secondmatrix);

void   printAveragePPM     (Array<ArrayInt>& matrix);
int    getStartIndex       (HumdrumFile& infile, double startbeat);
int    getStopIndex        (HumdrumFile& infile, double stopbeat);
double difference          (Array<double>& tonicscores, int max);
void   scaleColor          (int& R, int& G, int& B, double cfraction, 
                              double afraction);
void   scaleColorInterpolate      
                           (int& R, int& G, int& B, int sec, double afraction);
void   rgb2hsi             (double *hsi, double *rgb);
void   hsi2rgb             (double* rgb, double* hsi);
int    secondbest          (Array<double>& tonicscores, int max);
int    secondbestmax       (Array<double>& tonicscores, int max);
void   generateLineCoefficients   (HumdrumFile& infile);
void   generatePictureKS   (HumdrumFile& infile);
void   generatePictureCS   (HumdrumFile& infile);

// command line options:
Options   options;              // database for command-line arguments
int       debugQ       = 0;     // used with --debug option
int       compoundQ    = 1;     // used with -c option
double    empirical1   = -4;    // used with --e1 option
double    empirical2   = -3;    // used with --e2 option
double    sx           = 0.578; // used with --sx option
double    sy           = 1.0;   // used with --sx option
int       ppmQ         = 0;     // used with --ppm option
int       imageWidth     = 1000;  // used with -m option
int       imageHeight    = 500;   // used with -t option
int       linearQ      = 0;     // used with -l option
int       monochromeQ  = 0;     // used with -b option
int       matrixQ      = 0;     // used with -x option
int       mirrorQ      = 0;     // used with --mirror option
int       maxdivisions = -1;    // used with -d option
int       averageQ     = 0;     // used with -a option
int       algorithm    = 0;     // used with --algorithm option
int       clarityQ     = 0;     // used with -c option
int       ambigQ       = 0;     // used with -a option
int       secondQ      = 0;     // used with --second option
int       ksalgorithmQ = 0;     // used with -k option
int       rhythmQ      = 1;     // used with -q option
int       rootQ        = 0;     // used with --root option
int       interpolateQ = 0;     // used with -i option
int       gradientQ    = 0;     // used with -g option
int       lineQ        = 0;     // used with --line option
int       linenum      = 0;     // used with --line option
double    sfactor      = 5.0;   // range scaling factor


///////////////////////////////////////////////////////////////////////////

// Pitch to Color translations

int red[40] = {
	0, 9, 18, 0, 63, 63, 63, 73, 82, 0, 218, 237, 255, 255, 255,
	255, 255, 255, 218, 182, 0, 45, 54, 63, 63, 63, 0, 109, 118,
	127, 145, 164, 0, 255, 255, 255, 255, 255, 73, 36
   };

int green[40] = {
	255, 246, 237, 0, 123, 109, 95, 86, 77, 0, 9, 4, 0, 18, 36, 218,
	237, 255, 255, 255, 0, 209, 200, 191, 177, 164, 0, 50, 41, 31,
	27, 22, 0, 91, 109, 127, 145, 164, 255, 255
   };

int blue[40] = {
	0, 36, 73, 0, 255, 255, 255, 255, 255, 0, 73, 36, 0, 0, 0, 0,
	0, 0, 0, 0, 0, 182, 218, 255, 255, 255, 0, 255, 255, 255, 219,
	182, 0, 0, 0, 0, 0, 0, 0, 0
   };


///////////////////////////////////////////////////////////////////////////

double duration = 0.0;
int totalrows = 0;
Array<ArrayInt>   matrix;         // storage for root analysis results
Array<ArrayFloat> claritymatrix;  // storage for root clarity
Array<ArrayFloat> ambigmatrix;    // storage for root ambiguity
Array<ArrayInt>   secondmatrix;   // storage for second best root
Array<double>     parameters(3);


int main(int argc, char* argv[]) {
   HumdrumFile infile;           // input Humdrum Format file

   // process the command-line options
   checkOptions(options, argc, argv);

   // figure out the number of input files to process
   int numinputs = options.getArgCount();

   for (int i=0; i<numinputs || i==0; i++) {
      infile.clear();
      // if no command-line arguments read data file from standard input
      if (numinputs < 1) {
         infile.read(cin);
      } else {
         infile.read(options.getArg(i+1));
      }
      infile.analyzeRhythm("4");   // force the beat to be quarter notes for now

      duration = infile.getTotalDuration();
      if (debugQ) {
         cout << "Duration of input file is " << duration 
              << " quarter notes" << endl;
      }
      totalrows = imageHeight;

      if (lineQ) {
         generateLineCoefficients(infile);
         continue;
      } else if (ksalgorithmQ) {
         generatePictureKS(infile);
      } else {
         generatePictureCS(infile);
      }

      if (rootQ) {
         printRoots(matrix);
      } else if (clarityQ) {
         printCPPM(matrix, claritymatrix);
      } else if (ambigQ) {
         printAPPM(matrix, ambigmatrix);
      } else if (interpolateQ) {
         printIAPPM(matrix, ambigmatrix, secondmatrix);
      } else if (gradientQ) {
         printCAPPM(matrix, claritymatrix, ambigmatrix);
      } else if (secondQ) {
         printPPM(secondmatrix);
      } else {
         printPPM(matrix);
      }
   }

   return 0;
}


///////////////////////////////////////////////////////////////////////////

//////////////////////////////
//
// initializeMatrices -- set all values in the matrices to -1.
//

void initializeMatrices(void) {
   int row, column;

   matrix.setSize(totalrows);
   secondmatrix.setSize(totalrows);
   claritymatrix.setSize(totalrows);
   ambigmatrix.setSize(totalrows);
   for (row=0; row<totalrows; row++) {
      matrix[row].setSize(imageWidth);
      matrix[row].allowGrowth(0);
      secondmatrix[row].setSize(imageWidth);
      secondmatrix[row].allowGrowth(0);
      claritymatrix[row].setSize(imageWidth);
      claritymatrix[row].allowGrowth(0);
      ambigmatrix[row].setSize(imageWidth);
      ambigmatrix[row].allowGrowth(0);
      for (column=0; column<totalrows; column++) {
         matrix[row][column]         = -1;
         secondmatrix[row][column]   = -1;
         claritymatrix[row][column]  = -1;
         ambigmatrix[row][column]    = -1;
      }
   }
}



//////////////////////////////
//
// generatePictureCS --
//

void generatePictureCS(HumdrumFile& infile) {
   int row, column;
   int rrow     = -1;
   int lastrrow = -2;

   initializeMatrices();

   Array<double> tonicscores(40);  // the scores for each possible key
   tonicscores.zero();             // initialize the scores to zero

   double startbeat;               // absolute start beat of analysis window
   double stopbeat;                // absolute ending beat of analysis window

   double centerbeat;              // absolute center of the analysis window
   double windowhalfwidth;         // 1/2 of the analysis window duration

   for (row=0; row<imageHeight; row++) {
      if (linearQ) {
         rrow = imageHeight - row;
      } else {
         rrow = (int)(imageHeight * 
                (1.0-log10((double)(imageHeight-row))/log10(imageHeight+1.0)));
         rrow = imageHeight - rrow + 1;
      }
      if (lastrrow == rrow) {
         matrix[row] = matrix[row-1];
         secondmatrix[row] = secondmatrix[row-1];
         ambigmatrix[row] = ambigmatrix[row-1];
         claritymatrix[row] = claritymatrix[row-1];
         continue;
      }
      for (column=0; column<imageWidth; column++) {
         centerbeat = duration * (1 + 2 * column) / (2.0 * imageWidth);
         windowhalfwidth = duration * (1 - (double)(rrow-1)/imageHeight)/2.0;
               
         startbeat = centerbeat - windowhalfwidth;
         stopbeat  = centerbeat + windowhalfwidth;

         if (fabs(startbeat) < 0.001) {
            startbeat = 0.0;
         }
         if (fabs(stopbeat-duration) < 0.001) {
            startbeat = duration;
         }
         if (startbeat < 0.0 || stopbeat > duration) {
            // if the window falls anywhere outside the music put empty pixel
            continue;
         }

         matrix[row][column] = infile.measureChordRoot(tonicscores,
               parameters, startbeat, stopbeat, algorithm);
         if (matrix[row][column] != -1) {
            claritymatrix[row][column] = tonicscores[matrix[row][column]];
            ambigmatrix[row][column] = 
                  difference(tonicscores, matrix[row][column]);
            secondmatrix[row][column] = 
                  secondbest(tonicscores, matrix[row][column]);
         }
      } 
      lastrrow = rrow;
   }
}



//////////////////////////////
//
// printRoots -- print the roots instead of a picture
//

void printRoots(Array<ArrayInt>& matrix) {
   int i, j;
   char buffer[128];
   for (i=matrix.getSize()-1; i>=0; i--) {
      for (j=0; j<matrix[i].getSize(); j++) {
         if (matrix[i][j] < 0) {
            cout << "x   ";
         } else {
            Convert::base40ToKern(buffer, matrix[i][j] + 2 + 3 * 40);
            strcat(buffer, "     ");
            buffer[4] = '\0';
            cout << buffer;
         }
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// generatePictureKS -- generate Krumhansl-Schmuckler key picture
//

void generatePictureKS(HumdrumFile& infile) {
   initializeMatrices();

   double startbeat;          // absolute start beat of analysis
   double stopbeat;           // absolute ending beat of analysis
   double laststartbeat = -1.0;
   double laststopbeat = -1.0;
   Array<double> scores(24);  // storage for root scores
   double centerbeat;
   double beatwidth;
   int laststartindex;
   int laststopindex;
   int startindex;
   int stopindex;
   int lastroot = -1;
   int jj = -1;
   int lastjj = -2;
   int j, k;
   for (j=0; j<imageHeight; j++) {
      laststartbeat = -1.0;
      laststopbeat = -1.0;
      lastroot = -1;
      if (!linearQ) {
         jj = (int)(imageHeight * (1.0 - 
               log10((double)(imageHeight-j))/log10((double)(imageHeight+1))));
         jj = imageHeight - jj + 1;
      } else {
         jj = imageHeight - j;
      }
      if (jj == lastjj) {
         for (k=0; k<imageWidth; k++) {
            matrix[j][k] = matrix[j-1][k];
            secondmatrix[j][k] = secondmatrix[j-1][k];
            // claritymatrix[j][k] = claritymatrix[j-1][k];
            // ambigmatrix[j][k] = ambigmatrix[j-1][k];
         }
      } else {
         startindex = 0;
         stopindex = 0;
         laststartindex = -1;
         laststopindex = -1;
         for (k=0; k<imageWidth; k++) {
            centerbeat = (k+0.5) * duration / imageWidth;
//            beatwidth = duration * (1.0 - (double)jj/(imageHeight+1));
            beatwidth = duration * (j+1)/imageHeight;
            startbeat = centerbeat - fabs(beatwidth/2.0);
            stopbeat = centerbeat + fabs(beatwidth/2.0) - 0.001;
            if (startbeat < 0.0) {
               stopbeat += -startbeat;
               startbeat = 0.0;
               matrix[j][k] = -1;
               continue;
            }
            if (stopbeat > duration) {
               startbeat -= stopbeat - duration;
               stopbeat = duration;
               matrix[j][k] = -1;
               continue;
            }
            if (startbeat < 0.0) {
               startbeat = 0.0;
            }
            if (startbeat == laststartbeat && stopbeat == laststopbeat) {
               matrix[j][k] = lastroot;
               continue;
            }
            startindex = getStartIndex(infile, startbeat);
            stopindex = getStopIndex(infile, stopbeat);
            if (debugQ) {
               cout << "Window beats: " << startbeat 
                    << "(" << startindex << ")"
                    << "\tto " << stopbeat 
                    << "(" << stopindex << ")"
                    << "\trow = " << j 
                    << "\tcol = " << k 
                    << endl;
            }

            matrix[j][k] = infile.analyzeKeyKS(scores, startindex,
                  stopindex, rhythmQ);
            matrix[j][k] = (Convert::base12ToBase40(matrix[j][k])-2+40)% 40;
            if (matrix[j][k] != -1) {
               // claritymatrix[j][k] = tonicscores[matrix[j][k]];
               // ambigmatrix[j][k] = 
               //       difference(tonicscores, matrix[j][k]);
               secondmatrix[j][k] = secondbestmax(scores, matrix[j][k]);
            }
            laststartindex = startindex;
            laststopindex = stopindex;
         } 
      }
      lastjj = jj;
   }
}



//////////////////////////////
//
// generateLineCoefficients -- list 40 sets of coefficients for a
//    particular line in a plot
//

void generateLineCoefficients(HumdrumFile& infile) {
      int j, k;
      Array<ArrayDouble> matrix(imageWidth);
      for (k=0; k<imageWidth; k++) {
         matrix[k].setSize(40);
         matrix[k].allowGrowth(0);
         for (j=0; j<40; j++) {
            matrix[k][j] = -1;
         }
      }

      int analysisstart = -1;
      int analysisstop = -1;

      Array<double> parameters(3);
      parameters[0] = sx;
      parameters[1] = empirical1;
      parameters[2] = empirical2;
  
      double startbeat;     // absolute start beat of analysis
      double stopbeat;      // absolute ending beat of analysis
      double laststartbeat = -1.0;
      double laststopbeat = -1.0;
      double centerbeat;
      double beatwidth;
      int laststartindex;
      int laststopindex;
      int startindex;
      int stopindex;
      int lastroot = -1;
      int jj = -1;
      for (j=linenum; j<=linenum; j++) {
         laststartbeat = -1.0;
         laststopbeat = -1.0;
         lastroot = -1;
         if (!linearQ) {
            jj = (int)(imageHeight * (1.0 - 
                log10((double)(imageHeight-j))/log10((double)(imageHeight+1))));
            jj = imageHeight - jj + 1;
         } else {
            jj = imageHeight - j;
         }

         startindex = 0;
         stopindex = 0;
         laststartindex = -1;
         laststopindex = -1;
         for (k=0; k<imageWidth; k++) {
            centerbeat = k * duration / imageWidth;
            beatwidth = duration * (1.0 - (double)jj/(imageHeight+1));
            startbeat = centerbeat - fabs(beatwidth/2.0);
            stopbeat = centerbeat + fabs(beatwidth/2.0);
            if (startbeat < 0.0) {
               stopbeat += -startbeat;
               startbeat = 0.0;
               matrix[k].zero();
               continue;
            }
            if (stopbeat > duration) {
               startbeat -= stopbeat - duration;
               stopbeat = duration;
               if (analysisstop == -1) {
                  analysisstop = k - 1;
               }
               matrix[k].zero();
               continue;
            }
            if (analysisstart == -1) {
               analysisstart = k;
            }
            if (startbeat < 0.0) {
               startbeat = 0.0;
            }
            startindex = getStartIndex(infile, startbeat);
            stopindex = getStopIndex(infile, stopbeat);
            int flag = infile.measureChordRoot(matrix[k],
                          parameters, startindex, stopindex, algorithm);
            if (flag == -1) {
               matrix[k].zero();
            }
         laststartindex = startindex;
         laststopindex = stopindex;
      } 
   }

   int m, n;
   for (m=0; m<40; m++) {
      cout << "root" << m << " = Graphics[{RGBColor["
           << red[m]/256.0 << ", "
           << green[m]/256.0 << ", "
           << blue[m]/256.0 << "], "
           << "Line[{";
      for (n=analysisstart; n<=analysisstop; n++) {
         cout << " {" << n;
         cout << ", " << 1.0/matrix[n][m] << "}";
         if (n < analysisstop) {
            cout << ",";
         }
      }
      cout << "}]}];\n";
   }
}



//////////////////////////////
//
// secondbest -- returns the second best root 
//

int secondbest(Array<double>& tonicscores, int max) {
   int max2 = 0;
   if (max == 0) max2 = 1;
   for (int i=1; i<tonicscores.getSize(); i++) {
      if (i == max) {
         continue;
      }
      // actually looking for a minimum right now
      if (tonicscores[i] < tonicscores[max2]) {
         max2 = i;
      }
   }

   return max2;
}



//////////////////////////////
//
// secondbestmax -- returns the second best root 
//

int secondbestmax(Array<double>& tonicscores, int max) {
   int max2 = 0;
   if (max == 0) max2 = 1;
   for (int i=1; i<tonicscores.getSize(); i++) {
      if (i == max) {
         continue;
      }
      // actually looking for a minimum right now
      if (tonicscores[i] > tonicscores[max2]) {
         max2 = i;
      }
   }

   return max2;
}



//////////////////////////////
//
// difference -- returns the difference between the highest score 
//    and the second highest score.
//

double difference(Array<double>& tonicscores, int max) {
   int max2 = 0;
   if (max == 0) max2 = 1;
   for (int i=1; i<40; i++) {
      if (i == max) {
         continue;
      }
      // actually looking for a minimum right now
      if (tonicscores[i] < tonicscores[max2]) {
         max2 = i;
      }
   }

   return tonicscores[max2] - tonicscores[max];
}



//////////////////////////////
//
// printCPPM --  print only the clarity 
//

void printCPPM(Array<ArrayInt>& matrix, 
      Array<ArrayFloat>& claritymatrix) {
   int row, column;
   int R, G, B;
   int maxrow = imageHeight;
   int maxcolumn = imageWidth;
   int i, j;
   int root;

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "256\n";

   double minc, maxc;
   minc = maxc = 0.0;
   int minit;

   for (row=maxrow-1; row>=0; row--) {
      i = row;
      // find the max and min in claritymatrix
      minit = 0;
      for (column=0; column<maxcolumn; column++) {
         j = column;
         if (matrix[i][j] < 0) {
            continue;
         }
         if (minit == 0) {
            maxc = minc = claritymatrix[i][j];
            minit = 1;
            continue;
         }
         if (claritymatrix[i][j] < minc)  minc = claritymatrix[i][j];
         if (claritymatrix[i][j] > maxc)  maxc = claritymatrix[i][j];
      }
      if (minc == maxc) {
         minc -= 0.5;
         maxc += 0.5;
      }

      double cfraction = 0;
      for (column=0; column<maxcolumn; column++) {
         j = column;
         root = matrix[i][j];
         if (root >= 0) {
            cfraction = 1.0 - (claritymatrix[i][j] - minc)/(maxc-minc);
            cfraction = pow(cfraction, 1/sfactor);
            R = (int)(cfraction * 256 + 0.5);
            G = (int)(cfraction * 256 + 0.5);
            B = (int)(cfraction * 256 + 0.5);
            cout << R << " " << G << " " << B << "  ";
         } else {
            cout << " 255 0 0  ";
         }
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// printAPPM --  print only the ambiguity 
//

void printAPPM(Array<ArrayInt>& matrix, 
      Array<ArrayFloat>& ambigmatrix) {
   int row, column;
   int R, G, B;
   int maxrow = imageHeight;
   int maxcolumn = imageWidth;
   int i, j;
   int root;

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "256\n";

   double minc, maxc;
   minc = maxc = 0.0;
   int minit;

   for (row=maxrow-1; row>=0; row--) {
      i = row;
      // find the max and min in ambigmatrix
      minit = 0;
      for (column=0; column<maxcolumn; column++) {
         j = column;
         if (matrix[i][j] < 0) {
            continue;
         }
         if (minit == 0) {
            maxc = minc = ambigmatrix[i][j];
            minit = 1;
            continue;
         }
         if (ambigmatrix[i][j] < minc)  minc = ambigmatrix[i][j];
         if (ambigmatrix[i][j] > maxc)  maxc = ambigmatrix[i][j];
      }
      if (minc == maxc) {
         minc -= 0.5;
         maxc += 0.5;
      }

      double cfraction = 0;
      for (column=0; column<maxcolumn; column++) {
         j = column;
         root = matrix[i][j];
         if (root >= 0) {
            cfraction = (ambigmatrix[i][j] - minc)/(maxc-minc);
            cfraction = pow(cfraction, 1/sfactor);
            R = (int)(cfraction * 256 + 0.5);
            G = (int)(cfraction * 256 + 0.5);
            B = (int)(cfraction * 256 + 0.5);
            cout << R << " " << G << " " << B << "  ";
         } else {
            cout << " 0 0 255  ";
         }
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// printIAPPM -- 
//

void printIAPPM(Array<ArrayInt>& matrix, 
      Array<ArrayFloat>& ambigmatrix, 
      Array<ArrayInt>& secondmatrix) {
   int row, column;
   int R, G, B;
   int maxrow = imageHeight;
   int maxcolumn = imageWidth;
   int i, j;
   int root;

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "256\n";

   double mina, maxa;
   mina = maxa = 0.0;
   int minit;

   for (row=maxrow-1; row>=0; row--) {
      i = row;
      // find the max and min in ambigmatrix
      minit = 0;
      for (column=0; column<maxcolumn; column++) {
         j = column;
         if (matrix[i][j] < 0) {
            continue;
         }
         if (minit == 0) {
            maxa = mina = ambigmatrix[i][j];
            minit = 1;
            continue;
         }
         if (ambigmatrix[i][j] < mina)    mina = ambigmatrix[i][j];
         if (ambigmatrix[i][j] > maxa)    maxa = ambigmatrix[i][j];
      }
      if (mina == maxa) {
         mina -= 0.5;
         maxa += 0.5;
      }

      for (column=0; column<maxcolumn; column++) {
         j = column;
         root = matrix[i][j];
         if (root >= 0) {
            R = red[root]; 
            G = green[root];
            B = blue[root];
            scaleColorInterpolate(R, G, B, secondmatrix[i][j],
                  (ambigmatrix[i][j] - mina)/(maxa-mina));
            cout << R << " " << G << " " << B << "  ";
         } else {
            cout << " 0 0 0  ";
         }
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// printCAPPM -- 
//

void printCAPPM(Array<ArrayInt>& matrix, 
      Array<ArrayFloat>& claritymatrix, 
      Array<ArrayFloat>& ambigmatrix) {
   int row, column;
   int R, G, B;
   int maxrow = imageHeight;
   int maxcolumn = imageWidth;
   int i, j;
   int root;

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "256\n";

   double minc, maxc;
   double mina, maxa;
   minc = maxc = mina = maxa = 0.0;
   int minit;

   for (row=maxrow-1; row>=0; row--) {
      i = row;
      // find the max and min in claritymatrix
      minit = 0;
      for (column=0; column<maxcolumn; column++) {
         j = column;
         if (matrix[i][j] < 0) {
            continue;
         }
         if (minit == 0) {
            maxc = minc = claritymatrix[i][j];
            maxa = mina = ambigmatrix[i][j];
            minit = 1;
            continue;
         }
         if (claritymatrix[i][j] < minc)  minc = claritymatrix[i][j];
         if (claritymatrix[i][j] > maxc)  maxc = claritymatrix[i][j];
         if (ambigmatrix[i][j] < mina)    mina = ambigmatrix[i][j];
         if (ambigmatrix[i][j] > maxa)    maxa = ambigmatrix[i][j];
      }
      if (minc == maxc) {
         minc -= 0.5;
         maxc += 0.5;
      }
      if (mina == maxa) {
         mina -= 0.5;
         maxa += 0.5;
      }

      for (column=0; column<maxcolumn; column++) {
         j = column;
         root = matrix[i][j];
         if (root >= 0) {
            R = red[root]; 
            G = green[root];
            B = blue[root];
            scaleColor(R, G, B, (claritymatrix[i][j] - minc)/(maxc-minc),
                  (ambigmatrix[i][j] - mina)/(maxa-mina));
            cout << R << " " << G << " " << B << "  ";
         } else {
            cout << " 0 0 0  ";
         }
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// scaleColorInterpolate -- interpolate between second best and best
//      root colors
//

void scaleColorInterpolate(int& R, int& G, int& B, int sec, double afraction) {
   double rgbA[3];
   rgbA[0] = R/256.0;
   rgbA[1] = G/256.0;
   rgbA[2] = B/256.0;

   double hsiA[3] = {0};
   rgb2hsi(hsiA, rgbA);

   double rgbB[3];
   rgbB[0] = red[sec]/256.0;
   rgbB[1] = green[sec]/256.0;
   rgbB[2] = blue[sec]/256.0;

   double hsiB[3] = {0};
   rgb2hsi(hsiB, rgbB);

   double hsiC[3] = {0};
   double rgbC[3] = {0};

   // interpolate the hsi values by the ambiguity fraction
   double difference = fabs(hsiA[0] - hsiB[0]);
   double pointA = hsiA[0];
   double pointM;
   if (difference < 0.5) {
      pointM = (pointA + hsiB[0])/2;
   } else if (pointA < 0.5) {
      pointM = (pointA + hsiB[0] - 1.0)/2;
   } else {
      pointM = (pointA + hsiB[0] + 1.0)/2;
   }

   if (pointA < pointM) {
      hsiC[0] = (pointM - pointA) * (1.0 - afraction) + pointA;
   } else {
      hsiC[0] = (pointA - pointM) * afraction + pointM;
   }
   // hsiC[0] = (pointA + pointB)/2;
   // hsiC[0] = pointM;
   if (hsiC[0] > 1.0) hsiC[0] += -1.0;
   if (hsiC[0] < 0.0) hsiC[0] += +1.0;

   // hsiC[1] = hsiA[1] * (1.0 - afraction) + hsiB[1] * afraction;
   //if (afraction < 0.5) {
   //   hsiC[2] = (hsiA[2] + hsiB[2])/2;
   //} else {
   //   hsiC[2] = hsiA[2];
   // }

   hsiC[1] = (hsiA[1] + hsiB[1])/2;
   hsiC[2] = (hsiA[2] + hsiB[2])/2;

   if (hsiC[0] < 0) {
      hsiC[0] += 1.0;
   }

   if (hsiC[0] < 0.0) {
      hsiC[0] += 1.0;
   }
   if (hsiC[0] > 1.0) {
      hsiC[0] -= 1.0;
   }

   if (hsiC[1] < 0.0) {
      hsiC[1] += 1.0;
   }
   if (hsiC[1] > 1.0) {
      hsiC[1] -= 1.0;
   }

   if (hsiC[2] < 0.0) {
      hsiC[2] += 1.0;
   }
   if (hsiC[2] > 1.0) {
      hsiC[2] -= 1.0;
   }

   hsi2rgb(rgbC, hsiC);

   if (rgbA[0] < rgbB[0]) {
      rgbC[0] = (rgbB[0] - rgbA[0])/2 * pow((1.0 - afraction), sfactor) + rgbA[0];
   } else {
      rgbC[0] = rgbA[0] - (rgbA[0] - rgbB[0])/2 * pow((1.0 - afraction), sfactor);
   }

   if (rgbA[1] < rgbB[1]) {
      rgbC[1] = (rgbB[1] - rgbA[1])/2 * pow((1.0 - afraction), sfactor) + rgbA[1];
   } else {
      rgbC[1] = rgbA[1] - (rgbA[1] - rgbB[1])/2 * pow((1.0 - afraction), sfactor);
   }

   if (rgbA[2] < rgbB[2]) {
      rgbC[2] = (rgbB[2] - rgbA[2])/2 * pow((1.0 - afraction), sfactor) + rgbA[2];
   } else {
      rgbC[2] = rgbA[2] - (rgbA[2] - rgbB[2])/2 * pow((1.0 - afraction), sfactor);
   }

   // rgbC[1] = (rgbA[1] - rgbB[1]) * afraction + rgbA[0];
   // rgbC[2] = (rgbA[2] - rgbB[2]) * afraction + rgbA[0];

   R = (int)(rgbC[0] * 256.0 + 0.5);
   G = (int)(rgbC[1] * 256.0 + 0.5);
   B = (int)(rgbC[2] * 256.0 + 0.5);
}



//////////////////////////////
//
// scaleColor -- add clarity and ambiguity color shifting to root
//      analysis
//

void scaleColor(int& R, int& G, int& B, double cfraction, double afraction) {
   double rgb[3];
   rgb[0] = R/256.0;
   rgb[1] = G/256.0;
   rgb[2] = B/256.0;
   double hsi[3] = {0};
   rgb2hsi(hsi, rgb);

   // scale the saturation by the clarity fraction
   hsi[1] = hsi[1] * pow(1.0 - cfraction, 1/2.0);
   if (hsi[1] < 0.0) {
      hsi[1] = 0.0;
   }
   if (hsi[1] > 1.0) {
      hsi[1] = 1.0;
   }

   // scale the intensity by the ambiguous fraction
   hsi[2] = hsi[2] * pow(afraction, 1/sfactor);
   if (hsi[2] < 0.0) {
      hsi[2] = 0.0;
   }
   if (hsi[2] > 1.0) {
      hsi[2] = 1.0;
   }

   hsi2rgb(rgb, hsi);
   R = (int)(rgb[0] * 256.0 + 0.5);
   G = (int)(rgb[1] * 256.0 + 0.5);
   B = (int)(rgb[2] * 256.0 + 0.5);
}



//////////////////////////////
//
// printPPM -- 
//

void printPPM(Array<ArrayInt>& matrix) {
   int row, column;
   int maxrow = imageHeight;
   int maxcolumn = imageWidth;
   int i, j;
   int root;

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "256\n";

   for (row=maxrow-1; row>=0; row--) {
      i = row;
      // i = (int)(maxrow * (1.0 - log10((double)(maxrow-row)) /
      //         log10((double)(maxrow+1))));
      // i = maxrow - i + 1;
      for (column=0; column<maxcolumn; column++) {
         j = column;
         root = matrix[i][j];
         if (root >= 0) {
            cout << red[root]   << " "
                 << green[root] << " "
                 << blue[root]  << "  ";
         } else {
            if (ksalgorithmQ) {
               cout << " 255 255 255 ";
            } else {
               cout << " 0 0 0 ";
            }
         }
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// checkOptions -- validate and process command-line options.
//

void checkOptions(Options& opts, int argc, char* argv[]) {
   opts.define("C|compound=b",  "don't try to use compound meters");   
   opts.define("r|rhythm=i:4",  "rhythm to measure root at");   
   opts.define("e1|empirical1=d:-4",  "empirical rhythm value 1");
   opts.define("e2|empirical2=d:-3",  "empirical rhythm value 2");
   opts.define("sx=d:0.578", "scale factor for the x-pitch space axis");
   opts.define("sy=d:1.0", "scale factor for the y-pitch space axis");
   opts.define("algorithm=i:0",  "0 = default algorithm, 1+ = new algorithms");
   opts.define("roots|root=b", "display root note analysis");

   opts.define("ppm=b", "generate a ppm graphic lanscape image");
   opts.define("d|dimension=s:10x10","pixel width and height of picture");

   opts.define("k|ks=b", "generate Krumhansl-Schmuckler key picture");
   opts.define("c|clarity=b", "display clarity gradations");
   opts.define("a|ambig=b", "display ambiguity gradations");
   opts.define("second=b", "display second best roots");
   opts.define("i|interpolate=b", "ambiguity gradations via interploation");
   opts.define("g|gradient=b", "display gradations of tonal clarity");
   opts.define("line=i:0", "print out root weightings for a horizontal line");

   opts.define("q|quarter=b", "Calculate KS algorithm without durations");
   opts.define("s|sfactor=d:5.0", "gradation scaling factor");
   opts.define("l|linear=b", "linear plot instead of logarithmic");
   opts.define("b|monochrome=b", "use black and white instead of color");
   opts.define("x|matrix=b", "generate a matrix rather than a table");
   opts.define("mirror=b","generate mirror matrix rather than cyclical matrix");
   opts.define("divisions=i:-1", "the maximum number of divisions");
   opts.define("average=b","generate average key picture");

   opts.define("debug=b",  "trace input parsing");   
   opts.define("author=b",  "author of the program");   
   opts.define("version=b", "compilation information"); 
   opts.define("example=b", "example usage"); 
   opts.define("h|help=b",  "short description"); 
   opts.process(argc, argv);
   
   // handle basic options:
   if (opts.getBoolean("author")) {
      cout << "Written by Craig Stuart Sapp, "
           << "craig@ccrma.stanford.edu, Dec 2000" << endl;
      exit(0);
   } else if (opts.getBoolean("version")) {
      cout << argv[0] << ", version: 6 Dec 2000" << endl;
      cout << "compiled: " << __DATE__ << endl;
      cout << MUSEINFO_VERSION << endl;
      exit(0);
   } else if (opts.getBoolean("help")) {
      usage(opts.getCommand());
      exit(0);
   } else if (opts.getBoolean("example")) {
      example();
      exit(0);
   }

   algorithm = opts.getInteger("algorithm");

   debugQ = opts.getBoolean("debug");

   if (opts.getBoolean("compound")) {
      compoundQ = 0;
   } else {
      compoundQ = 0;
   }

   empirical1 = opts.getDouble("empirical1");
   empirical2 = opts.getDouble("empirical2");
   sx = opts.getDouble("sx");
   sy = opts.getDouble("sy");

   ppmQ         = opts.getBoolean("ppm");

   sscanf(opts.getString("dimension"), "%dx%d", &imageWidth, &imageHeight);

   linearQ      = opts.getBoolean("linear");
   matrixQ      = opts.getBoolean("matrix");
   mirrorQ      = opts.getBoolean("mirror");
   maxdivisions = opts.getInteger("divisions");
   averageQ     = opts.getBoolean("average");
   clarityQ     = opts.getBoolean("clarity");
   ambigQ       = opts.getBoolean("ambig");
   secondQ      = opts.getBoolean("second");
   interpolateQ = opts.getBoolean("interpolate");
   gradientQ    = opts.getBoolean("gradient");
   lineQ        = opts.getBoolean("line");
   linenum      = opts.getInteger("line");
   sfactor      = opts.getDouble("sfactor");
   ksalgorithmQ = opts.getBoolean("ks");
   rhythmQ      = !opts.getBoolean("quarter");
   rootQ        = opts.getBoolean("roots");

   if (opts.getBoolean("monochrome")) {
      makeGrayscale();
   }

   parameters[0] = sx;              // alpha
   parameters[1] = empirical1;      // delta
   parameters[2] = empirical2;      // lambda
}



//////////////////////////////
//
// example -- example usage of the maxent program
//

void example(void) {
   cout <<
   "                                                                        \n"
   << endl;
}



//////////////////////////////
//
// usage -- gives the usage statement for the quality program
//

void usage(const char* command) {
   cout <<
   "                                                                        \n"
   << endl;
}



//////////////////////////////
//
// makeGrayscale --
//

void makeGrayscale(void) {
   for (int i=0; i<40; i++) {
      red[i]   = (int)((i+1.0)/40.0);
      green[i] = red[i];
      blue[i]  = red[i];
   }
}


//////////////////////////////
//
// getStartIndex --
//

int getStartIndex(HumdrumFile& infile, double startbeat) {
   int index = 1;
   while (index < infile.getNumLines() - 1) {
      if (startbeat <= infile[index].getAbsBeat() + 0.001 && 
          startbeat > infile[index-1].getAbsBeat() - 0.001) {
         return index;
      }
      index++;
   }

   return infile.getNumLines() - 1;
}



//////////////////////////////
//
// getStopIndex --
//

int getStopIndex(HumdrumFile& infile, double stopbeat) {
   int index = 1;
   while (index < infile.getNumLines()) {
      if (stopbeat <= infile[index].getAbsBeat() + 0.001 && 
          stopbeat > infile[index-1].getAbsBeat() - 0.001) {
         return index;
      }
      index++;
   }

   return infile.getNumLines() - 1;
}



//////////////////////////////
//
// rgb2hsi -- convert from RGB color space to HSI color space.
//

#ifndef PI
   #define PI M_PI
#endif

void rgb2hsi(double *hsi, double *rgb) {
   double& R = rgb[0];
   double& G = rgb[1];
   double& B = rgb[2];
   double& H = hsi[0];
   double& S = hsi[1];
   double& I = hsi[2];

   double min = R;
   if (G < min) min = G;
   if (B < min) min = B;

   I = (R+G+B)/3.0;
   S = 1 - min/I;
   if (S == 0.0) {
      H = 0.0;
   } else {
      H = ((R-G)+(R-B))/2.0;
      H = H/sqrt((R-G)*(R-G) + (R-B)*(G-B));
      H = acos(H);
      if (B > G) {
         H = 2*PI - H;
      }
      H = H/(2*PI);
   }
}



//////////////////////////////
//
// hsi2rgb -- convert from HSI color space to RGB color space.
//

#ifndef PI
   #define PI M_PI
#endif

void hsi2rgb(double* rgb, double* hsi) {
   double& R = rgb[0];
   double& G = rgb[1];
   double& B = rgb[2];
   double& H = hsi[0];
   double& S = hsi[1];
   double& I = hsi[2];

   double r;
   double g;
   double b;

   if (H < 1.0/3.0) {
      b = (1-S)/3;
      r = (1+S*cos(2*PI*H)/cos(PI/3-2*PI*H))/3.0;
      g = 1 - (b + r);
   } else if (H < 2.0/3.0) {
      H = H - 1.0/3.0;
      r = (1-S)/3;
      g = (1+S*cos(2*PI*H)/cos(PI/3-2*PI*H))/3.0;
      b = 1 - (r+g);
   } else {
      H = H - 2.0/3.0;
      g = (1-S)/3;
      b = (1+S*cos(2*PI*H)/cos(PI/3-2*PI*H))/3.0;
      r = 1 - (g+b);
   }

   R = I * r * 3;
   G = I * g * 3;
   B = I * b * 3;
}