Small-time scanline VM

[p0nce]

Rationale

There is a lot of potential conversions, adaptation with colors.

Make tiny programs and opcodes in a "scanline VM" (that will process pixels identically as a row)
in order to be able to find some optimizations on the bytecode.

The individual opcodes are on the order of:

• convert rgba8 to float
• greyscale / luma
• exponentiation / sRGB mapping
  This will allow to support any input format to any output format, with a little runtime cost but smaller code in the end I think.
  The problem is all the scanlines functions we already have in Gamut.

Preliminary design

VM opcode design:

• 4 input ports, size of a scanline
• 4 output ports, size of a scanline
• 8 temporaries registers, size of a scanline
• can we avoid relocation with input/output pointers and length being passed as arguments? => yes, similar to shaders or GPGPU, input and output pointer are arguments, also item count. "Registers" are whole scanline allocated next to each other.
• math opcodes

A VM program is always about one scanline, possibly having 2 inputs or 2 outputs.

Ability to navigate: CIE XYZ, various RGB, various HSL spaces, power-laws, but also do any kind of blending.
Very fast startup.

[p0nce]

Example: sRGB to oklab

for (int n = 0; n < N; ++n)
    luma[n] = oklabLuminance(r, g, b);


float oklabLuminance(ubyte r, ubyte g, ubyte b) 
{
    static float srgbToLinear(float c) 
    {
        return (c <= 0.04045f) 
                ? (c / 12.92f) 
                : powf((c + 0.055f) / 1.055f, 2.4f);
    }

    // Step 1: Convert sRGB to linear sRGB
    float R = srgbToLinear(r / 255.0f);
    float G = srgbToLinear(g / 255.0f);
    float B = srgbToLinear(b / 255.0f);

    // Step 2: Convert linear sRGB (D65) to CIE XYZ (D50)
    float X = 0.4124564f * R + 0.3575761f * G + 0.1804375f * B;
    float Y = 0.2126729f * R + 0.7151522f * G + 0.0721750f * B;
    float Z = 0.0193339f * R + 0.1191920f * G + 0.9503041f * B;

    // Step 3: Convert linear XYZ (D50) to LMS
    float L =  0.4002f * X + 0.7076f * Y - 0.0808f * Z;
    float M = -0.2263f * X + 1.1653f * Y + 0.0457f * Z;
    float S =  0.0000f * X + 0.0000f * Y + 0.9182f * Z;

    // Step 4: Nonlinear transform (cube root)
    float l_ = cbrtf(L);
    float m_ = cbrtf(M);
    float s_ = cbrtf(S);

    // Step 4: Compute OKLab L
    float fL = 0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_;
    float fa = 1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_;
    float fb = 0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_; 

    return fL; // Perceptual luminance in range ~[0, 1]
}

