Fast scale-and-rotate blit.

scale_rotate.c

typedef unsigned char Pixel_t;

typedef struct _Point_t {
    int x, y;
} Point_t;

typedef struct _Rectangle_t {
    int x, y;
    int width, height;
} Rectangle_t;

typedef struct _Surface_t {
    Pixel_t *data;
    Pixel_t **data_rows;
    int width, height;
} Surface_t;

static inline int _iroundf(float x)
{
    return (int)floorf(x + 0.5f);
}

void scale_rotate(Surface_t *destination, Point_t position,
                  const Surface_t *source, Rectangle_t area,
                  float scale_x, float scale_y, float angle, float anchor_x, float anchor_y)
{
    const float w = (float)area.width;
    const float h = (float)area.height;
    const float sw = w * fabs(scale_x);
    const float sh = h * fabs(scale_y);

    const float sax = (w - 1.0f) * anchor_x; // Anchor points, relative to the source and destination areas.
    const float say = (h - 1.0f) * anchor_y;
    const float dax = (sw - 1.0f) * anchor_x;
    const float day = (sh - 1.0f) * anchor_y;

    const float sx = area.x;
    const float sy = area.y;
    const float dx = position.x;
    const float dy = position.y;

    const float c = cosf(angle);
    const float s = sinf(angle);

    // The counter-clockwise 2D rotation matrix is
    //
    //      |  c  -s |
    //  R = |        |
    //      |  s   c |
    //
    // In order to calculate the clockwise rotation matrix one can use the
    // similarities `cos(-a) = cos(a)` and `sin(-a) = -sin(a)` and get
    //
    //      |  c   s |
    //  R = |        |
    //      | -s   c |

    // Rotate the four corners of the scaled image to compute the rotate/scaled AABB.
    //
    // Note that we aren *not* adding `dst/dty` on purpose to rotate around the anchor point.
    const float aabb_x0 = -dax;
    const float aabb_y0 = -day;
    const float aabb_x1 = sw - 1.0f - dax;
    const float aabb_y1 = sh - 1.0f - day;

    const float x0 = c * aabb_x0 - s * aabb_y0;
    const float y0 = s * aabb_x0 + c * aabb_y0;

    const float x1 = c * aabb_x1 - s * aabb_y0;
    const float y1 = s * aabb_x1 + c * aabb_y0;

    const float x2 = c * aabb_x1 - s * aabb_y1;
    const float y2 = s * aabb_x1 + c * aabb_y1;

    const float x3 = c * aabb_x0 - s * aabb_y1;
    const float y3 = s * aabb_x0 + c * aabb_y1;

    int drawing_region_x0 = _iroundf(fmin(fmin(fmin(x0, x1), x2), x3) + dx);
    int drawing_region_y0 = _iroundf(fmin(fmin(fmin(y0, y1), y2), y3) + dy);
    int drawing_region_x1 = _iroundf(fmax(fmax(fmax(x0, x1), x2), x3) + dx);
    int drawing_region_y1 = _iroundf(fmax(fmax(fmax(y0, y1), y2), y3) + dy);

    if (drawing_region_x0 < clipping_region->x0) {
        drawing_region_x0 = clipping_region->x0;
    }
    if (drawing_region_y0 < clipping_region->y0) {
        drawing_region_y0 = clipping_region->y0;
    }
    if (drawing_region_x1 > clipping_region->x1) {
        drawing_region_x1 = clipping_region->x1;
    }
    if (drawing_region_y1 > clipping_region->y1) {
        drawing_region_y1 = clipping_region->y1;
    }

    const int width = drawing_region_x1 - drawing_region_x0 + 1;
    const int height = drawing_region_y1 - drawing_region_y0 + 1;
    if ((width <= 0) || (height <= 0)) { // Nothing to draw! Bail out!
        return;
    }

    const int sminx = area.x;
    const int sminy = area.y;
    const int smaxx = area.x + (int)area.width - 1;
    const int smaxy = area.y + (int)area.height - 1;

    const float M11 = c / scale_x;  // Since we are doing an *inverse* transformation, we combine rotation and *then* scaling (TRS -> SRT).
    const float M12 = s / scale_x;  // | 1/sx    0 | |  c s |
    const float M21 = -s / scale_y; // |           | |      |
    const float M22 = c / scale_y;  // |    0 1/sy | | -s c |

    const float tlx = drawing_region_x0 - dx; // Transform the top-left corner of the to-be-drawn rectangle to texture space.
    const float tly = drawing_region_y0 - dy; // (could differ from AABB x0 due to clipping, we need to compute it again)
    float ou = (tlx * M11 + tly * M12) + sax + sx; // Offset to the source texture quad.
    float ov = (tlx * M21 + tly * M22) + say + sy;

    if (scale_x < 0.0f) { // Fix flipping offset, relative to anchor. Presumably correct... :)
        const float factor = 0.5f - anchor_x;
        const float sign = (factor > 0.0f) - (factor < 0.0f);
        ou += factor * 2.0f * (float)area.width - sign / fabs(scale_x);
    }
    if (scale_y < 0.0f) {
        const float factor = 0.5f - anchor_y;
        const float sign = (factor > 0.0f) - (factor < 0.0f);
        ov += factor * 2.0f * (float)area.height - sign / fabs(scale_y);
    }

    for (int y = dminy; y <= dmaxy; ++y) {
        float u = ou;
        float v = ov;

        for (int x = dminx; x <= dmaxx; ++x) {
            int xx = _iroundf(u); // Round down, to preserve negative values as such (e.g. `-0.3` is `-1`)
            int xy = _iroundf(v);

            if (xx >= sminx && xx <= smaxx && xy >= sminy && xy <= smaxy) {
                const Pixel_t *src = source->data_rows[xy] + xx;
                Pixel_t *dst = destination->data_rows[y] + x;
                *dst = *src;
            }

            u += M11;
            v += M21;
        }

        ou += M12;
        ov += M22;
    }
}