"False Leveler"

About

False Leveler is a Processing program I created to do histogram matching with random, artificially-created destination histograms.

Histogram matching is typically used to adjust the colors of one image to those of another, so they can be compared, combined, etc. For example if you take two overlapping photos of a scene and want to stich them together to form a panorama, you'll want them to have the same color profile to help them look like one large image. You can use histogram matching to help do that matching. The interested reader is better off heading to Wikipedia for details rather than relying on whatever half-baked explanation I could cobble together right now.

Rather than having a reference image to match to, I wanted to use an arbitrary and dynamic color distribution. This would allow the color balance of an image to wander randomly over time, giving a look at the myriad possible images that underly and are implied by the one actual photograph taken.

To generate these distributions, I used Perlin noise to create parameters for a Kumaraswamy probability distribution. Kumaraswamy is great for this purpose because it's bounded above and below — at 0 and 1; it's trivial to shift this to [0,255] for use with colors — and it has closed-form, and easy to compute CDF and quantile functions. In addition it produces a fairly wide range of differently shaped PDFs, which gives me the chance to transform my images in to a wider section of color space. For technical information on the Kumaraswamy distribution, see Jones, 2009 [1].

$$ f(x; \alpha, \beta, \delta, \lambda) = \ldots$$
$$\lambda \alpha \beta (\delta+1) x^{\alpha-1}(1-x^\alpha)^{\beta-1}\left(1-(1-x^\alpha)^\beta\right)^{\lambda-1}\left(1-\left(1-(1-x^\alpha)^\beta\right)^\lambda\right)^\delta $$

and CDF by

$$ F(x; \alpha, \beta, \delta, \lambda) = 1-\left(1-\left(1-(1-x^\alpha)^\beta\right)^\lambda\right)^{\delta+1} $$

Gallery

Okay. Enough background. Let's get to some examples. Here are a couple of other samples generated from photos of the Carnival of Venice. Both of these are set to have the colors change slower than Carnival 01 (sup.).

Rust 01 is based on a photo of some oxidized metal scrap, is interesting because unlike the photos above it's impossible to know what the "right" colors are for this image.

This one is generated from an HDR photo, which gives another interesting effect.

Eye 05 is derived from a shot from Rankin's gorgeous "Eyescapes" series. These photos are pretty abstract and alien looking to begin with, and I like how shifting the colors around enhances that.

Finally, here are a couple derived from Monet's "Haystacks at the End of Summer, Morning Effect" (1891) and "San Giorgio Maggiore at Dusk" (1908). I really love the way he would return to the same subjects and paint them again and again in different light conditions. Part of the inspiration for this project was wondering what all the other possible chromatic variations on a single scene could be, and wanting a system that would show them to me.

I think these are probably the least successful of the renders I've tried so far, but I think these and other Monets have good potential as an application for False Leveler.

I would be remiss if I did not note that the beautiful music for all of these comes from Carlos Márquez's recordings of Dirk Massen's piano compositions.

Misc. Discussion

The top image is the photo of Carnivale from the beginning of the page which shows both the original photo's histogram (the solid lines) and the Kumaraswamy distribution that's responsible for this frame (dashed line). Below that is another frame, but this time showing the cumulative histogram from the source (solid lines) and generated image (dotted lines).

One other thing to note is that since I'm using 8 bits per color channel, the matching is not being done by percentile or real-valued quantile, but by 1/256th-ile.

This last bit is non-trivial: it allows me to take advantage of dynamic programming to avoid having to compute the value of the inverse CDF transform for every single pixel. I can memoize the results as I go, so I only need to calculate it a maximum of 256 times per frame instead of once per pixel per frame. Doing that gave me several orders of magnitude speed-up, pushing the frame rate on my test images from something like 2fps to 60+.

You may notice some rendering artifacts when you run this code. It seems depend a lot on the input image. The above videos all used jpegs as input, which certainly doesn't help. Another part of the problem is that I didn't put much effort in to bounding the parameter space. Depending on the whims of my RNG, I might get paramaters for the probability distribution anywhere between $$2^{-2.5}$$ and $$2^{2.5}$$. That works well for most images I've tried but I'm sure I could do a more intelligent check based on the particular mapping function to set bounds on the parameters per input image.

Most of the information you'd need to transform a very dark, unsaturated image in to a very bright, saturated one has been thrown out by the initial encoding and compression of the source. All the invisible-to-the-naked-eye adjustments made by a lossy compression scheme suddenly become very visible when you throw the whole thing halfway across color space.

I suppose the Art World-approved way to cover my ass on something like this would be to claim that the artifacts are left in intentionally as a commentary on the artificiality of my process, but that's bollocks and we all know it. The truth is that I have better and more interesting problems to work on that fiddling with parameters to polish this sort of thing to a high gloss.

