Converting cartesian pixels to polar pixels (I've largely revamped this entire question, though the motivation remains the same.)

revised question

I want to convert a raster of cartesian pixels into polar pixels.  Is there a sensible algorithm for doing this?  For example, how do I compute the value of the shaded (polar) pixel in the image below, given the value of the three (cartesian) pixels that it overlaps?

 

original question

Is there a reasonable way to compute the area of the intersection of a square and an annular section, as shown in the orange section below?



The motivation: I have a raster of square pixels, and I'm converting it to "polar pixels" -- I want to find out the contribution of each cartesian pixel to each polar pixel.
I have implemented the cartesian-to-polar-conversion and have used different interpolation methods:

1) nearest neighbor

2) a subsampling approach, which averages 81 subpixel locations

3) bilinear interpolation

The 2nd row in the image below shows detail magnifications of the output for the three approaches:



Here is the GLSL shader code for approach 1 and 3:

precision mediump float;
varying vec2 tc; // texture coordinate of current output pixel
uniform sampler2D inputImg; // input image
const float PI = 3.141592653589793238462643383;
void main() {
  float phi = tc.x; // phi is x-direction in output
  float r = tc.y;   // radius is y-direction in output
  float xx = r * cos(2.0 * PI * phi);  // unit circle in range [-1.0, 1.0]
  float yy = r * sin(2.0 * PI * phi);
  float x = xx * 0.5 + 0.5; // input texture coordinate in range [0.0, 1.0];
  float y = yy * 0.5 + 0.5;
  gl_FragColor = texture2D(inputImg, vec2(x,y));
}


This is the GLSL shader code for approach 2:

precision mediump float;
varying vec2 tc; // texture coordinate of current output pixel
uniform sampler2D inputImg; // input image
uniform int outputWidth; // width of output image
uniform int outputHeight; // height of output image
const float PI = 3.141592653589793238462643383;
const int samplesPerSide = 4;
void main() {
  int samplesPerDirection = samplesPerSide * 2 + 1;
  // compute size of subsample step in texture coordinates
  float sampleStepX = 1.0 / float(outputWidth * samplesPerDirection); 
  float sampleStepY = 1.0 / float(outputHeight * samplesPerDirection);
  vec4 sampleSum = vec4(0); // init sum to zero
  for(int i = -samplesPerSide; i &lt;= samplesPerSide; i++) {
    for(int j = -samplesPerSide; j &lt;= samplesPerSide; j++) {
      float phi = tc.x + float(i) * sampleStepX; // phi is x-direction in output
      float r = tc.y + float(j) * sampleStepY;   // radius is y-direction in output
      float xx = r * cos(2.0 * PI * phi);  // unit circle in range [-1.0, 1.0]
      float yy = r * sin(2.0 * PI * phi);
      float x = xx * 0.5 + 0.5; // input texture coordinate in range [0.0, 1.0];
      float y = yy * 0.5 + 0.5;
      sampleSum += texture2D(inputImg, vec2(x,y));
    }
  }
  gl_FragColor += sampleSum / float(samplesPerDirection * samplesPerDirection);
}


Approach 2, which is approx. factor 81 slower than 1 and 3, is probably the one you are looking for but I also like the bilinear interpolation result.

With a WebGL capable browser you can try out these implementations here:

https://www.gsn-lib.org/index.html#projectName=CartesianToPolar&amp;graphName=CartesianToPolar