Getting the average HSB from pixels

Thank you Jeremy! I will reply in chunks because there are a number of unrelated questions.
(1) The Rosetta Code examples are very interesting. BUT the Processing contribution is wrong. Some of the other contributions are also wrong. As the introduction says, the average of 350 and 10 should be zero. The Processing example in Rosetta Code gives -7.80.
How does that website work? Should we replace the current version?

Mean angle from Rosetta Code

Wow! Just one small thing. Do you need to add an extra line of code saying:
a = (a < 0) ? a += 360 : 0;

a is the answer to meanHSB[0], etc.

1 Like

I could not agree more. I am using this in a larger sketch. A sample is taken from an area (ā€œneighbourhood blockā€) and then averaged. This will be the code for that averaging.

1 Like

Yes, it is an open wiki ā€“ if you find errors or problems, just edit that section (not the whole page) and replace them!

Funny, I didnā€™t see a Processing version when I looked earlier, only Java. I wonder if I missed it, or if it is new. Edit: looks like it was added May 4 2021, so I missed it.

Can I post your code there? Or do you want to?

Good catch. There is a meaningful negative angle, but not a meaningful negative hue.

Iā€™ll modify it above.

I was going to say ā€œyesā€ā€¦ butā€¦

After taking another look at the MeanAngle Rosetta page, mean angle is a more complex problem than I thought. The example questions get at common corner cases on either general floating point precision or the behavior of common atan2() implementations.

The most common answers given are 0.0, -90.0, and 20.0 ā€“ but even though most languages agree, the second answer is clearly a strange artifact, when a more meaningful answer would be either 0 or NaN. A comment on the XPL0 answer explains:

The second example is interesting because it computes ATan2(0.,0.), which is undefined in mathematics (like division by zero), but the floating point processor in IBM-PC-type computers defines it as zero. The reason -90 gets displayed instead is due to small rounding errors (and another extra-mathematical concept, -0). In fact almost any angle can result because of slight rounding errors as Y and X both approach zero.

Some of the answers special-case any very small values as equivalent to atan2(0,0) and manually return something else. But it seems like the Julia answer gets it right ā€“ in Julia they implement two built-in methods, sind() and cosd(), which use a piecewise design to return good sin and cos values for degrees:

In particular this eliminates floating point error on common values such as:

  • sin(180) degrees: 0.0, not -8.742278E-8
  • sin(360) degrees: 0.0, not 1.7484555E-7

I feel like the need for a good piece-wise solution comes up over and over with winding. Anywayā€¦ I did a quick line-by-line translation of the Julia ā†’ processing:

// function sind(x::Real)
   float sind(float x) {
//   if isinf(x)
//     return throw(DomainError(x, "`x` cannot be infinite."))
     if (x==Float.POSITIVE_INFINITY || x==Float.NEGATIVE_INFINITY) {
       throw new IllegalArgumentException("x cannot be infinite");
     }
//   elseif isnan(x)
//     return oftype(x,NaN)
//   end
     else if (Float.isNaN(x)) {
       throw new IllegalArgumentException("x cannot be NaN");
     }
//   rx = copysign(float(rem(x,360)),x)
     float rx = x % 360;
//   arx = abs(rx)
     float arx = abs(rx);
//   if rx == zero(rx)
//     return rx
     if (rx == 0.0 || rx == -0.0 ) {
       return rx;
     }
//   elseif arx < oftype(rx,45)
//     return sin_kernel(deg2rad_ext(rx))
     else if (arx < 45) {
       return sin(radians(rx));
     }
//   elseif arx <= oftype(rx,135)
//     y = deg2rad_ext(oftype(rx,90) - arx)
//     return copysign(cos_kernel(y),rx)
     else if (arx <= 135) {
       float y = abs(cos(radians(90 - arx)));
       if (rx < 0) return -y;
       return y;
     }
//   elseif arx == oftype(rx,180)
//     return copysign(zero(rx),rx)
     else if (arx == 180) {
       if (rx < 0) return -0.0;
       return 0.0;
     }
//   elseif arx < oftype(rx,225)
//     y = deg2rad_ext((oftype(rx,180) - arx)*sign(rx))
//     return sin_kernel(y)
     else if (arx < 255) {
       float y = radians(180.0 - arx);
       if(rx < 0) y = y * -1;
       return sin(y);
     }
//   elseif arx <= oftype(rx,315)
//     y = deg2rad_ext(oftype(rx,270) - arx)
//     return -copysign(cos_kernel(y),rx)
     else if (arx <= 315) {
       float y = abs(cos(radians(270.0 - arx)));
       if (rx < 0) return y;
       return -y;
     }
//   else
//     y = deg2rad_ext(rx - copysign(oftype(rx,360),rx))
//     return sin_kernel(y)
//   end
     else {
       float y;
       if (rx < 0) y = radians(rx - 360);
       else y = radians(rx + 360);
       return sin(y);
     }
   }