would translate in opcodes:

    // as simplification most opcodes could happen in 4 channels, float
    OP_UBYTE_TO_FLOAT
    OP_POWER_LAW(INPUT#0, TEMP#0, SRGB_TO_LINEAR_RGB) 
    OP_DIVIDE_BY_255 (TEMP#0, TEMP#0)
    OP_LOAD_MATRIX_3x3(#MATRIX#0, RGB_D65_TO_CIE_XYZ_D50)
    OP_LOAD_MATRIX_3x3(#MATRIX#1, CIE_XYZ_D50_TO_LMS_D65)
    OP_MULT_MATRIX_3x3(TEMP#0, TEMP#0, #MATRIX#0)
    OP_MULT_MATRIX_3x3(TEMP#0, TEMP#0, #MATRIX#1)
    OP_LOAD_MATRIX_3x3(#MATRIX#0, LMS_TO_OKLAB)
    OP_MULT_MATRIX_3x3(TEMP#0, TEMP#0, #MATRIX#0)
    OP_CBRT
    OP_PACK_FIRST_CHANNEL(OUTPUT#0, TEMP#0)

And then it's easy to implement opcode fusion whenever a loop is too slow, for that the bytecode is iterated until a fixed point is attained.

• Matrix ports have predefined constants
• Buffers have predefined constants (for inputs and outputs)

Not sure about fixed opcode size or not.

[p0nce]

Requirements

• Should help with gamut scanline conversions
• Should help with scanline conversions in canvasity and dplug:canvas
• Should help with scanline conversions in stb_image_resize2.

Typical function:

void stbir__decode_uint8_linear_scaled(CoderHelper H)(float * decodep,
                                                      int width_times_channels,
                                                      const(void)* inputp)
{
  @restrict float* decode = decodep;
  float * decode_end = cast(float*) decode + width_times_channels;
  const(char)* input = cast(const(char)*)inputp;

  const(char)* end_input_m16 = input + width_times_channels - 16;
  if ( width_times_channels >= 16 )
  {
    decode_end -= 16;
    for(;;)
    {
      stbir__simdi i; stbir__simdi8 o0,o1;
      stbir__simdf8 of0, of1;
      STBIR_NO_UNROLL(decode);
      stbir__simdi_load( i, input );
      stbir__simdi8_expand_u8_to_u32( o0, o1, i );
      stbir__simdi8_convert_i32_to_float( of0, o0 );
      stbir__simdi8_convert_i32_to_float( of1, o1 );
      stbir__simdf8_mult( of0, of0, STBIR_max_uint8_as_float_inverted8);
      stbir__simdf8_mult( of1, of1, STBIR_max_uint8_as_float_inverted8);
      stbir__decode_simdf8_flip!H( of0 );
      stbir__decode_simdf8_flip!H( of1 );
      stbir__simdf8_store( decode + 0, of0 );
      stbir__simdf8_store( decode + 8, of1 );
      decode += 16;
      input += 16;
      if ( decode <= decode_end )
        continue;
      if ( decode == ( decode_end + 16 ) )
        break;
      decode = decode_end; // backup and do last couple
      input = end_input_m16;
    }
    return;
  }

Fot this to work it needs to support:

• L aka R

• LA aka RA

• AL aka AR

• RGB

• RGBA

• BGRA

• ARGB

• ABGR
  Note: there is one opcode currently lol

• optimized sRGB direct and inverse (stbir__srgb_uchar_to_linear_float)

• operations it does can support every colorspace in CSS, though it's a user library thing

[p0nce]

Decision 1 : Should internal registers and input/ouput pointers be tagged with their current type?

A - type tag at runtime corresponding to "pixel types" in Gamut.
or B - no type tag

Pros:

• can report errors on wrong programs
• an opcode can change meaning with various tags, small cost probably
• more safety when debugging bytecode programs
• opcode fusion looks easier with that dynamic type?

Cons:

• processors don't work like that
• doesn't seem to be needed, there is no user creating bytecode programs.
• more annoying API.
• since type tags are domain specific, the VM becomes harder to reused for other purposes, which are out of scope.

Use case:

• encounter opcode SRGB-to-LINEAR. Quite different for RGBA8, RGB8, and RGBA8 with companded alpha (rare).
• Without type tags the opcode fusions will be more expensive.
• With type tag the same opcode can be used.

Type tag could also be: (Number of channel per items + base data type)

=> yes, sounds like a good idea

[p0nce]

Decision 2 What's in the type tag? Semantic for colors or not?

Two general possibilities:

• A. PixelType-like, contains number of channel, base data type, and a color order semantic (where is alpha), premultiplied => N given to the VM is the number of pixels.

• B. Number of channel per items + base data type + location of alpha if any + premultiplied if any + (colorspace?) => N given to the VM is the number of (channel per items x base data type) constituents.

What to include
=> B since it can support many datatypes. Can do trivial audio work. Can do CMYK.

• Number of channel per items
• Base data type: ubyte (0..255), ushort(0..65535), float (0..1)
• Location of alpha, or none is alpha. Alpha not affected by power-curves.
• Premultiplied alpha bit.

[p0nce]

Mmmm it gets complicated.

Opcodes being polymorphic is a problem:

• can't know the types of output before executing
• need the types of input in order to optimize, compute memory size, etc.

The next design decision should be that opcodes aren't polymorphic and assume some types. unsure about this

EDIT on the other hand, need some amount of type-checking to reserve buffers close to each other, so if types and item count would travel the program it would allow some very nice optimizations.
Need a new program format that has: output type for each result. Output size for each result. scanline length (or maximum length?), and scanline type, travel the instructions. An optimized, attributed program could be produced like another IR.

or alternatively:

Each opcode know its output type, of which there isn't too many. Even if the inputs could be polymorphic themselves...

=> we will prolly need type-checking of some sort (can be a field in Instruction stream), do we allow a varying item count for each slice in the machine? Would be quite flexible. Item count would be embed in bytecode/type.

Provided we do type-checking, is item count static or also type-checked?

Pros

• Would support more slice operations, which makes the VM more general.
• Perhaps could trim the number of types, by removing the "channel" concept?
• Doesn't seem any slower.

Cons

• Harder to statically allocate registers before execution
• Conceptually more complicated
• Type structure will be 64-bit instead of 32, and Instruction will hold a Type.
• Removing the channel concept could be useful, but nothing is certain.

[p0nce]

Decision: remove typechecking, opcodes do not store a type. opcodes are fully monomorphic.

=> OK, but now instructions have "opcode modifiers" that will pass a higher slice lengh in order to share some loop generation (eg: all "float add" opcodes are now one, SVM_OP_FADD and can do 1 to 16 channels

[p0nce]

Need an assembler already.

This characterization of data types was nice to be a superset of what gamut and stb_image_resize2 require, but it turns out the ScanlineVM itself doesn't need to know image types.

So this is dead code, could be used to extend PixelType.

// FUTURE: for now this cannot represent BGR, only RGB order.
//         maybe encoder order in 2 bits: swap RB, 
//                                        and alpha-first vs alpha-last

/**
    Tells everything there is to know about a register or I/O in SVM.
    Maximum number of channels is 128.
    
    Note: for now, no info about RGB order.
     _______________________________________________________________
    | Bits          | Name     | Function                           |
    +---------------+----------+------------------------------------+
    |    0 to 2     | dataType | Data type of sample                |
    |               |          | eg: SVM_FLOAT                      |
    +---------------+----------+------------------------------------+
    |    3          | RESERVED | (Keep it at zero for now)          |
    +---------------+----------+------------------------------------+
    |    4          | hasAlpha | 1 if there is an alpha channel     |
    |               |          | eg: SVM_HAS_ALPHA                  |
    +---------------+----------+------------------------------------+
    |    5          | premul   | 1 if alpha is premultiplied        |
    |               |          | eg: SVM_PREMUL                     |
    +---------------+----------+------------------------------------+
    |    6 to 7     | RESERVED | (Keep it at zero for now)          |
    +---------------+----------+------------------------------------+
    |    8 to 14    | chans    | Samples in an item - 1 (0 to 127)  |
    |               |          | eg: SVM_1_CHANNEL  is  0           |
    |               |          |     SVM_2_CHANNELS     1           |
    |               |          |     SVM_3_CHANNELS     2           |
    |               |          |     SVM_4_CHANNELS     3           |
    +---------------+----------+------------------------------------+
    |   15          | RESERVED | (Keep it at zero for now)          |
    +---------------+----------+------------------------------------+
    |   16 to 23    | alphaLoc | Alpha location (0 to 255)          |
    +---------------+----------+------------------------------------+
    |   24 to 31    | RESERVED | (Keep it at zero for now)          |
    |_______________|__________|____________________________________|

 


    Now the **length** of each rows in terms of is given at execute
    time by the VM.
*/
alias SVMType = int;

/**
    Registers and I/O contains Items that consist of a row of 
    elements with a particular `SVMDataType`, arranged in a number
    of channels.

    Note that floating-point range is implicitely 0 to 1, and a 
    conversion would imply a scaling.
*/
alias SVMDataType = int;


/**
    Return data type for given `SVMType`.
*/
SVMDataType svmTypeDataType(SVMType type)
{
    return cast(SVMDataType)( type & 3 );
}


/**
    Return size of data type.
*/
int svmDataTypeSize(SVMDataType type)
{  
    switch (type)
    {
        case SVM_FLOAT: return 4;
        case SVM_BYTE: return 1;
        case SVM_SHORT: return 2;
        case SVM_INT: return 4;
        default:
            assert(false);
    }
}

/**
    Return size of one item.
*/
int svmTypeSize(SVMDataType type)
{  
    return svmDataTypeSize(svmTypeDataType(type)) * svmTypeChannels(type);
}

/**
    Return number of channels for given `SVMType`.
*/
int svmTypeChannels(SVMType type)
{
    return ((type & 0x7f00) >> 8) + 1;
}

/**
    Return channel index that is alpha, if any.
    Result is only valid if `SVM_ALPHA` bit set.

    See_also: `svmTypeHasAlpha`.
*/
int svmTypeAlphaLocation(SVMType type)
{
    return (type & 0x7f00) >> 8;
}

/**
    Does this type has alpha?
*/
bool svmTypeHasAlpha(SVMType type)
{
    return (type & SVM_ALPHA) != 0;
}

/**
    Does this type has alpha premul?
*/
bool svmTypeHasPremul(SVMType type)
{
    return (type & SVM_PREMUL) != 0;
}

///ditto
enum : SVMDataType
{
    SVM_FLOAT            = 0x0000,
    SVM_BYTE             = 0x0001,
    SVM_SHORT            = 0x0002,
    SVM_INT              = 0x0003,

    // Note: No unsigned, instead unsigned opcodes.
}

/**
    This data has one alpha channel.
*/
enum int SVM_ALPHA = 0x0010;

/**
    This data alpha, if it exists, is premultiplied.
*/
enum int SVM_PREMUL = 0x0020;

///
enum int SVM_1_CHANNEL  = 0x0000;
///
enum int SVM_2_CHANNELS = 0x0100;
///
enum int SVM_3_CHANNELS = 0x0200;
///
enum int SVM_4_CHANNELS = 0x0300;
///
enum int SVM_ALPHA_IN_CHAN0 = 0x0_0000;
///
enum int SVM_ALPHA_IN_CHAN1 = 0x1_0000;
///
enum int SVM_ALPHA_IN_CHAN2 = 0x2_0000;
///
enum int SVM_ALPHA_IN_CHAN3 = 0x3_0000;


/**
   Equivalents to Gamut's pixel types, plus some more.
*/

enum int SVM_L8       = SVM_BYTE  | SVM_1_CHANNEL;
enum int SVM_L16      = SVM_SHORT | SVM_1_CHANNEL;
enum int SVM_LF32     = SVM_FLOAT | SVM_1_CHANNEL;
enum int SVM_LA8      = SVM_BYTE  | SVM_2_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN1;
enum int SVM_LA16     = SVM_SHORT | SVM_2_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN1;
enum int SVM_LAF32    = SVM_FLOAT | SVM_2_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN1;
enum int SVM_LAP8     = SVM_BYTE  | SVM_2_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN1 | SVM_PREMUL;
enum int SVM_LAP16    = SVM_SHORT | SVM_2_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN1 | SVM_PREMUL;
enum int SVM_LAPF32   = SVM_FLOAT | SVM_2_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN1 | SVM_PREMUL;
enum int SVM_RGB8     = SVM_BYTE  | SVM_3_CHANNELS;
enum int SVM_RGB16    = SVM_SHORT | SVM_3_CHANNELS;
enum int SVM_RGBF32   = SVM_FLOAT | SVM_3_CHANNELS;
enum int SVM_RGBA8    = SVM_BYTE  | SVM_3_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN3;
enum int SVM_RGBA16   = SVM_SHORT | SVM_3_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN3;
enum int SVM_RGBAF32  = SVM_FLOAT | SVM_3_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN3;
enum int SVM_RGBAP8   = SVM_BYTE  | SVM_3_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN3 | SVM_PREMUL;
enum int SVM_RGBAP16  = SVM_SHORT | SVM_3_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN3 | SVM_PREMUL;
enum int SVM_RGBAPF32 = SVM_FLOAT | SVM_3_CHANNELS | SVM_ALPHA | SVM_ALPHA_IN_CHAN3 | SVM_PREMUL;

[p0nce]

State of the VM prototype as of August 2025:

• 8 general registers, 4 in, 4 out
• unpublished
• no optimizations yet
• working assembler
• 17 opcodes, really easy to add one
• can express the "sRGB u8 to oklab f32" challenge above.

Sample code:

void convertScanline_RGB8_to_oklabf32_SVM(
                                          ubyte[3]* input,
                                          float[3]* output,
                                          int count)
{
    ScanlineVM vm;
    Vec!SVMInstruction program;
    string error;
    svmAssemble(`
                u8_to_f32 x3 r0, src0;
                srgb_to_lin x3 r0, r0;
                mul_mat3f_v3f r0, src1, r0; 
                mul_mat3f_v3f r0, src2, r0;
                cbrtf x3 r0, r0;
                mul_mat3f_v3f dst0, src3, r0;
                `, program, error);
    if (error) writeln(error);
    assert(!error);

    float[9] linear_sRGB_to_CIE_XYZ_D50 =
    [
        0.4124564f, 0.3575761f, 0.1804375f,
        0.2126729f, 0.7151522f, 0.0721750f,
        0.0193339f, 0.1191920f, 0.9503041f
    ];

    float[9] CIE_XYZ_D50_to_LMS =
    [
         0.4002f, 0.7076f, -0.0808f,
        -0.2263f, 1.1653f,  0.0457f,
         0.0000f, 0.0000f,  0.9182f
    ];

    float[9] LMS_TO_OKLAB =
    [
        0.2104542553f,  0.7936177850f, -0.0040720468f,
        1.9779984951f, -2.4285922050f,  0.4505937099f,
        0.0259040371f,  0.7827717662f, -0.8086757660f
    ];    
    vm.inputArg(0, input[0..count]);
    vm.inputArg(1, linear_sRGB_to_CIE_XYZ_D50);
    vm.inputArg(2, CIE_XYZ_D50_to_LMS);
    vm.inputArg(3, LMS_TO_OKLAB);
    vm.outputArg(0, output[0..count]);
    vm.itemCount(count);
    vm.execute(program);
}

Some opcodes can operate not on scanline but individual values, in this case the start of a register is used.
No gain has been mesured yet. The end game is to compile and optimize the bytecode once, and execute on each scanline with only some arguments changing. The input arguments are guaranteed unwriteable for this reason.
Large opportunities for optimized opcodes, working on whole scanlines at once.

[p0nce]

Crazy thing is that now that it works, the ScanlineVM routine above is 20% faster without any bytecode optimization or SIMD opcode implementation, in both x86 and arm. niiice
Despite traversing scanlines 6 times instead of 1, (EDIT: probably slower that a full SIMD solution competely specialized, but not sure. We don't want to optimize all these loops anyway)

[p0nce]

Implementing segmenting / aka batch size ! Should reduce memory requirements by a lot, and perhaps improve batch size, much like buffer max size in plugins.

[p0nce]

Conversion oklab back and forth for ~1Mpixel, best of 3
aarch64 Windows LDC 1.40.1

Test 1 GP register only use their first items

Results

10.885 (baseline, no buffer split)
10.888 (batch 2048)
10.809 (batch 1024)
10.788 (batch 512)
10.538 (batch 256) ***WINNER***
10.730 (batch 128)
11.018 (batch 64)
11.452 (batch 32)
11.764 (batch 16)

[p0nce]

Now that everything works:

11.344 us (no buffer split) (11kb allocated)
10.640 us (batch 256)
10.798 us (batch 128) (1.4kb allocated)
10.864 us (batch 64)
10.999 us (batch 32) (384 bytes allocated)

[p0nce]

Current meta: CTFE assembler, working buffer (no size given tho), and it's pretty fast to "boot" as a "VM"
Default batch size is 128

SVM programs have no loops, run each items independently. There is no cross-lane thing, if it was there it would kinda break batching, so need to think about that.

[p0nce]

Here is a loop that doesn't work, example program:

/// Compute a sRGB image luminance.
/// source MUST be PixelType.rgba8
/// dest MUST be PixelType.l32f
void imageLuma(ref Image dest, 
               ref Image source)
{
    assert(dest.width  == source.width);
    assert(dest.height == source.height);
    assert(source.type == PixelType.rgba8);    
    assert(dest.type   == PixelType.lf32);

    int W = dest.width;
    int H = dest.height;

    version(Use_SVM)
    {
        import svm;

        ScanlineVM vm;
        vm.itemCount(W);
        vm.batchSize(512);

        static immutable float[3] LUMINANCE_WEIGHTS = 
        [
            0.2126f * 255, 0.7152f * 255, 0.0722f * 255
        ];
        vm.inputArg(1, LUMINANCE_WEIGHTS.ptr);

        for (int y = 0; y < H; ++y)
        {
            RGBA8* inp  = cast(RGBA8*) source.scanptr(y);
            float* outp = cast(float*) dest.scanptr(y);

            vm.inputArg(0, inp);
            vm.outputArg(0, outp);            

            enum program = svmAsm(`
                drop_alpha_u8 r0, src0;
                u8_to_f32 x3 r1, r0;
                chan_mix_v3f dst0, src1, r1;`);
            vm.execute(program);
        }
    }
    else
    {
        for (int y = 0; y < H; ++y)
        {
            RGBA8* inp  = cast(RGBA8*) source.scanptr(y);
            float* outp = cast(float*) dest.scanptr(y);

            for (int x = 0; x < W; ++x)
            {
                RGBA8 px = inp[x];
                float r = px.r;
                float g = px.g;
                float b = px.b;
                float lum = 0.2126f * r + 0.7152f * g + 0.0722f;
                outp[x] = lum;
            }
        }
    }
}

and performance look like this (whole image, best of 5 runs)

batch 1024 => 471 us
batch 512 => 475 us
batch 256 => 500 us
batch 128 => 576 us

original: 194 us

which makes sense this the opcodes aren't all optimized yet and also traversing the data thrice is naturally suboptimal in this case. I expect since the beginning a lot of new opcodes will be needed over time.

[p0nce]

All this is supposed to improve build times, often speed, and executable size.

For this one loop computing a 3x3 convolution:

bool useSVM = true;
    if (useSVM)
    {
        import svm;
        ScanlineVM vm;
        vm.itemCount(W);
        vm.batchSize(1024); 
        //char[10000] buf; slower!
        //vm.workBuf(buf);
        static immutable float[3] HZ_L0 = [  -0.1875 ,  0,  0.1875  ];
        static immutable float[3] HZ_L1 = [   -0.625 ,  0,  0.625  ];
        static immutable float[3] VT_L0 = [  -0.1875, -0.625, -0.1875 ];
        static immutable float[3] VT_L2 = [   0.1875,  0.625,  0.1875 ];
        vm.inputArg(3, HZ_L0);
        vm.inputArg(4, HZ_L1);
        vm.inputArg(6, VT_L0);
        vm.inputArg(7, VT_L2);

        for (int y = 0; y < H; ++y)
        {
            float[] row     = cast(float[]) source.scanline(y);
            float[] rowUp   = cast(float[]) source.scanline(max2(y - 1, 0));
            float[] rowDown = cast(float[]) source.scanline(min2(y + 1, H - 1));
            float[] outx    = cast(float[]) destX.scanline(y);
            float[] outy    = cast(float[]) destY.scanline(y);
            vm.inputArg(0, rowUp);
            vm.inputArg(1, row);
            vm.inputArg(2, rowDown);
            static immutable program 
                = svmAsm(`
                          convolve_v3f dst0, src0, src3;
                          convolve_v3f r1, src1, src4;
                          convolve_v3f r2, src2, src3;
                          accumf dst0, r1, r2;
                          convolve_v3f r1, src0, src6;
                          convolve_v3f r2, src2, src7;
                          addf dst1, r1, r2;
                         `);
            vm.outputArg(0, outx);
            vm.outputArg(1, outy);
            vm.execute(program);
        }
    }
    else
    {
        // A B C           -3 0  3           -3 -10 -3
        // D E F          -10 0 10 for horz   0   0  0
        // G H I           -3 0  3            3  10  3
        for (int y = 0; y < H; ++y)
        {
            float[] row     = cast(float[]) source.scanline(y);
            float[] rowUp   = cast(float[]) source.scanline(max2(y - 1, 0));
            float[] rowDown = cast(float[]) source.scanline(min2(y + 1, H - 1));
            float[] outx    = cast(float[]) destX.scanline(y);
            float[] outy    = cast(float[]) destY.scanline(y);

            for (int x = 0; x < W; ++x)
            {
                int xm = max2(x-1, 0);
                int xp = min2(x+1, W-1);
                float A =   rowUp[xm]; float B =   rowUp[x]; float C =   rowUp[xp];
                float D =     row[xm]; float E =     row[x]; float F =     row[xp];
                float G = rowDown[xm]; float J = rowDown[x]; float I = rowDown[xp];
                float dx = 0.1875f * (C + I - A - G) + 0.625f * (F - D);
                float dy = 0.1875f * (G + I - A - C) + 0.625f * (J - B);

                outx[x] += dx;
                outx[x] += dy;

            }
      }
 }

• Best timing with SVM is 2 ms 65 us vs 1 ms 850 us for normal code.