Future Work

I'd like to do more with HDR images. They've already got such an atypical appearance that pushing the colors further from normal is easier to get away with — and more dramatic — than when you start with a standard photo. I have one example above with an HDR input, but most of the ones I tried turned out &hellips; weird. I think they need different parameters than normal photos. Maybe a little "flatter" histograms should be generated? I'd have to play around with it.

It would be good to adjust the paramters to the source image, rather than using the same destination histograms regardless of source. I should do an MLE of the KwKw distribution which best fits the source, and then perturb the values from there, rather than mapping every image to .

I could also adjust the rate of change on everything based on the actual frameRate so that a live version would be more consistent from one machine to another.

I'm imagining this displayed in a lobby. It looks like a photograph, but if you sit there long enough you notice the hues are different over time, but if you stare at it you can't see the changes at the margin. If I was really ambitious I would crack open a psychophysics text and figure out what the limen for color change is, then back-solve to determine appropriate parameters to put the output right around the threshold of static-vs-dyanmic. In reality, I'm going to settle for cranking down some of the step sizes to about 0.000001 and that should do the trick.

I've messed about a little with doing this in LAB rather than RGB color space. It worked, but the histograms of photos look totally different in LAB, so the paramaters would have to change drastically, and I haven't gotten around to finding ones that work. I suppose this would be a non-issue if I had that MLE capability I mentioned...

Finally, you could also use False Leveler to false-color B&W/intensity images, but I haven't played around with that much. There's still some cruft left over in the code below from when I tested that a bit.

Code

Here's the code. It's more than a bit messy. Sorry. I figured better to get it posted publicly-but-dirty than not at all. You can download a tarball of it here.

Each color channel gets four paramaters to define it's own KwKw distribution. aParam corresponds to $$\alpha^{-1}$$,  bParam to $$\beta^{-1}$$,  etc. The inverses there are used because that's what's needed for the quantile function, which is being calculated way more often than the PDF or CDF themselves.

Increasing aStep (&c.) will make the colors change more quickly, and vice-versa. The only other knob that you would want to turn immediately is expon, in the setup() function. Higher values will allow more dramatically skewed distribution. That's about all, so go ahead and dive in.

/****************************************************************
 * False Leveler
 * Artificial Histogram Matching
 * Jared Sylvester, 2013
 *
 * www.jsylvest.com | jared@jsylvest.com | @jsylvest
 * vimeo.com/jsylvest/videos
 *
 * Released under GPL v3 License.
 *   Feel free to use and modify this code as you wish
 *   but be a mensch and throw a little credit my way
 *
 * This code is very much still a work in progress.
 *   I take no liability, implicit or otherwise, etc. etc.
 ****************************************************************/

/* For information on the generalized Kw distribution:
 http://arxiv.org/pdf/1004.0911v1.pdf
 NB: 4-param version with closed form CDF & inverse CDF */

String inputImage = "eye5.jpg";

final int DISTRO_UNIFORM   = 1001;
final int DISTRO_KUMARA    = 1002; // Kumaraswamy http://en.wikipedia.org/wiki/Kumaraswamy_distribution
final int DISTRO_GKUMARA4P = 1003;
int distro = DISTRO_GKUMARA4P;

boolean saveFrames = true;

boolean useLAB = false;

PImage img;

int histSize;
float[] srcHistR;
float[] srcHistG;
float[] srcHistB;
float[] srcHist;
int n;

float[] srcHistL;
float[] srcHistA;

float[] aParam = new float[4];
float[] bParam = new float[4];
float[] cParam = new float[4];
float[] dParam = new float[4];
float[] aOff = new float[4];
float[] bOff = new float[4];
float[] cOff = new float[4];
float[] dOff = new float[4];
float[] aStep = new float[4];
float[] bStep = new float[4];
float[] cStep = new float[4];
float[] dStep = new float[4];

boolean showHisto = false;
boolean colorHisto = true;

boolean bwSrc = false;
boolean bwNewMethod = true;

boolean smoothPDFs = true;

float[][] memo;

boolean showOrigHisto = false;
boolean showDestHisto = false;
boolean showOrigPDF =  true;
boolean showDestPDF =  true;

double[][] srcPDF;

float minVal =    0;
float maxVal =  110;

