integral based circle anti-aliasing test

src/standalone_srgb_antialiasing_experiments.c

@@ -1,4 +1,5 @@
 #include <Windows.h>
+#include <math.h>
 
 ////////////////////////////////
 //- Experiment Configuration
@@ -47,6 +48,8 @@ typedef double F64;
 #define PtrOffsetOf(T,m) (void*)(&Member(T,m))
 #define OffsetOf(T,m) U64FromPtr(&Member(T,m))
 
+#define Clamp(a,x,b) (((x)<(a))?(a):((b)<(x))?(b):(x))
+
 #define allocate_array(T,c) (T*)(malloc(sizeof(T)*(c)))
 
 
@@ -91,6 +94,248 @@ static void opengl_draw_texture_geometry(Vertex *v, U64 count, U32 texture);
 static void opengl_draw_rect_geometry(Vertex *v, U64 count);
 
 
+
+////////////////////////////////
+//- Circle Generator
+
+static F32
+integral_indef(F32 r, F32 x){
+  Assert(fabs(x) <= fabs(r));
+  F32 t = asinf(x/r);
+  F32 result = r*r*0.5f*(t + 0.5f*sinf(2*t));
+  return(result);
+}
+
+static F32
+integral_def(F32 r, F32 min_x, F32 max_x){
+  F32 max_v = integral_indef(r, max_x);
+  F32 min_v = integral_indef(r, min_x);
+  F32 result = max_v - min_v;
+  return(result);
+}
+
+static B32
+approx_equal(F32 a, F32 b){
+  B32 result = (a - 0.0078125 <= b && b <= a + 0.0078125);
+  return(result);
+}
+
+static void
+sanity_test_integrals(void){
+  {
+    F32 a = integral_def(1.f, -1.f, 1.f);
+    Assert(approx_equal(a, 3.141593f*0.5f));
+  }
+  {
+    F32 a = integral_def(2.f, -2.f, 2.f);
+    Assert(approx_equal(a, 3.141593f*2.f));
+  }
+  {
+    F32 a = integral_def(2.f, 0.f, 2.f);
+    Assert(approx_equal(a, 3.141593f));
+  }
+  {
+    F32 a0 = integral_def(3.f, 0.f, 1.f);
+    F32 a1 = integral_def(3.f, 1.f, 2.f);
+    F32 a2 = integral_def(3.f, 2.f, 3.f);
+    Assert(approx_equal(a0 + a1 + a2, 3.141593f*9.f/4.f));
+  }
+}
+
+static F32
+intersection_area(F32 cir_r, F32 unit_box_x0, F32 unit_box_y0){
+  F32 rr = cir_r*cir_r;
+  
+  F32 rel_y0 = unit_box_y0;
+  F32 rel_y1 = unit_box_y0 + 1.f;
+  F32 rel_x0 = unit_box_x0;
+  F32 rel_x1 = unit_box_x0 + 1.f;
+  
+  F32 y0y0 = rel_y0*rel_y0;
+  F32 y1y1 = rel_y1*rel_y1;
+  
+  F32 x0x0 = rel_x0*rel_x0;
+  F32 x1x1 = rel_x1*rel_x1;
+  
+  // check each corner for hitting the circle
+  B32 corner_bl = ((y0y0 + x0x0) <= rr);
+  B32 corner_tl = ((y1y1 + x0x0) <= rr);
+  B32 corner_br = ((y0y0 + x1x1) <= rr);
+  B32 corner_tr = ((y1y1 + x1x1) <= rr);
+  
+  // calculate area
+  F32 a = 0.f;
+#define COMB(a,b,c,d) ((a)|((b)<<1)|((c)<<2)|((d)<<3))
+  switch (COMB(corner_bl, corner_tl, corner_br, corner_tr)){
+    // full miss
+    case 0x0: a = 0.f; break;
+    // full hit
+    case 0xF: a = 1.f; break;
+    
+    // 1 corner hit
+    case COMB(1,0,0,0): case COMB(0,1,0,0):
+    case COMB(0,0,1,0): case COMB(0,0,0,1):
+    {
+      F32 intersection_y = 0.f;
+      if (corner_bl || corner_br){
+        intersection_y = rel_y0;
+      }
+      else{
+        intersection_y = rel_y1;
+      }
+      
+      F32 intersection_x = sqrtf(rr - intersection_y*intersection_y);
+      if (intersection_x < rel_x0 || rel_x1 < intersection_x){
+        intersection_x = -intersection_x;
+      }
+      
+      F32 x_min = 0.f;
+      F32 x_max = 0.f;
+      if (corner_bl || corner_tl){
+        x_min = rel_x0;
+        x_max = intersection_x;
+      }
+      else{
+        x_min = intersection_x;
+        x_max = rel_x1;
+      }
+      
+      F32 h = fabs(intersection_y);
+      a = integral_def(cir_r, x_min, x_max) - h*(x_max - x_min);
+    }break;
+    
+    // 2 corner hits (horizontal edge)
+    case COMB(1,0,1,0): case COMB(0,1,0,1):
+    {
+      F32 h = 0.f;
+      if (corner_bl){
+        h = rel_y0;
+      }
+      else{
+        h = -rel_y1;
+      }
+      a = integral_def(cir_r, rel_x0, rel_x1) - h;
+    }break;
+    
+    // 2 corner hits (vertical edge)
+    case COMB(1,1,0,0): case COMB(0,0,1,1):
+    {
+      F32 w = 0.f;
+      if (corner_bl){
+        w = rel_x0;
+      }
+      else{
+        w = -rel_x1;
+      }
+      a = integral_def(cir_r, rel_y0, rel_y1) - w;
+    }break;
+    
+    // 3 corner hits
+    case COMB(0,1,1,1): case COMB(1,0,1,1):
+    case COMB(1,1,0,1): case COMB(1,1,1,0):
+    {
+      F32 intersection_y = 0.f;
+      if (!corner_bl || !corner_br){
+        intersection_y = rel_y0;
+      }
+      else{
+        intersection_y = rel_y1;
+      }
+      
+      F32 intersection_x = sqrtf(rr - intersection_y*intersection_y);
+      if (intersection_x < rel_x0 || rel_x1 < intersection_x){
+        intersection_x = -intersection_x;
+      }
+      
+      F32 w_extra = 0.f;
+      F32 x_min = 0.f;
+      F32 x_max = 0.f;
+      if (!corner_bl || !corner_tl){
+        x_min = rel_x0;
+        x_max = intersection_x;
+        w_extra = rel_x1 - intersection_x;
+      }
+      else{
+        x_min = intersection_x;
+        x_max = rel_x1;
+        w_extra = intersection_x - rel_x0;
+      }
+      
+      F32 h = fabs(intersection_y) - 1.f;
+      a = integral_def(cir_r, x_min, x_max) - h*(x_max - x_min) + w_extra;
+    }break;
+    
+    // these cases should be impossible
+    default:
+    {
+      a = -1.f;
+    }break;
+  }
+#undef COMB
+  
+  return(a);
+}
+
+static void
+sanity_test_circle_area(void){
+  {
+    F32 cr = 5.f;
+    F32 cx = 5.f;
+    F32 cy = 5.f;
+    F32 sum_area = 0.f;
+    for (U32 y = 0; y < 10; y += 1){
+      F32 box_y0 = (F32)y;
+      for (U32 x = 0; x < 10; x += 1){
+        F32 box_x0 = (F32)x;
+        F32 a = intersection_area(cr, box_x0 - cx, box_y0 - cy);
+        Assert(0.f <= a && a <= 1.f);
+        sum_area += a;
+      }
+    }
+    F32 formula_area = 3.141593f*cr*cr;
+    Assert(approx_equal(sum_area, formula_area));
+  }
+  
+  {
+    F32 cr = 8.f;
+    F32 cx = 8.1f;
+    F32 cy = 8.7f;
+    F32 sum_area = 0.f;
+    for (U32 y = 0; y < 20; y += 1){
+      F32 box_y0 = (F32)y;
+      for (U32 x = 0; x < 20; x += 1){
+        F32 box_x0 = (F32)x;
+        F32 a = intersection_area(cr, box_x0 - cx, box_y0 - cy);
+        Assert(0.f <= a && a <= 1.f);
+        sum_area += a;
+      }
+    }
+    F32 formula_area = 3.141593f*cr*cr;
+    Assert(approx_equal(sum_area, formula_area));
+  }
+}
+
+static void
+render_circle_in_buffer(F32 cir_r, F32 cir_x, F32 cir_y,
+                        U8 *buf, U32 buf_side_length){
+  for (U32 y = 0; y < buf_side_length; y += 1){
+    F32 box_y0 = (F32)y;
+    for (U32 x = 0; x < buf_side_length; x += 1){
+      F32 box_x0 = (F32)x;
+      
+      // calculate area of the intersection (circle & box)
+      F32 a = intersection_area(cir_r, box_x0 - cir_x, box_y0 - cir_y);
+      Assert(0.f <= a && a <= 1.f);
+      
+      U32 vraw = (U32)(256.f*a);
+      U8 v = (vraw <= 255)?((U8)vraw):(U8)255;
+      buf[x + buf_side_length*y] = v;
+    }
+  }
+}
+
+
+
 ////////////////////////////////
 //- GL Definitions
 