// end

ā€¦and the final Processing function looks like this:

float sind(float x) {
  if (x==Float.POSITIVE_INFINITY || x==Float.NEGATIVE_INFINITY) {
    throw new IllegalArgumentException("x cannot be infinite");
  }
  else if (Float.isNaN(x)) {
    throw new IllegalArgumentException("x cannot be NaN");
  }
  float rx = x % 360;
  float arx = abs(rx);
  if (rx == 0.0 || rx == -0.0 ) {
    return rx;
  }
  else if (arx < 45) {
    return sin(radians(rx));
  }
  else if (arx <= 135) {
    float y = abs(cos(radians(90 - arx)));
    if (rx < 0) return -y;
    return y;
  }
  else if (arx == 180) {
    if (rx < 0) return -0.0;
    return 0.0;
  }
  else if (arx < 255) {
    float y = radians(180.0 - arx);
    if(rx < 0) y = y * -1;
    return sin(y);
  }
  else if (arx <= 315) {
    float y = abs(cos(radians(270.0 - arx)));
    if (rx < 0) return y;
    return -y;
  }
  else {
    float y;
    if (rx < 0) y = radians(rx - 360);
    else y = radians(rx + 360);
    return sin(y);
  }
}

Test it with:

float[] tests = {0, 180, 360, -360,
                 Float.POSITIVE_INFINITY,
                 Float.NEGATIVE_INFINITY,
                 Float.NaN
                };
void setup(){
  for (float test : tests) {
    try {
      float good = sind(test);
      float bad = sin(radians(test));
      println(good, bad);
    } catch (IllegalArgumentException e) {
      println(e);
    }
  }
}

Output:

0.0 0.0
0.0 -8.742278E-8
0.0 1.7484555E-7
-0.0 -1.7484555E-7
java.lang.IllegalArgumentException: x cannot be infinite
java.lang.IllegalArgumentException: x cannot be infinite
java.lang.IllegalArgumentException: x cannot be NaN

So that looks great.

I believe that if getMeanAngle() uses sind() and cosd() instead of sin() and cos() it will give better answers. But Iā€™ll finish that up a little later.

1 Like

Brilliant! I am looking forward to seeing the finished thing.

That is what I was almost doing (but in a different way) with

  // find the angle from the arc tangent
  if (meansin < 0 && meancos > 0) {
    ha = 360 + degrees(atan(meansin/meancos));
  } else if (meancos < 0) {
    ha = 180 + degrees(atan(meansin/meancos));
  } else if (meansin > 0 && meancos > 0) {
    ha = degrees(atan(meansin/meancos));
  }

But perhaps I was using arctangent instead of the two-argument variant of arctangent. I wonder whether catching exceptions like this works with atan or atan2ā€¦ or neither!

1 Like