void setup() {
  //noiseSeed(0);
  //randomSeed(0);
  img = loadImage(inputImage);
  size(img.width, img.height);
  n = img.width*img.height;
  frameRate(60);
  //noLoop();
  noiseDetail(5);
  ellipseMode(RADIUS);
  if ( !useLAB ) {
    histSize = 256;
    srcHistR = new float[histSize];
    srcHistG = new float[histSize];
    srcHistB = new float[histSize];
    srcHist = new float[histSize];
  } 
  else {
    histSize = 300;
    srcHistL = new float[histSize];
    srcHistA = new float[histSize];
    srcHistB = new float[histSize];
  }
  memo = new float[histSize][4];
  println(nfc(n)+" pixels");
  loadPixels();
  img.loadPixels();
  float r = 0, g = 0, b = 0, v = 0;
  srcPDF = new double[histSize][3];

  if ( useLAB ) {
    float[] lab = new float[3];
    int labL, labA, labB;
    for ( int i = 0; i < n; i++ ) {
      r = red(img.pixels[i]);
      g = green(img.pixels[i]);
      b = blue(img.pixels[i]);
      lab = rgb2xyz(r, g, b);
      labL = round(map(lab[0], minVal, maxVal-12, 0, histSize));
      labA = round(map(lab[1], minVal, maxVal-8, 0, histSize));
      labB = round(map(lab[2], minVal, maxVal, 0, histSize));
      for ( int j = labL; j < histSize; j++ ) { 
        srcHistL[j] += 1;
      }
      for ( int j = labA; j < histSize; j++ ) { 
        srcHistA[j] += 1;
      }
      for ( int j = labB; j < histSize; j++ ) { 
        srcHistB[j] += 1;
      }
      srcPDF[labL][0] += 1;
      srcPDF[labA][1] += 1;
      srcPDF[labB][2] += 1;
      if ( smoothPDFs ) {
        srcPDF[constrain(round(labL-1), 0, histSize-1)][0] += 1.0;
        srcPDF[constrain(round(labL+1), 0, histSize-1)][0] += 1.0;
        srcPDF[constrain(round(labA-1), 0, histSize-1)][1] += 1.0;
        srcPDF[constrain(round(labA+1), 0, histSize-1)][1] += 1.0;
        srcPDF[constrain(round(labB-1), 0, histSize-1)][2] += 1.0;
        srcPDF[constrain(round(labB+1), 0, histSize-1)][2] += 1.0;
      }
    }
  } 
  else {

    for ( int i = 0; i < n; i++ ) {
      r = red(img.pixels[i]);
      g = green(img.pixels[i]);
      b = blue(img.pixels[i]);
      if ( bwNewMethod && bwSrc ) {
        r = red(colormap(red(img.pixels[i])/255.0));
        g = green(colormap(green(img.pixels[i])/255.0));
        b = blue(colormap(blue(img.pixels[i])/255.0));
      }
      for ( int j = int(r); j < histSize; j++ ) {
        srcHistR[j] += 1;
      }
      for ( int j = int(g); j < histSize; j++ ) {
        srcHistG[j] += 1;
      }
      for ( int j = int(b); j < histSize; j++ ) {
        srcHistB[j] += 1;
      }
      v = round((r+g+b)/3.0);
      for ( int j = int(v); j < histSize; j++ ) {
        srcHist[j] += 1;
      }
      srcPDF[round(r)][0] += 1;
      srcPDF[round(g)][1] += 1;
      srcPDF[round(b)][2] += 1;
      if ( smoothPDFs ) {
        srcPDF[constrain(round(r-1), 0, histSize-1)][0] += 1.0;
        srcPDF[constrain(round(r+1), 0, histSize-1)][0] += 1.0;
        srcPDF[constrain(round(g-1), 0, histSize-1)][1] += 1.0;
        srcPDF[constrain(round(g+1), 0, histSize-1)][1] += 1.0;
        srcPDF[constrain(round(b-1), 0, histSize-1)][2] += 1.0;
        srcPDF[constrain(round(b+1), 0, histSize-1)][2] += 1.0;
      }
    }
  }

  double divisor = (double)(n);
  if ( smoothPDFs ) {
    divisor = (double)(3.0*n);
  }
  for ( int i = 0; i < histSize; i++ ) {
    if ( useLAB ) {
      srcHistL[i] = srcHistL[i] / float(n);
      srcHistA[i] = srcHistL[i] / float(n);
      srcHistB[i] = srcHistL[i] / float(n);
    } 
    else {
      srcHistR[i] = srcHistR[i] / float(n);
      srcHistG[i] = srcHistG[i] / float(n);
      srcHistB[i] = srcHistB[i] / float(n);
      srcHist[i]  = srcHist[i]  / float(n);
    }
    srcPDF[i][0] = srcPDF[i][0] / divisor;
    srcPDF[i][1] = srcPDF[i][1] / divisor;
    srcPDF[i][2] = srcPDF[i][2] / divisor;
  }
  for ( int i = 0; i < 4; i++ ) {
    aParam[i] = 1.0;
    bParam[i] = 1.0;
    cParam[i] = 1.0;
    dParam[i] = 1.0;
    aOff[i] = random(0.0, 1000.0);
    bOff[i] = random(0.0, 1000.0);
    cOff[i] = random(0.0, 1000.0);
    dOff[i] = random(0.0, 1000.0);
    /* increase these values to make colors shift faster & vice-versa */
    aStep[i] = 0.005;
    bStep[i] = 0.005; 
    cStep[i] = 0.005; 
    dStep[i] = 0.005;
  }
}

