/**
RASTER_FP

A floating-point anti-aliasing renderer.
Graham Asher, August 2016.

Most of this code is derived from Raph Levien's font-rs code in the Rust language,
which is licensed under the Apache License, version 2.0.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

#include "raster_fp.h"
#include "assert.h"
#include "malloc.h"
#include "string.h"
#include "math.h"
#include "stdlib.h"

/** Creates a point by linear interpolation at the fraction T of the distance between aP0 and aP1. */
static RasterFP_Point Lerp(float aT,RasterFP_Point aP0,RasterFP_Point aP1)
    {
    RasterFP_Point p;
    p.m_x = aP0.m_x + aT * (aP1.m_x - aP0.m_x);
    p.m_y = aP0.m_y + aT * (aP1.m_y - aP0.m_y);
    return p;
    }

/** Initialises a RasterFP structure: the constructor, in C++ terms. */
void RasterFP_Create(RasterFP* aRasterFP)
    {
    assert(aRasterFP);
    memset(aRasterFP,0,sizeof(RasterFP));
    }

/**
Starts the rasterizing process.
Call this, then make a series of calls to RasterFP_DrawLine, RasterFP_DrawQuadratic and RasterFP_DrawCubic to render a shape.
The parameters aOriginX, aOriginY, aWidth and aHeight define a clipping rectangle. All rendering is clipped to this rectangle,
and the final bitmap represents the rectangle.

This function preserves the previous accumulation buffer (m_a) if possible, or creates a new one if more memory is needed.
To free up the memory, call RasterFP_Destroy before calling RasterFP_StartRasterizing.
*/
void RasterFP_StartRasterizing(RasterFP* aRasterFP,int aOriginX,int aOriginY,int aWidth,int aHeight)
    {
    assert(aRasterFP);
    assert(aWidth >= 0 && aHeight >= 0);

    aRasterFP->m_origin_x = aOriginX;
    aRasterFP->m_origin_y = aOriginY;
    aRasterFP->m_w = aWidth;
    aRasterFP->m_h = aHeight;

    int size = aRasterFP->m_w * aRasterFP->m_h + 4;
    if (aRasterFP->m_a_size < size)
        {
        free(aRasterFP->m_a);
        aRasterFP->m_a = malloc(sizeof(float) * size);
        if (aRasterFP->m_a == NULL)
            {
            aRasterFP->m_a_size = 0;
            return;
            }
        aRasterFP->m_a_size = size;
        }
    memset(aRasterFP->m_a,0,sizeof(float) * size);
    }

/** Cleans up resources used by a RasterFP structure: the destructor, in C++ terms.*/
void RasterFP_Destroy(RasterFP* aRasterFP)
    {
    assert(aRasterFP);
    free(aRasterFP->m_a);
    memset(aRasterFP,0,sizeof(RasterFP));
    }