For the sake of completeness, I am sharing my three values version. I think it works. If so, one could use arrays instead of a variable for each value, a function, and a different way to catch exceptions. But see above: https://discourse.processing.org/t/getting-the-average-hsb-from-pixels/31965/27

int a1 = 240;
int a2 = 30;
int a3 = 200;

float sin1;
float sin2;
float sin3;

float cos1;
float cos2;
float cos3;

float meansin;
float meancos;

float ha;
float ho;



void setup() {

  // mean of sine of each angle
  sin1 = sin(radians(a1));
  sin2 = sin(radians(a2));
  sin3 = sin(radians(a3));

  meansin = (sin1 + sin2 + sin3)/3.0;

  // mean of cosine of each angle
  cos1 = cos(radians(a1));
  cos2 = cos(radians(a2));
  cos3 = cos(radians(a3));

  meancos = (cos1 + cos2 + cos3)/3.0;

  // find the angle from the arc tangent (unfiltered)
  ho = degrees(atan2(meansin, meancos));
  println("raw "+ho, ho+180, ho+360);

  // find the angle from the arc tangent
  if (meansin < 0 && meancos > 0) {
    ha = 360 + degrees(atan2(meansin, meancos));
    println(ha % 360, "mean sine less than zero, mean cosine greater than zero");
  } else if (meancos < 0) {
    ha = 180 + degrees(atan2(meansin, meancos));
    println(ha, "mean cosine less than zero");
  } else if (meansin > 0 && meancos > 0) {
    ha = degrees(atan2(meansin, meancos));
    println(ha, "mean sine greater than zero, mean cosine greater than zero");
  }

  //ha = (ha < 0) ? ha += 360 : 0;
  //println(ha);
}
1 Like

Inspecting the cases on meanAngle, Iā€™m now a bit suspicious about its utility for average color ā€“ at least in edge cases.

The essential issue is what is the desired ā€œcorrectā€ behavior when averaging two angles that are 180 degrees from each other. The incorrect implementations often give -90 or 90 (usually purple, sometimes green) any time they encounter pairs of angles at 180 degrees. The more ā€œcorrectā€ implementations always give 0 ā€“ in our domain, that would mean that two complimentary hues always produce red. Green + Purple? Red. Orange and Blue? Red.

colorMode(HSB, 360, 100, 100);
noStroke();
fill(45, 100, 100);
rect(0,0,width/2,height);
fill(225, 100, 100);
rect(width/2,0,width,height);

Screen Shot 2021-08-30 at 9.13.00 AM

What should the correct mean hue be for an image that is half orange and half blue? Should it be NaN? Should it always be 0? or should it be at perpendicular to the source angles ā€“ and if so, on which side of the circle? 135 (green) or 315 (magenta)?

So depending on how we implement mean angle, here are four potentially correct geometric answer colors, with black serving as NaN:

Screen Shot 2021-08-30 at 9.19.33 AM

1 Like

Well, I can offer three different answers.

(1) Choose the average that other tools give. GIMP gives magenta. My picker tool made in Processing gives 135 Green. But more people use GIMP so one could go for consistency with the most popular tools.

(2) The circular mean is useful for calendar dates, times, angles, hue, and many other cyclic quantities. One could choose an answer that made consistent sense, generally, across these different domains. Or, at least, is not nonsense for dates, times, angles.

(3) But probably the best answer is black, because no answer is better than a wrong answer. Green is obviously wrong, Magenta is obviously wrong. No colour is saying there is no correct answer.

1 Like

All great answers.

What Iā€™m wondering next is ā€“ supposing we have an image with three pixels ā€“ orange, blue, pink. Right now what Iā€™m seeing in meanAngle tests (modified to return 0 or NaN) is that the orange and blue ā€œmutually annihilateā€ ā€“ whatever the compliments are, they reduce to nothing, so the mean hue is whatever is left (e.g. the mean hue of orange-blue-pink is pink, and the mean hue of green-magenta-pink isā€¦ pink).