void draw() {
  update();
  float v1, v2;
  float r1 = 0, g1 = 0, b1 = 0;
  float r2 = 0, g2 = 0, b2 = 0;
  int labL1, labA1, labB1;
  float labL2, labA2, labB2;
  float[] lab = new float[3];
  float[] rgb = new float[3];

  color c1 = color(0, 0, 0);
  color c2 = color(0, 0, 0);
  for ( int i = 0; i < histSize; i++ ) {
    memo[i][0] = -1;
    memo[i][1] = -1;
    memo[i][2] = -1;
    memo[i][3] = -1;
  }
  for ( int x = 0; x < width; x++ ) {
    for ( int y = 0; y < height; y++ ) {
      c1 = img.pixels[y*width+x];
      if ( !useLAB ) {
        if ( !bwNewMethod && bwSrc ) {
          v1 = (red(c1)+green(c1)+blue(c1))/3.0;
          if ( memo[round(v1)][4] == -1 ) {
            v2 = transform(v1, srcHist, aParam[3], bParam[3], cParam[3], dParam[3]);
            memo[round(v1)][4] = v2;
          } 
          else {
            v2 = memo[round(v1)][4];
          }
          c2 = color(v2);
        } 
        else {
          if ( bwSrc ) {
            c1 = colormap((red(c1)+green(c1)+blue(c1))/(3.0*255.0));
          }
          r1 = red(c1);   // c1 >> 16 & 0xFF
          g1 = green(c1); // c1 >> 8 & 0xFF
          b1 = blue(c1);  // c1 & 0xFF

          if ( memo[round(r1)][0] == -1 ) {
            r2 = transform(r1, srcHistR, aParam[0], bParam[0], cParam[0], dParam[0]);
            memo[round(r1)][0] = r2;
          } 
          else {
            r2 = memo[round(r1)][0];
          }

          if ( memo[round(g1)][1] == -1 ) {
            g2 = transform(g1, srcHistG, aParam[1], bParam[1], cParam[1], dParam[1]);
            memo[round(g1)][1] = g2;
          } 
          else {
            g2 = memo[round(g1)][1];
          }

          if ( memo[round(b1)][2] == -1 ) {
            b2 = transform(b1, srcHistB, aParam[2], bParam[2], cParam[2], dParam[2]);
            memo[round(b1)][2] = b2;
          } 
          else {
            b2 = memo[round(b1)][2];
          }

          c2 = color(r2, g2, b2);
        }
      } 
      else {
        lab = rgb2xyz( red(c1), green(c1), blue(c1) );
        labL1 = round(map(lab[0], minVal, maxVal-12, 0, histSize));
        labA1 = round(map(lab[1], minVal, maxVal-8, 0, histSize));
        labB1 = round(map(lab[2], minVal, maxVal, 0, histSize));
        if ( memo[labL1][0] == -1 ) {
          labL2 = transform(labL1, srcHistL, aParam[0], bParam[0], cParam[0], dParam[0]);
          memo[labL1][0] = labL2;
        } 
        else {
          labL2 = memo[labL1][0];
        }
        if ( memo[labA1][1] == -1 ) {
          labA2 = transform(labA1, srcHistL, aParam[0], bParam[0], cParam[0], dParam[0]);
          memo[labA1][1] = labA2;
        } 
        else {
          labA2 = memo[labA1][1];
        }
        if ( memo[labB1][2] == -1 ) {
          labB2 = transform(labB1, srcHistL, aParam[0], bParam[0], cParam[0], dParam[0]);
          memo[labB1][2] = labB2;
        } 
        else {
          labB2 = memo[labB1][2];
        }
        rgb = xyz2rgb(labL2, labA2, labB2);
        c2 = color(rgb[0], rgb[1], rgb[2]);
      }
      pixels[y*width+x] = c2;
    }
  }
  if (mousePressed == true) {
    image(img, 0, 0);
  } 
  else {
    updatePixels();
  }
  ///println(red(c1)+", "+green(c1)+", "+blue(c1));
  ///println(srcHistR[round(red(c1))]);

  if ( saveFrames ) {
    saveFrame("out/frame-####.jpg");
  }
  if ( showHisto ) {
    if ( showOrigHisto ) { 
      plotOrigHisto();
    }
    if ( showDestHisto ) { 
      plotDestHisto();
    }
    if ( showOrigPDF ) { 
      plotOrigPDF();
    }
    if ( showDestPDF ) { 
      plotDestPDF();
    }
  }
  fill(0);
  text(nf(frameRate, 2, 1), 10, height-20);
}

