OxiGeo

Implementation code

// crates/oxigeo-wasm/src/warp.rs:312-377 — verbatim, the code running above
for out_row in 0..height {
    let world_y = f64::from(out_row).mul_add(-pixel_y, bbox[3] - pixel_y * 0.5);
    row.clear();
    for out_col in 0..width {
        let world_x = f64::from(out_col).mul_add(pixel_x, bbox[0] + pixel_x * 0.5);
        row.push(Coordinate::new(world_x, world_y));
    }

    projected.clear();
    // Point by point, never transform_batch: the batch path silently drops
    // +lat_0 on transverse-Mercator pairs, and one undefined point must not
    // cost the whole row either. See this module's header.
    //
    // (No backticks in this comment on purpose: it is inside the range the
    // page quotes verbatim, and the page holds it in a JS template literal.)
    for coord in &row {
        projected.push(
            inverse
                .transform(coord)
                .unwrap_or_else(|_| Coordinate::new(f64::NAN, f64::NAN)),
        );
    }

    for (out_col, point) in projected.iter().enumerate() {
        if !point.x.is_finite() || !point.y.is_finite() {
            continue;
        }
        let (px, py) = raster.transform.inverse(point.x, point.y);
        if !(px >= 0.0 && py >= 0.0) {
            continue;
        }
        // Both are non-negative and finite here, so the truncation to a
        // pixel index is the intended floor and cannot go negative.
        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
        let (column, srow) = (px as usize, py as usize);
        if column >= src_width || srow >= src_height {
            continue;
        }
        let Some(first) = raster.sample(0, column, srow) else {
            continue;
        };
        if raster.is_nodata(first) {
            continue;
        }
        let (red, green, blue) = if bands >= 3 {
            let (Some(g), Some(b)) = (
                raster.sample(1, column, srow),
                raster.sample(2, column, srow),
            ) else {
                continue;
            };
            (
                raster.stretch(0, first),
                raster.stretch(1, g),
                raster.stretch(2, b),
            )
        } else {
            let gray = raster.stretch(0, first);
            (gray, gray, gray)
        };
        let offset = ((out_row as usize) * (width as usize) + out_col) * 4;
        if let Some(pixel) = rgba.get_mut(offset..offset + 4) {
            pixel.copy_from_slice(&[red, green, blue, 255]);
        }
    }
}

This is the code running above, right now.