// Richard Harter's version of the 2D Median
// 2007.11.20 Steve Witham sw@remove.this.tiac.net
// This plays with the idea of a multi-dimensional equivalent of
// the median.  
//
// My definition was: First define the "sum of directions" function for any
// point as the sum of unit vectors pointing at all the input points.
// (Right on top of an input point, leave the function undefined.)
// 
// The "median" is a pixel around which the sum-of-directions vectors
// for neighboring pixels "point in opposite directions".  In other words, 
// the sum may not be zero here, but there's either a zero or a 
// discontinuity that cuts through zero, right near here.
//
// Richard Harter's variation of this is: define f(x), where x is a point,
// as the sum of absolute distances from x to all input points.
// A value of x that is a minimum of f() is a candidate median;
// if not unique maybe average the candidates.
//
// I think my sum-of-unit-vectors function is the derivative of Richard's
// sum-of-absolute-values function, so my "solution" and his minimum
// should coincide.
//
// Comments in the code below explain the visuals.
//
// Whether this median is ever non-unique for non-colinear sets of points
// is something I wanted to see with this applet.
// For colinear points it can be a line segment.  E.g. defined this way,
// the median of { 0, 1, 2, 3 } is the segment between 1 and 2 (inclusive).
//
// Another idea similar to a median, from statistics, is the "convex peel,"
// which you get by removing the convex hull points in layers until you 
// have either one point inside (call it the median) 
// or none (in which case you average the innermost hull points).
//
// WEIRDNESS IN THE CODE:
//   I use arrays of two elements to represent points or complex numbers, 
//   and arrays of [width][height][2] for arrays of points or complex numbers,
//   rather than use a point or complex number class.

boolean[][] isPoint;
int [][] pointsXY;
float [][][] dirTotalXY;
float [][] absTotal;  // abs of total of normalized difference vectors
float [][] totalAbs;  // total of abs's of difference vectors
int npointsXY;
int arraySize;
int array_width, zoomX;

final int stateHalt = 0;
final int stateSetup = 1;
final int stateSettingUp = 2;
final int stateDraw = 3;
int drawState = stateSetup;


void mousePressed() {
  drawState = stateSetup;
  frameRate( 4 );
  loop();
}


void setup() {
  arraySize = 250;
  array_width = arraySize;
  zoomX = int( array_width * 1.2 );
  size( zoomX + array_width, arraySize );  // Sets height and width for dims.
  drawState = stateSetup;
  textFont( loadFont( "Helvetica-Bold-12.vlw" ) );
  frameRate( 4 );
  loop();
}

 
void real_setup() 
{ 
  pointsXY = new int[200][2];
  isPoint = new boolean[array_width][height];
  dirTotalXY = new float[array_width][height][2];
  absTotal = new float[array_width][height];
  totalAbs = new float[array_width][height];
  int x, y;
   
  // Pick a set of input points without duplicates
  // (which could be dealt with but would complicate things):
  npointsXY = int( random( 4 ) + 3 );
  for( int i = 0;  i < npointsXY;  i++ ) {
    do {
      x = pointsXY[ i ][0] = int( random(array_width-3) ) + 2;
      y = pointsXY[ i ][1] = int( random(height-3) ) + 2;
    } while( isPoint[ x ][ y ] );  // No two points at the same place.
    isPoint[ x ][ y ] = true;
  }
  loop();
}