float transform(float src, float[] hist, float a, float b, float c, float d) {
  float result = hist[round(src)];
  if ( distro == DISTRO_UNIFORM ) {
    // do nothing
  } 
  else if ( distro == DISTRO_KUMARA ) {
    result = pow( 1.0 - pow(1.0-result, b), a );
  } 
  else if ( distro == DISTRO_GKUMARA4P ) {
    result = pow(1.0-pow(1.0-pow(1.0-pow(1.0-result, 1.0/(d+1)), c), b), a);
  }
  result = result * 255.0;
  return result;
}

// http://www.sandia.gov/~kmorel/documents/ColorMaps/
color c0 = color( 0.217*255, 0.525*255, 0.910*255 );
color c1 = color( 0.865*255, 0.865*255, 0.865*255 );
color c2 = color( 0.677*255, 0.492*255, 0.093*255 );

color colormap(float src) {
  color res;
  if ( src < 0.5 ) {
    res = lerpColor(c0, c1, map(src, 0.0, 0.5, 0.0, 1.0));
  } 
  else {
    res = lerpColor(c1, c2, map(src, 0.5, 1.0, 0.0, 1.0));
  }
  return res;
}

void update() {
  if ( distro == DISTRO_KUMARA ) {
    float expon = 4.0;
    for ( int i = 0; i < 4; i++ ) {
      aOff[i] = aOff[i] + aStep[i];
      bOff[i] = bOff[i] + bStep[i];
      if ( !useLAB ) {
        aParam[i] = pow(2.0, expon*( noise(aOff[i]) - 0.5 ));
        bParam[i] = pow(2.0, expon*( noise(bOff[i]) - 0.5 ));
      } 
      else {
        aParam[i] = 1.0+noise(aOff[i]);
        bParam[i] = pow(2.0, 8.0*( noise(bOff[i]) - 0.75 ));
        /*
        aParam[i] = 1.0/map(noise(aOff[i]), 0.0, 1.0, 1.7, 2.3);
         bParam[i] = 1.0/map(noise(bOff[i]), 0.0, 1.0, 1.7, 6.3);
         */
      }
    }
  }
  if ( distro == DISTRO_GKUMARA4P ) {
    /* increase expon to allow more extreme distribution shapes */
    float expon = 2.5;
    for ( int i = 0; i < 4; i++ ) {
      aOff[i] = aOff[i] + aStep[i];
      bOff[i] = bOff[i] + bStep[i];
      cOff[i] = cOff[i] + cStep[i];
      dOff[i] = dOff[i] + dStep[i];
      aParam[i] = pow(2.0, expon*( noise(aOff[i]) - 0.5 ));
      bParam[i] = pow(2.0, expon*( noise(bOff[i]) - 0.5 ));
      cParam[i] = pow(2.0, expon*( noise(cOff[i]) - 0.5 ));
      dParam[i] = pow(2.0, expon*( noise(dOff[i]) - 0.5 ));
    }
  }
}

void keyPressed() {
  if ( key == 't' ) {
    showHisto = !showHisto;
  } 
  else if ( key == 'T' ) {
    colorHisto = !colorHisto;
  }
}

These functions are just for drawing the histogram curves (for example in the figure at the top of this page). They're very handy for debugging/understanding but by default are off. Pressing t while the script is running will toggle them on. The default is to show the PDF view; pressing y will switch to CDF.

void plotOrigHisto() {
  pushStyle();
  stroke(0);
  strokeWeight(2);
  for( int i = 0; i < 256; i += 1 ) {
    if( colorHisto ) {
      stroke(255,0,0);
      line( map(i, 0, 255, 0, width), map(srcHistR[i], 0, 1.0, height, 0), map(i+1, 0, 255, 0, width), map(srcHistR[i+1], 0, 1.0, height, 0) );
      stroke(0,255,0);
      line( map(i, 0, 255, 0, width), map(srcHistG[i], 0, 1.0, height, 0), map(i+1, 0, 255, 0, width), map(srcHistG[i+1], 0, 1.0, height, 0) );
      stroke(0,0,255);
      line( map(i, 0, 255, 0, width), map(srcHistB[i], 0, 1.0, height, 0), map(i+1, 0, 255, 0, width), map(srcHistB[i+1], 0, 1.0, height, 0) );
    } else {
      stroke(0);
      line( map(i, 0, 255, 0, width), map(srcHist[i], 0, 1.0, height, 0), map(i+1, 0, 255, 0, width), map(srcHist[i+1], 0, 1.0, height, 0) );
    }
  }
  popStyle();
}