@@ -434,6 +679,24 @@ static GLuint canvas_framebuffer = 0;
 
 static GLuint test3_texture = 0;
 
+#define       CIRCLE5_RADIUS 5.f
+#define       CIRCLE5_X 6.f
+#define       CIRCLE5_Y 6.f
+#define       CIRCLE5_SIDE 12
+static GLuint circle5_texture  = 0;
+
+#define       CIRCLE25_RADIUS 25.f
+#define       CIRCLE25_X 30.f
+#define       CIRCLE25_Y 30.f
+#define       CIRCLE25_SIDE 60
+static GLuint circle25_texture  = 0;
+
+#define       CIRCLE125_RADIUS 125.f
+#define       CIRCLE125_X 150.f
+#define       CIRCLE125_Y 150.f
+#define       CIRCLE125_SIDE 300
+static GLuint circle125_texture  = 0;
+
 static void
 opengl_render_init(void){
   //- setup basic shader program
@@ -535,6 +798,81 @@ opengl_render_init(void){
     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
   }
   
+  //- circle5 texture
+  {
+    U8 buf[CIRCLE5_SIDE*CIRCLE5_SIDE] = {0};
+    render_circle_in_buffer(CIRCLE5_RADIUS, CIRCLE5_X, CIRCLE5_Y,
+                            buf, CIRCLE5_SIDE);
+    
+    U8 img[CIRCLE5_SIDE*CIRCLE5_SIDE*3] = {0};
+    for (U32 i = 0; i < CIRCLE5_SIDE*CIRCLE5_SIDE; i += 1){
+      U8 c = buf[i];
+      img[i*3 + 0] = c;
+      img[i*3 + 1] = c;
+      img[i*3 + 2] = c;
+    }
+    
+    glGenTextures(1, &circle5_texture);
+    glBindTexture(GL_TEXTURE_2D, circle5_texture);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, CIRCLE5_SIDE, CIRCLE5_SIDE, 0,
+                 GL_RGB, GL_UNSIGNED_BYTE, img);
+    
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
+  }
+  
+  //- circle25 texture
+  {
+    U8 buf[CIRCLE25_SIDE*CIRCLE25_SIDE] = {0};
+    render_circle_in_buffer(CIRCLE25_RADIUS, CIRCLE25_X, CIRCLE25_Y,
+                            buf, CIRCLE25_SIDE);
+    
+    U8 img[CIRCLE25_SIDE*CIRCLE25_SIDE*3] = {0};
+    for (U32 i = 0; i < CIRCLE25_SIDE*CIRCLE25_SIDE; i += 1){
+      U8 c = buf[i];
+      img[i*3 + 0] = c;
+      img[i*3 + 1] = c;
+      img[i*3 + 2] = c;
+    }
+    
+    glGenTextures(1, &circle25_texture);
+    glBindTexture(GL_TEXTURE_2D, circle25_texture);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, CIRCLE25_SIDE, CIRCLE25_SIDE, 0,
+                 GL_RGB, GL_UNSIGNED_BYTE, img);
+    
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
+  }
+  
+  //- circle125 texture
+  {
+    U8 buf[CIRCLE125_SIDE*CIRCLE125_SIDE] = {0};
+    render_circle_in_buffer(CIRCLE125_RADIUS, CIRCLE125_X, CIRCLE125_Y,
+                            buf, CIRCLE125_SIDE);
+    
+    U8 img[CIRCLE125_SIDE*CIRCLE125_SIDE*3] = {0};
+    for (U32 i = 0; i < CIRCLE125_SIDE*CIRCLE125_SIDE; i += 1){
+      U8 c = buf[i];
+      img[i*3 + 0] = c;
+      img[i*3 + 1] = c;
+      img[i*3 + 2] = c;
+    }
+    
+    glGenTextures(1, &circle125_texture);
+    glBindTexture(GL_TEXTURE_2D, circle125_texture);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, CIRCLE125_SIDE, CIRCLE125_SIDE, 0,
+                 GL_RGB, GL_UNSIGNED_BYTE, img);
+    
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
+  }
+  
   //- check for setup errors
   if (glGetError() != 0){
     error_message("error in opengl renderer initialization");
@@ -638,6 +976,10 @@ WinMain(HINSTANCE hInstance,
         HINSTANCE hPrevInstance,
         LPSTR lpCmdLine,
         int nShowCmd){
+  sanity_test_integrals();
+  sanity_test_circle_area();
+  
+  // do graphics tests
   wgl_helper_make_graphics_window(hInstance);
   
   U64 frame_counter = 0;
@@ -675,6 +1017,7 @@ WinMain(HINSTANCE hInstance,
     F32 test4_y = graphics_h - 160.f;
     F32 test5_y = graphics_h - 340.f;
     F32 test6_y = graphics_h - 450.f;
+    F32 test7_y = graphics_h - 480.f;
     
     // render to canvas frame buffer
     {
@@ -976,6 +1319,90 @@ WinMain(HINSTANCE hInstance,
             opengl_draw_rect_geometry(v, 6);
           }
         }
+        
+        // 7: circles rendered with analytical antialiasing via integration
+        {
+          F32 y = test7_y;
+          F32 x = grad_x0;
+          
+          // circle5
+          {
+            F32 y1 = y;
+            F32 y0 = y - CIRCLE5_SIDE;
+            F32 x0 = x;
+            F32 x1 = x + CIRCLE5_SIDE;
+            
+            x += CIRCLE5_SIDE + 10;
+            
+            Vertex v[6];
+            v[0].px  = x0;   v[0].py  = y0;
+            v[0].uvx = 0.f;  v[0].uvy = 0.f;
+            v[1].px  = x0;   v[1].py  = y1;
+            v[1].uvx = 0.f;  v[1].uvy = 1.f;
+            v[2].px  = x1;   v[2].py  = y0;
+            v[2].uvx = 1.f;  v[2].uvy = 0.f;
+            v[3].px  = x0;   v[3].py  = y1;
+            v[3].uvx = 0.f;  v[3].uvy = 1.f;
+            v[4].px  = x1;   v[4].py  = y0;
+            v[4].uvx = 1.f;  v[4].uvy = 0.f;
+            v[5].px  = x1;   v[5].py  = y1;
+            v[5].uvx = 1.f;  v[5].uvy = 1.f;
+            
+            opengl_draw_texture_geometry(v, 6, circle5_texture);
+          }
+          
+          // circle25
+          {
+            F32 y1 = y;
+            F32 y0 = y - CIRCLE25_SIDE;
+            F32 x0 = x;
+            F32 x1 = x + CIRCLE25_SIDE;
+            
+            x += CIRCLE25_SIDE + 10;
+            
+            Vertex v[6];
+            v[0].px  = x0;   v[0].py  = y0;
+            v[0].uvx = 0.f;  v[0].uvy = 0.f;
+            v[1].px  = x0;   v[1].py  = y1;
+            v[1].uvx = 0.f;  v[1].uvy = 1.f;
+            v[2].px  = x1;   v[2].py  = y0;
+            v[2].uvx = 1.f;  v[2].uvy = 0.f;
+            v[3].px  = x0;   v[3].py  = y1;
+            v[3].uvx = 0.f;  v[3].uvy = 1.f;
+            v[4].px  = x1;   v[4].py  = y0;
+            v[4].uvx = 1.f;  v[4].uvy = 0.f;
+            v[5].px  = x1;   v[5].py  = y1;
+            v[5].uvx = 1.f;  v[5].uvy = 1.f;
+            
+            opengl_draw_texture_geometry(v, 6, circle25_texture);
+          }
+          
+          // circle125
+          {
+            F32 y1 = y;
+            F32 y0 = y - CIRCLE125_SIDE;
+            F32 x0 = x;
+            F32 x1 = x + CIRCLE125_SIDE;
+            
+            x += CIRCLE125_SIDE + 10;
+            
+            Vertex v[6];
+            v[0].px  = x0;   v[0].py  = y0;
+            v[0].uvx = 0.f;  v[0].uvy = 0.f;
+            v[1].px  = x0;   v[1].py  = y1;
+            v[1].uvx = 0.f;  v[1].uvy = 1.f;
+            v[2].px  = x1;   v[2].py  = y0;
+            v[2].uvx = 1.f;  v[2].uvy = 0.f;
+            v[3].px  = x0;   v[3].py  = y1;
+            v[3].uvx = 0.f;  v[3].uvy = 1.f;
+            v[4].px  = x1;   v[4].py  = y0;
+            v[4].uvx = 1.f;  v[4].uvy = 0.f;
+            v[5].px  = x1;   v[5].py  = y1;
+            v[5].uvx = 1.f;  v[5].uvy = 1.f;
+            
+            opengl_draw_texture_geometry(v, 6, circle125_texture);
+          }
+        }
       }
       
     }