Getting the average HSB from pixels

Processing Coding Questions
27

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