void plotDestHisto() {
  pushStyle();
  strokeWeight(1);
  stroke(0);
  float r, g, b;
  float aP, bP, v;
  /*
  // draw markers equally across height
  for( float i = 0; i <= 1.0; i += 0.025 ) {
    if( colorHisto ) {
      fill(255,0,0);
      r = pow( 1.0 - pow(1.0-i, bParam[0]), aParam[0] );
      ellipse( map(r, 0, 1.0, 0, width), map(i, 0, 1.0, height, 0), 3, 3 );
      fill(0,255,0);
      g = pow( 1.0 - pow(1.0-i, bParam[1]), aParam[1] );
      ellipse( map(g, 0, 1.0, 0, width), map(i, 0, 1.0, height, 0), 3, 3 );
      fill(0,0,255);
      b = pow( 1.0 - pow(1.0-i, bParam[2]), aParam[2] );
      ellipse( map(b, 0, 1.0, 0, width), map(i, 0, 1.0, height, 0), 3, 3 );
    } else {
      aP = (aParam[0]+aParam[1]+aParam[2])/3.0;
      bP = (bParam[0]+bParam[1]+bParam[2])/3.0;
      fill(255);
      v = pow( 1.0 - pow(1.0-i, bP), aP );
      ellipse( map(v, 0, 1.0, 0, width), map(i, 0, 1.0, height, 0), 3, 3 );
    }
  }
  */
  // draw markers equally across width
  for( float i = 0; i <= 1.0; i += 0.025 ) {
    if( colorHisto ) {
      fill(255,0,0);
      r = 1.0 - pow( 1.0 - pow(i,1/aParam[0]), 1/bParam[0]);
      ellipse( map(i, 0, 1.0, 0, width)-3, map(r, 0, 1.0, height, 0), 3, 3 );
      fill(0,255,0);
      g = 1.0 - pow( 1.0 - pow(i,1/aParam[1]), 1/bParam[1]);
      ellipse( map(i, 0, 1.0, 0, width), map(g, 0, 1.0, height, 0), 3, 3 );
      fill(0,0,255);
      b = 1.0 - pow( 1.0 - pow(i,1/aParam[2]), 1/bParam[2]);
      ellipse( map(i, 0, 1.0, 0, width)+3, map(b, 0, 1.0, height, 0), 3, 3 );
    } else {
      aP = (aParam[0]+aParam[1]+aParam[2])/3.0;
      bP = (bParam[0]+bParam[1]+bParam[2])/3.0;
      fill(255);
      v = 1.0 - pow( 1.0 - pow(i,1/aP), 1/bP);
      ellipse( map(i, 0, 1.0, 0, width), map(v, 0, 1.0, height, 0), 3, 3 );
    }
  }
  popStyle();
}

void plotOrigPDF() {
  pushStyle();
  float x1 = 0, y1 = 0, x2 = 0, y2 = 0;
  float s = 0.01;
  if( useLAB ) {
    s = 0.01;
  }
  strokeWeight(2);
  for( int i = 0; i < histSize-1; i += 1 ) {
    x1 = map(i,0,histSize-1,0,width);
    x2 = map(i+1,0,histSize-1,0,width);
    y1 = map((float)(srcPDF[i][0]), 0, s, height, 0);
    y2 = map((float)(srcPDF[i+1][0]), 0, s, height, 0);
    stroke(255,0,0);
    line(x1,y1,x2,y2);
    y1 = map((float)(srcPDF[i][1]), 0, s, height, 0);
    y2 = map((float)(srcPDF[i+1][1]), 0, s, height, 0);
    stroke(0,200,0);
    line(x1,y1,x2,y2);
    y1 = map((float)(srcPDF[i][2]), 0, s, height, 0);
    y2 = map((float)(srcPDF[i+1][2]), 0, s, height, 0);
    stroke(0,0,255);
    line(x1,y1,x2,y2);
  }
  popStyle();
}

