[freetype2] gsoc-anurag-2023 7d1e2c4b1 05/11: [dense] Add drawing functions to rasterizer

[Werner Lemberg]

branch: gsoc-anurag-2023
commit 7d1e2c4b111125afa3964435ba426e516d835efb
Author: Anurag Thakur <anurag105csec21@bpitindia.edu.in>
Commit: Anurag Thakur <anurag105csec21@bpitindia.edu.in>

    [dense] Add drawing functions to rasterizer
    
    * src/dense/ftdense.c: (dense_render_line, dense_render_quadratic,
    dense_render_cubic, dense_render_glyph, dense_raster_render, Lerp):
    New Functions
---
 src/dense/ftdense.c | 268 +++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 264 insertions(+), 4 deletions(-)

diff --git a/src/dense/ftdense.c b/src/dense/ftdense.c
index 80eab48a2..d2a4cde84 100644
--- a/src/dense/ftdense.c
+++ b/src/dense/ftdense.c
@@ -32,6 +32,15 @@ typedef struct dense_TRaster_
 
 } dense_TRaster, *dense_PRaster;
 
+/* Linear interpolation between P0 and P1 */
+static FT_Vector
+Lerp( float T, FT_Vector P0, FT_Vector P1 )
+{
+  FT_Vector p;
+  p.x = P0.x + T * ( P1.x - P0.x );
+  p.y = P0.y + T * ( P1.y - P0.y );
+  return p;
+}
 
 static int
 dense_move_to( const FT_Vector* to, dense_worker* worker )
@@ -54,11 +63,112 @@ dense_line_to( const FT_Vector* to, dense_worker* worker )
 }
 
 void
-dense_render_line( dense_worker* worker, TPos tox, TPos toy )
+dense_render_line( dense_worker* worker, FT_Pos tox, FT_Pos toy )
 {
-  return;
+  float from_x = worker->prev_x;
+  float from_y = worker->prev_y;
+  if ( from_y == toy )
+    return;
+
+
+  from_x /= 256.0;
+  from_y /= 256.0;
+  float to_x = tox / 256.0;
+  float to_y = toy / 256.0;
+
+
+  float dir;
+  if ( from_y < to_y )
+    dir = 1;
+  else
+  {
+    dir = -1;
+    FT_SWAP(from_x, to_x );
+    FT_SWAP(from_y, to_y );
+  }
+
+  // Clip to the height.
+  if ( from_y >= worker->m_h || to_y <= 0 )
+    return;
+
+  float dxdy = ( to_x - from_x ) / (float)( to_y - from_y );
+  if ( from_y < 0 )
+  {
+    from_x -= from_y * dxdy;
+    from_y = 0;
+  }
+  if ( to_y > worker->m_h )
+  {
+    to_x -= ( to_y - worker->m_h ) * dxdy;
+    to_y = (float)worker->m_h;
+  }
+
+  float  x       = from_x;
+  int    y0      = (int)from_y;
+  int    y_limit = (int)ceil( to_y );
+  float* m_a     = worker->m_a;
+
+  for ( int y = y0; y < y_limit; y++ )
+  {
+    int   linestart = y * worker->m_w;
+    float dy        = fmin( y + 1.0f, to_y ) - fmax( (float)y, from_y );
+    float xnext     = x + dxdy * dy;
+    float d         = dy * dir;
+
+    float x0, x1;
+    if ( x < xnext )
+    {
+      x0 = x;
+      x1 = xnext;
+    }
+    else
+    {
+      x0 = xnext;
+      x1 = x;
+    }
+
+    /*
+    It's possible for x0 to be negative on the last scanline because of
+    floating-point inaccuracy That would cause an out-of-bounds array access at
+    index -1.
+    */
+    float x0floor = x0 <= 0.0f ? 0.0f : (float)floor( x0 );
+
+    int   x0i    = (int)x0floor;
+    float x1ceil = (float)ceil( x1 );
+    int   x1i    = (int)x1ceil;
+    if ( x1i <= x0i + 1 )
+    {
+      float xmf = 0.5f * ( x + xnext ) - x0floor;
+      m_a[linestart + x0i] += d - d * xmf;
+      m_a[linestart + ( x0i + 1 )] += d * xmf;
+    }
+    else
+    {
+      float s   = 1.0f / ( x1 - x0 );
+      float x0f = x0 - x0floor;
+      float a0  = 0.5f * s * ( 1.0f - x0f ) * ( 1.0f - x0f );
+      float x1f = x1 - x1ceil + 1.0f;
+      float am  = 0.5f * s * x1f * x1f;
+      m_a[linestart + x0i] += d * a0;
+      if ( x1i == x0i + 2 )
+        m_a[linestart + ( x0i + 1 )] += d * ( 1.0f - a0 - am );
+      else
+      {
+        float a1 = s * ( 1.5f - x0f );
+        m_a[linestart + ( x0i + 1 )] += d * ( a1 - a0 );
+        for ( int xi = x0i + 2; xi < x1i - 1; xi++ )
+          m_a[linestart + xi] += d * s;
+        float a2 = a1 + ( x1i - x0i - 3 ) * s;
+        m_a[linestart + ( x1i - 1 )] += d * ( 1.0f - a2 - am );
+      }
+      m_a[linestart + x1i] += d * am;
+    }
+    x = xnext;
+  }
 }
 