/** Draws a straight line from aP0 to aP1. */
void RasterFP_DrawLine(RasterFP* aRasterFP,RasterFP_Point aP0,RasterFP_Point aP1)
    {
    assert(aRasterFP);
    if (aP0.m_y == aP1.m_y)
        return;

    aP0.m_x -= aRasterFP->m_origin_x;
    aP0.m_y -= aRasterFP->m_origin_y;
    aP1.m_x -= aRasterFP->m_origin_x;
    aP1.m_y -= aRasterFP->m_origin_y;

    float dir;
    if (aP0.m_y < aP1.m_y)
        dir = 1;
    else
        {
        dir = -1;
        RasterFP_Point temp = aP0;
        aP0 = aP1;
        aP1 = temp;
        }

    // Clip to the height.
    if (aP0.m_y >= aRasterFP->m_h || aP1.m_y <= 0)
        return;

    float dxdy = (aP1.m_x - aP0.m_x) / (aP1.m_y - aP0.m_y);
    if (aP0.m_y < 0)
        {
        aP0.m_x += aP0.m_y * dxdy;
        aP0.m_y = 0;
        }
    if (aP1.m_y > aRasterFP->m_h)
        {
        aP1.m_x -= (aP1.m_y - aRasterFP->m_h) * dxdy;
        aP1.m_y = (float)aRasterFP->m_h;
        }

    // Clip off parts of lines to the right of the raster.
    if (aP0.m_x >= aRasterFP->m_w && aP1.m_x >= aRasterFP->m_w)
        return;
    if (aP0.m_x < aRasterFP->m_w && aP1.m_x > aRasterFP->m_w)
        {
        aP1.m_y -= (aP1.m_x - aRasterFP->m_w) / dxdy;
        aP1.m_x = (float)aRasterFP->m_w;
        }
    else if (aP0.m_x > aRasterFP->m_w && aP1.m_x < aRasterFP->m_w)
        {
        aP0.m_y += (aP0.m_x - aRasterFP->m_w) / dxdy;
        aP0.m_x = (float)aRasterFP->m_w;
        }

    // Handle parts of the line left of x = 0 by clipping or recursive calls.
    if (aP0.m_x < 0 && aP1.m_x < 0)
        {
        aP0.m_x = aP1.m_x = 0;
        dxdy = 0;
        }
    else if (aP0.m_x < 0 && aP1.m_x > 0)
        {
        RasterFP_Point p, q;
        p.m_x = q.m_x = 0;
        p.m_y = aP0.m_y;
        q.m_y = aP0.m_y - aP0.m_x / dxdy;
        if (dir > 0)
            RasterFP_DrawLine(aRasterFP,p,q);
        else
            RasterFP_DrawLine(aRasterFP,q,p);
        aP0 = q;
        }
    else if (aP0.m_x > 0 && aP1.m_x < 0)
        {
        RasterFP_Point p,q;
        p.m_x = q.m_x = 0;
        p.m_y = aP1.m_y;
        q.m_y = aP1.m_y - aP1.m_x / dxdy;
        if (dir > 0)
            RasterFP_DrawLine(aRasterFP,q,p);
        else
            RasterFP_DrawLine(aRasterFP,p,q);
        aP1 = q;
        }

    float x = aP0.m_x;
    int y0 = (int)aP0.m_y;
    int y_limit = (int)ceil(aP1.m_y);
    float* m_a = aRasterFP->m_a;
    for (size_t y = y0; y < y_limit; y++)
        {
        size_t linestart = y * aRasterFP->m_w;
        float dy = min(y + 1.0f,aP1.m_y) - max((float)y,aP0.m_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;
            }

        float x0floor = (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;
        }
    }

/** Draws a quadratic spline with the control points aP0, aP1, and aP2. */
void RasterFP_DrawQuadratic(RasterFP* aRasterFP,RasterFP_Point aP0,RasterFP_Point aP1,RasterFP_Point aP2)
    {
    assert(aRasterFP);

    /*
    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.
    */
    float devx = aP0.m_x - aP1.m_x - aP1.m_x + aP2.m_x;
    float devy = aP0.m_y - aP1.m_y - aP1.m_y + aP2.m_y;
    float devsq = devx * devx + devy * devy;

    if (devsq < 0.333f)
        {
        RasterFP_DrawLine(aRasterFP,aP0,aP2);
        return;
        }

    float tol = 3.0;
    int n = (int)floor(sqrt(sqrt(tol * devsq)));
    RasterFP_Point p = aP0;
    float nrecip = 1.0f / (n + 1.0f);
    float t = 0.0;
    for (int i = 0; i < n; i++)
        {
        t += nrecip;
        RasterFP_Point next = Lerp(t,Lerp(t,aP0,aP1),Lerp(t,aP1,aP2));
        RasterFP_DrawLine(aRasterFP,p,next);
        p = next;
        }

    RasterFP_DrawLine(aRasterFP,p,aP2);
    }

/** Draws a cubic spline with the control points aP0, aP1, aP2 and aP3. */
void RasterFP_DrawCubic(RasterFP* aRasterFP,RasterFP_Point aP0,RasterFP_Point aP1,RasterFP_Point aP2,RasterFP_Point aP3)
    {
    assert(aRasterFP);

    float devx = aP0.m_x - aP1.m_x - aP1.m_x + aP2.m_x;
    float devy = aP0.m_y - aP1.m_y - aP1.m_y + aP2.m_y;
    float devsq0 = devx * devx + devy * devy;
    devx = aP1.m_x - aP2.m_x - aP2.m_x + aP3.m_x;
    devy = aP1.m_y - aP2.m_y - aP2.m_y + aP3.m_y;
    float devsq1 = devx * devx + devy * devy;
    float devsq = max(devsq0,devsq1);

    if (devsq < 0.333f)
        {
        RasterFP_DrawLine(aRasterFP,aP0,aP3);
        return;
        }

    float tol = 3.0;
    int n = (int)floor(sqrt(sqrt(tol * devsq)));
    RasterFP_Point p = aP0;
    float nrecip = 1.0f / (n + 1.0f);
    float t = 0.0;
    for (int i = 0; i < n; i++)
        {
        t += nrecip;
        RasterFP_Point a = Lerp(t,Lerp(t,aP0,aP1),Lerp(t,aP1,aP2));
        RasterFP_Point b = Lerp(t,Lerp(t,aP1,aP2),Lerp(t,aP2,aP3));
        RasterFP_Point next = Lerp(t,a,b);
        RasterFP_DrawLine(aRasterFP,p,next);
        p = next;
        }

    RasterFP_DrawLine(aRasterFP,p,aP3);
    }

/** Gets the bitmap created by the drawing functions. aBitmap must point to a buffer containing at least aRasterFP->m_w * aRasterFP->m_h bytes. */
void RasterFP_GetBitmap(RasterFP* aRasterFP,unsigned char* aBitmap)
    {
    assert(aRasterFP);
    assert(aBitmap);

    const float* source = aRasterFP->m_a;
    unsigned char* dest = aBitmap;
    unsigned char* dest_end = aBitmap + aRasterFP->m_w * aRasterFP->m_h;
    float value = 0;
    while (dest < dest_end)
        {
        value += *source++;
        if (value > 0)
            {
            int n = (int)(value * 255 + 0.5);
            if (n > 255)
                n = 255;
            *dest = (unsigned char)n;
            }
        dest++;
        }
    }