But we (probably?) donā€™t want an image with 999 orange pixels, 999 blue pixels, and two pink pixels to have an average hue of pinkā€¦

Perhaps the GIMP example is best. Make the code work for most cases. Sometimes it will throw up strange results like magenta as the average of orange and blue. The user looks at it goes yuk, then looks at what they were trying to average and knows from real world experience that you canā€™t mix orange and blue. When there are many different coloured pixels in a sample, the odd results will hopefully just look muddy brown.

An update! The code (below) gives the average without using sines, cosines, or two-argument arctangents. The solution is quite simple.
When any hue is over 180, subtract 360.

    if (hue(c) > 180) {
      h = h - 360;
    }

When the result is less than zero, add 360.

  if (h < 0) {
    h += 360;
  }

I think this solution is somewhere in the code from @quark and @glv but I had not looked closely because @quark thought his solution would only work consistently for two values. I have tested my code with three or more values in any order and cannot find a problem.
As for the orange and blue problem posed by @jeremydouglass, my code gives magenta. This is consistent with GIMP.

First the simple code to do the averaging. And then the same code, but with the extra functionality of moving the mouse to choose the sample area.

PImage img1; // the image to be investigated

void setup() {
  colorMode(HSB, 360, 100, 100);
  size(690, 660);
  noStroke();

//451
  img1 = loadImage("01.jpg");
} 

void draw() {
  background(0, 0, 100);

  // display the image under investigation
  image(img1, 0, 0); 


  // after averaged by function
  fill(sampleFromImage(img1)); 
  rect(0, img1.height+20, width, 189);
}




color sampleFromImage(PImage img1) {

  img1.loadPixels();
  int h = 0, s = 0, b = 0;
  for (int i=0; i<img1.pixels.length; i++) {
    color c = img1.pixels[i];

    if (hue(c) > 180) {
      h = h - 360;
    }



    h += hue(c);
    s += saturation(c);
    b += brightness(c);
  }

  h /= img1.pixels.length;
  s /= img1.pixels.length;
  b /= img1.pixels.length;


  if (h < 0) {
    h += 360;
  }

  fill(0);
  text("HSB "+h+" "+s+" "+b, 20, img1.height+14);

  return color(h, s, b);
}

With extra code for choosing the sample area:

// ( )
PImage img1; // the image to be investigated
PImage img2; // the detail of the image to be averaged

// ( )
int samplesize = 120; // an even number
int half = samplesize/2;


void setup() {
  colorMode(HSB, 360, 100, 100);
  size(1920, 1200);


  img1 = loadImage("01.jpg");
} 

void draw() {
  background(0, 0, 100);

  // (1) display the image under investigation
  image(img1, 0, 0); 

  // (2) the sample, but in a new image
  img2 = img1.get(mouseX-half, mouseY-half, samplesize, samplesize); 

  // (3) after averaged by function
  fill(sampleFromImage(img2)); 
  rect(img1.width-200, img1.height, 200, 200);

  // (4) show border of sample
  rectMode(CENTER);
  stroke(0);
  noFill();
  rect(mouseX, mouseY, samplesize, samplesize); 
  rectMode(CORNER); // revert to default
  noStroke(); // revert to default
}




color sampleFromImage(PImage img2) {

  img2.loadPixels();
  int h = 0, s = 0, b = 0;
  for (int i=0; i<img2.pixels.length; i++) {
    color c = img2.pixels[i];

    if (hue(c) > 180) {
      h = h - 360;
    }



    h += hue(c);
    s += saturation(c);
    b += brightness(c);
  }

  h /= img2.pixels.length;
  s /= img2.pixels.length;
  b /= img2.pixels.length;


  if (h < 0) {
    h += 360;
  }

  fill(0);
  text("HSB "+h+" "+s+" "+b, 20, img1.height+220);

  return color(h, s, b);
}
2 Likes