+
 static int
 dense_conic_to( const FT_Vector* control,
                 const FT_Vector* to,
@@ -73,7 +183,56 @@ dense_render_quadratic( dense_worker*    worker,
                         FT_Vector* control,
                         FT_Vector* to )
 {
+  /*
+  Calculate devsq as the square of four times the
+  distance from the control point to the midpoint of the curve.
+  This is the place at which the curve is furthest from the
+  line joining the control points.
+
+  4 x point on curve = p0 + 2p1 + p2
+  4 x midpoint = 4p1
+
+  The division by four is omitted to save time.
+  */
+
+  FT_Vector aP0 = { DOWNSCALE( worker->prev_x ), DOWNSCALE( worker->prev_y ) };
+  FT_Vector aP1 = { control->x, control->y };
+  FT_Vector aP2 = { to->x, to->y };
+
+  float devx  = aP0.x - aP1.x - aP1.x + aP2.x;
+  float devy  = aP0.y - aP1.y - aP1.y + aP2.y;
+  float devsq = devx * devx + devy * devy;
+
+  if ( devsq < 0.333f )
+  {
+    dense_line_to( &aP2, worker );
     return;
+  }
+
+  /*
+  According to Raph Levien, the reason for the subdivision by n (instead of
+  recursive division by the Casteljau system) is that "I expect the flatness
+  computation to be semi-expensive (it's done once rather than on each 
potential
+  subdivision) and also because you'll often get fewer subdivisions. Taking a
+  circular arc as a simplifying assumption, where I get n, a recursive approach
+  would get 2^ceil(lg n), which, if I haven't made any horrible mistakes, is
+  expected to be 33% more in the limit".
+  */
+
+  const float tol = 3.0f;
+  int         n   = (int)floor( sqrt( sqrt( tol * devsq ) ) )/8;
+  FT_Vector p      = aP0;
+  float     nrecip = 1.0f / ( n + 1.0f );
+  float     t      = 0.0f;
+  for ( int i = 0; i < n; i++ )
+  {
+    t += nrecip;
+    FT_Vector next = Lerp( t, Lerp( t, aP0, aP1 ), Lerp( t, aP1, aP2 ) );
+    dense_line_to(&next, worker );
+    p              = next;
+  }
+
+  dense_line_to( &aP2, worker );
 }
 
 static int
@@ -92,7 +251,43 @@ dense_render_cubic( dense_worker* worker,
                     FT_Vector*    control_2,
                     FT_Vector*    to )
 {
-  return;
+  FT_Vector aP0 = { DOWNSCALE( worker->prev_x ), DOWNSCALE( worker->prev_y ) };
+  FT_Vector aP1 = { control_1->x, control_1->y };
+  FT_Vector aP2 = { control_2->x, control_2->y };
+  FT_Vector aP3 = { to->x, to->y };
+
+  float     devx   = aP0.x - aP1->x - aP1->x + aP2->x;
+  float     devy   = aP0.y - aP1->y - aP1->y + aP2->y;
+  float     devsq0 = devx * devx + devy * devy;
+  devx             = aP1->x - aP2->x - aP2->x + aP3->x;
+  devy             = aP1->y - aP2->y - aP2->y + aP3->y;
+  float devsq1     = devx * devx + devy * devy;
+  float devsq      = fmax( devsq0, devsq1 );
+
+  if ( devsq < 0.333f )
+  {
+    dense_render_line( worker, aP3->x, aP3->y );
+    return;
+  }
+
+  const float tol    = 3.0f;
+  int         n      = (int)floor( sqrt( sqrt( tol * devsq ) ) ) / 8;
+  FT_Vector   p      = aP0;
+  float       nrecip = 1.0f / ( n + 1.0f );
+  float       t      = 0.0f;
+  for ( int i = 0; i < n; i++ )
+  {
+    t += nrecip;
+    FT_Vector a    = Lerp( t, Lerp( t, aP0, *aP1 ), Lerp( t, *aP1, *aP2 ) );
+    FT_Vector b    = Lerp( t, Lerp( t, *aP1, *aP2 ), Lerp( t, *aP2, *aP3 ) );
+    FT_Vector next = Lerp( t, a, b );
+    dense_render_line( worker, next.x, next.y );
+    worker->prev_x = next.x;
+    worker->prev_y = next.y;
+    p              = next;
+  }
+
+  dense_line_to( &aP3, worker );
 }
 
 static int
@@ -152,13 +347,78 @@ dense_render_glyph( dense_worker* worker, const 
FT_Bitmap* target )
 {
   FT_Error error = FT_Outline_Decompose( &( worker->outline ),
                                          &dense_decompose_funcs, worker );
+  // Render into bitmap
+  const float* source = worker->m_a;
+
+  unsigned char* dest     = target->buffer;
+  unsigned char* dest_end = target->buffer + worker->m_w * worker->m_h;
+  float          value    = 0.0f;
+  while ( dest < dest_end )
+  {
+    value += *source++;
+    if ( value > 0.0f )
+    {
+      int n = (int)( fabs( value ) * 255.0f + 0.5f );
+      if ( n > 255 )
+        n = 255;
+      *dest = (unsigned char)n;
+    }
+    else
+      *dest = 0;
+    dest++;
+  }
+
+  free(worker->m_a);
   return error;
 }
 
 static int
 dense_raster_render( FT_Raster raster, const FT_Raster_Params* params )
 {
-  return 0;
+  const FT_Outline* outline    = (const FT_Outline*)params->source;
+  FT_Bitmap*  target_map = params->target;
+
+  dense_worker worker[1];
+
+  if ( !raster )
+    return FT_THROW( Invalid_Argument );
+
+  if ( !outline )
+    return FT_THROW( Invalid_Outline );
+
+  worker->outline = *outline;
+
+  if ( !target_map )
+    return FT_THROW( Invalid_Argument );
+
+  /* nothing to do */
+  if ( !target_map->width || !target_map->rows )
+    return 0;
+
+  if ( !target_map->buffer )
+    return FT_THROW( Invalid_Argument );
+
+  worker->m_origin_x = 0;
+  worker->m_origin_y = 0;
+  worker->m_w = target_map->pitch;
+  worker->m_h = target_map->rows;
+
+  int size = worker->m_w * worker->m_h + 4;
+
+  worker->m_a      = malloc( sizeof( float ) * size );
+  worker->m_a_size = size;
+
+  memset( worker->m_a, 0, ( sizeof( float ) * size ) );
+  /* exit if nothing to do */
+  if ( worker->m_w <= worker->m_origin_x || worker->m_h <= worker->m_origin_y )
+  {
+    return 0;
+  }
+
+  // Invert the pitch to account for different +ve y-axis direction in dense 
array
+  // (maybe temporary solution)
+  target_map->pitch *= -1;
+  return dense_render_glyph( worker, target_map );
 }
 
 FT_DEFINE_RASTER_FUNCS(