void plotDestPDF() {
  pushStyle();
  strokeCap(SQUARE);
  float a = 1.0;
  float b = 1.0;
  float c = 1.0;
  float d = 1.0;
  float x1 = 0, y1 = 0, x2 = 0, y2 = 0;
  strokeWeight(2);
  float pdfMax = 4.0;
  float step = 0.01;
  for( float i = step/2.0; i < 1.0-step; i += step ) {
    x1 = i;
    x2 = i + step*0.67;
    if( colorHisto ) {
      stroke(192, 0, 0);
      a = 1.0/aParam[0];
      b = 1.0/bParam[0];
      if( distro == DISTRO_KUMARA ) {
        y1 = min( a*b*pow(x1,a-1.0)*pow(1-pow(x1,a),b-1.0), 4000 );
        y2 = min( a*b*pow(x2,a-1.0)*pow(1-pow(x2,a),b-1.0), 4000 );
      } else if( distro == DISTRO_GKUMARA4P ) {
        c = 1.0/cParam[0];
        d = dParam[0];
        y1 = min( a*b*c*(d+1)*pow(x1,a-1.0)*pow(1-pow(x1,a),b-1.0)*pow(1-pow(1-pow(x1,a),b),c-1)*pow(1-pow(1-pow(1-pow(x1,a),b),c),d), 4000 );
        y2 = min( a*b*c*(d+1)*pow(x2,a-1.0)*pow(1-pow(x2,a),b-1.0)*pow(1-pow(1-pow(x2,a),b),c-1)*pow(1-pow(1-pow(1-pow(x2,a),b),c),d), 4000 );
      }
      line( map(x1, 0, 1.0, 0, width), map(y1, 0, pdfMax, height, 0), map(x2, 0, 1.0, 0, width), map(y2, 0, pdfMax, height, 0) );
      stroke(0, 192, 0);
      a = 1.0/aParam[1];
      b = 1.0/bParam[1];
      if( distro == DISTRO_KUMARA ) {
        y1 = min( a*b*pow(x1,a-1.0)*pow(1-pow(x1,a),b-1.0), 4000 );
        y2 = min( a*b*pow(x2,a-1.0)*pow(1-pow(x2,a),b-1.0), 4000 );
      } else if( distro == DISTRO_GKUMARA4P ) {
        c = 1.0/cParam[1];
        d = dParam[1];
        y1 = min( a*b*c*(d+1)*pow(x1,a-1.0)*pow(1-pow(x1,a),b-1.0)*pow(1-pow(1-pow(x1,a),b),c-1)*pow(1-pow(1-pow(1-pow(x1,a),b),c),d), 4000 );
        y2 = min( a*b*c*(d+1)*pow(x2,a-1.0)*pow(1-pow(x2,a),b-1.0)*pow(1-pow(1-pow(x2,a),b),c-1)*pow(1-pow(1-pow(1-pow(x2,a),b),c),d), 4000 );
      }
      line( map(x1, 0, 1.0, 0, width), map(y1, 0, pdfMax, height, 0), map(x2, 0, 1.0, 0, width), map(y2, 0, pdfMax, height, 0) );
      stroke(0, 0, 192);
      a = 1.0/aParam[2];
      b = 1.0/bParam[2];
      if( distro == DISTRO_KUMARA ) {
        y1 = min( a*b*pow(x1,a-1.0)*pow(1-pow(x1,a),b-1.0), 4000 );
        y2 = min( a*b*pow(x2,a-1.0)*pow(1-pow(x2,a),b-1.0), 4000 );
      } else if( distro == DISTRO_GKUMARA4P ) {
        c = 1.0/cParam[2];
        d = dParam[2];
        y1 = min( a*b*c*(d+1)*pow(x1,a-1.0)*pow(1-pow(x1,a),b-1.0)*pow(1-pow(1-pow(x1,a),b),c-1)*pow(1-pow(1-pow(1-pow(x1,a),b),c),d), 4000 );
        y2 = min( a*b*c*(d+1)*pow(x2,a-1.0)*pow(1-pow(x2,a),b-1.0)*pow(1-pow(1-pow(x2,a),b),c-1)*pow(1-pow(1-pow(1-pow(x2,a),b),c),d), 4000 );
      }
      line( map(x1, 0, 1.0, 0, width), map(y1, 0, pdfMax, height, 0), map(x2, 0, 1.0, 0, width), map(y2, 0, pdfMax, height, 0) );
    } else { 
      stroke(32);
      a = (1.0/(aParam[0]+aParam[1]+aParam[2]))/3.0;
      b = (1.0/(bParam[0]+bParam[1]+bParam[2]))/3.0;
      y1 = min( a*b*pow(x1,a-1.0)*pow(1-pow(x1,a),b-1.0), 4000 );
      y2 = min( a*b*pow(x2,a-1.0)*pow(1-pow(x2,a),b-1.0), 4000 );
      line( map(x1, 0, 1.0, 0, width), map(y1, 0, pdfMax, height, 0), map(x2, 0, 1.0, 0, width), map(y2, 0, pdfMax, height, 0) );
    }
  }
  popStyle();
}

This is a set of function to convert between RGB and LAB colorspaces. I deactivated the LAB feature for the time being, so this isn't really even necessary. But I want to return to that as an option in the future, so it's staying in.