void draw() 
{
  int medianXY[];
  int zoom = 5;
  int small = array_width / zoom;
  int left, top;
  
  if( drawState == stateSetup  ||  drawState == stateSettingUp )  { 
    background( #ddddcc );
    text( "Setting up...", width/3, height/3 + 12 );
  }
  if( drawState == stateSetup ) {
    drawState = stateSettingUp;
    real_setup();
    drawState = stateDraw;
    return;
  } else if( drawState == stateSettingUp ) {
    return;
  } else if( drawState != stateDraw ) {
    noLoop();
    return;
  }
  
  float zoomTotalXY[][][] = new float [array_width][height][2];
  float zoomAbsTotal[][] = new float [array_width][height];
  float zoomTotalAbs[][] = new float [array_width][height];
  
  colorMode( HSB, PI * 2.0, 1.0, 1.0 );  
  background( #ddddcc ); 
  
  sumDirections( dirTotalXY, absTotal, totalAbs, 0, 0, 1.0 );
  medianXY = draw_array( pointsXY, dirTotalXY, isPoint, absTotal, 
                          totalAbs, 0, 0, 0, 1.0 );

  stroke( 0, 0, 1 );
  noFill();
  left = medianXY[0] - small/2;
  top =  medianXY[1] - small/2;
  rect( left, top, small, small );
  line( left, top, zoomX, 0 );
  line( left, top + small, zoomX, height - 1 );
  sumDirections( zoomTotalXY, zoomAbsTotal, zoomTotalAbs, left, top, zoom );
  draw_array( pointsXY, zoomTotalXY, isPoint, zoomAbsTotal, 
              zoomTotalAbs, left, top, zoomX, zoom ); 
  drawState = stateHalt;
  noLoop();
}


void sumDirections( float dirTotalXY[][][], float absTotal[][],
                    float totalAbs[][],
                    int x_offset, int y_offset, float mag )
{
  float real, imag, absv;
  int i, j;
  float x, y;
  
  // Calculate the sum-of-directions function for the whole background:
  // Here, ( i, j ) are zoomed coordinates and ( x, y ) are original coordinates.
  for( i = 0;  i < array_width;  i++ ) {
    x = x_offset + i / mag;
    for( j = 0;  j < height;  j++ ) {
      y = y_offset + j / mag;
      real = imag = absv = 0.0;
      for( int k = 0;  k < npointsXY;  k++ ) {
          float d = dist( x, y, pointsXY[k][0], pointsXY[k][1] );
          if( d != 0 ) {
              real += ( pointsXY[k][0] - x ) / d;
              imag += ( pointsXY[k][1] - y ) / d;
          }
          absv += d;
      }
      dirTotalXY[i][j][0] = real;
      dirTotalXY[i][j][1] = imag;
      absTotal[i][j] = dist( 0, 0, real, imag );
      totalAbs[i][j] = absv;
    }
  }
} 


// This colors the picture and finds the median based on
// Richard Harter's min of sum of abs differences method.
int [] draw_array( int pointsXY[][], float dirTotalXY[][][],
                 boolean isPoint[][], float absTotal[][],
                 float totalAbs[][],
                 int x_offset, int y_offset, int draw_x, float mag )
{
  int result [] = new int [2];  // Returns x and y of estimated median.
  float median_x_total = 0.0;
  float median_y_total = 0.0;
  int n_median_guesses = 0;
  float minTotalAbs = -1.0;  // Marker for "unset."
  float contourZoom = 400 * mag / ( height * npointsXY );
  
  // Draw the colored background.
  // Hue represents the direction of the sum-of-directions vector.
  // Darkness bands are contours of the sum-of-lengths.
  noStroke();
  for (int x = 0; x < array_width; x=x+1) { 
    for (int y = 0; y < height; y=y+1) {
        if( minTotalAbs < 0  ||  totalAbs[x][y] < minTotalAbs ) {
          median_x_total = x;
          median_y_total = y;
          minTotalAbs = totalAbs[x][y];
          n_median_guesses = 1;
        } else if ( totalAbs[x][y] == minTotalAbs ) {
          // This case actually happens.
          median_x_total += x;
          median_y_total += y;
          n_median_guesses += 1;
        }
        float real = dirTotalXY[x][y][0];
        float imag = dirTotalXY[x][y][1];
        // Fuzzy contour bands based on totalAbs[][]:
        float strength = .5 + .5 * cos( totalAbs[x][y] * contourZoom );
        fill( atan2( real, imag ) + PI, 0.25 + strength * .5, 
                                        1.00 - strength * .25 );
        rect( draw_x + x, y, 1, 1 );
    } 
  }
  result[ 0 ] = int( .5 + median_x_total / n_median_guesses );
  result[ 1 ] = int( .5 + median_y_total / n_median_guesses );
  
  // Draw the input points:
  for( int k = 0;  k < npointsXY;  k++ ) {
      // Here, ( i, j ) are original coordinates and ( x, y ) are zoomed coordinates.
      int i = pointsXY[k][0];
      int j = pointsXY[k][1];
      int x = int( mag * ( i - x_offset ) );
      int y = int( mag * ( j - y_offset ) );
      if( x < 0  ||  x >= array_width  ||  y < 0  ||  y >= height )  continue;
      
        noStroke(); 
        fill( 0, 0, 0 );
        bigDot( draw_x + x, y );
  }
  
  // Whiten median point:
    noStroke();
    fill( 0, 0, 1 );
    diamond( draw_x + result[0], result[1] );
  return( result );
}


// Is (x,y) where one of the input points is in the zoomed view?
boolean isZoomPoint( int x, int y, int x_offset, int y_offset, float mag )
{
  float fi, fj;
  int i, j;

  // Here, ( i, j ) are original coordinates and ( x, y ) are zoomed coordinates.  
  i = int( fi = x_offset + x / mag );
  j = int( fj = y_offset + y / mag );
  return( i == fi  &&  j == fj  &&  isPoint[i][j] );
}


// The following routines exist because
//    1) fills of circles of a certain size were overflowing their stroke boundaries;
//    2) I wasn't sure of the alignment of centers of various circles (ellipses).


void diamond( int x, int y ) {
  rect( x, y - 2, 1, 5 );
  rect( x-1, y-1, 3, 3 );
  rect( x - 2, y, 5, 1 );
}


void bigDot( int x, int y ) {
  rect( x - 2, y - 3, 5, 7 );
  rect( x - 3, y - 2, 7, 5 );
}


void mediumDot( int x, int y ) {
  rect( x - 1, y - 2, 3, 5 );
  rect( x - 2, y - 1, 5, 3 );
}