float[] rgb2xyz( float R, float G, float B ) {
  // R, G, B all in 0 to 255

  float r = R/255.0;
  float g = G/255.0;
  float b = B/255.0;

  if (r > 0.04045) { 
    r = pow((r + 0.055) / 1.055, 2.4);
  }
  else { 
    r = r / 12.92;
  }
  if ( g > 0.04045) { 
    g = pow((g + 0.055) / 1.055, 2.4);
  }
  else { 
    g = g / 12.92;
  }
  if (b > 0.04045) { 
    b = pow((b + 0.055) / 1.055, 2.4);
  }
  else {  
    b = b / 12.92;
  }

  r = r * 100.0;
  g = g * 100.0;
  b = b * 100.0;

  //Observer. = 2°, Illuminant = D65
  float xyz[] = new float[3];
  xyz[0] = r * 0.4124 + g * 0.3576 + b * 0.1805;
  xyz[1] = r * 0.2126 + g * 0.7152 + b * 0.0722;
  xyz[2] = r * 0.0193 + g * 0.1192 + b * 0.9505;
  return xyz;
}

float[] xyz2lab( float X, float Y, float Z ) {
  final float REF_X = 95.047; // Observer= 2°, Illuminant= D65
  final float REF_Y = 100.000; 
  final float REF_Z = 108.883;
  
  float x = X / REF_X;
  float y = Y / REF_Y;
  float z = Z / REF_Z;
  
  float oneThird = 1.0/3.0;
  if ( x > 0.008856 ) { x = pow( x , oneThird ); }
  else { x = ( 7.787 * x ) + ( 16/116 ); }
  if ( y > 0.008856 ) { y = pow( y , oneThird ); }
  else { y = ( 7.787 * y ) + ( 16/116 ); }
  if ( z > 0.008856 ) { z = pow( z , oneThird ); }
  else { z = ( 7.787 * z ) + ( 16/116 ); }
  
  float lab[] = new float[3];
  lab[0] = ( 116.0 * y ) - 16.0;
  lab[1] = 500.0 * ( x - y );
  lab[2] = 200.0 * ( y - z );
  
  return lab;
}

float[] lab2xyz( float l, float a, float b ) {
  final float REF_X = 95.047; // Observer= 2°, Illuminant= D65
  final float REF_Y = 100.000; 
  final float REF_Z = 108.883;
  
  float y = (l + 16.0) / 116.0;
  float x = a / 500.0 + y;
  float z = y - b / 200.0;
  
  if ( pow( y , 3 ) > 0.008856 ) { y = pow( y , 3 ); }
  else { y = ( y - 16.0 / 116.0 ) / 7.787; }
  if ( pow( x , 3 ) > 0.008856 ) { x = pow( x , 3 ); }
  else { x = ( x - 16.0 / 116.0 ) / 7.787; }
  if ( pow( z , 3 ) > 0.008856 ) { z = pow( z , 3 ); }
  else { z = ( z - 16.0 / 116.0 ) / 7.787; }
 
  float[] xyz = new float[3];
  xyz[0] = REF_X * x;     
  xyz[1] = REF_Y * y; 
  xyz[2] = REF_Z * z; 
 
  return xyz;
}

float[] xyz2rgb( float X, float Y, float Z ) {
  //X from 0 to  95.047    (Observer = 2°, Illuminant = D65)
  //Y from 0 to 100.000
  //Z from 0 to 108.883
    
  float x = X / 100.0;      
  float y = Y / 100.0;      
  float z = Z / 100.0; 
  
  float r = x * 3.2406 + y * -1.5372 + z * -0.4986;
  float g = x * -0.9689 + y * 1.8758 + z * 0.0415;
  float b = x * 0.0557 + y * -0.2040 + z * 1.0570;
  
  if ( r > 0.0031308 ) { r = 1.055 * pow( r , ( 1.0 / 2.4 ) ) - 0.055; }
  else { r = 12.92 * r; }
  if ( g > 0.0031308 ) { g = 1.055 * pow( g , ( 1.0 / 2.4 ) ) - 0.055; }
  else { g = 12.92 * g; }
  if ( b > 0.0031308 ) { b = 1.055 * pow( b , ( 1.0 / 2.4 ) ) - 0.055; }
  else { b = 12.92 * b; }
  
  float[] rgb = new float[3];
  rgb[0] = round( r*255.0 );
  rgb[1] = round( g*255.0 );
  rgb[2] = round( b*255.0 );
  
  return rgb;
}

float[] rgb2lab( float R, float G, float B ) {
  float[] xyz = new float[3];
  xyz = rgb2xyz(R,G,B);
  float[] lab = new float[3];
  lab = xyz2lab(xyz[0], xyz[1], xyz[2]);
  return lab;
}

float[] lab2rgb( float l, float a, float b ) {
  float[] xyz = new float[3];
  xyz = lab2xyz(l,a,b);
  float[] rgb = new float[3];
  rgb = xyz2rgb(xyz[0], xyz[1], xyz[2]);
  return rgb;
}