Given a matrix,

myx myy myz myw

mzx mzy mzz mzw

mwx mwy mwz mww

A matrix vector multiply is defined as,

y' = x*myx + y*myy + z*myz + w*myw

z' = x*mzx + y*mzy + z*mzz + w*mzw

w' = x*mwx + y*mwy + z*mwz + w*mww

This doc also assumes one is using an infinite projection matrix defined as follows,

0 ? 0____ 0_____

0 0 (1-k) (k-1)n

0 0 1____ 0_____

Where

Note I'm using INT24 depth buffer here instead of inverted FP32 because INT24 still has better performance on some GPUs. While this is showing the process, one wouldn't want to do this in the shader as seen here because precision would be crap when view translation gets large (more on this later). Inputs,

d := depth fetched from INT24 depth buffer

n := constant near plane

k := constant insuring geometry does not hit far plane

First convert from texture XY to view XY. Note the DX convention is 1 at the top of the screen and -1 at the bottom (y is flipped compared to texture).

y' = 1 - 2*y

Next reverse the

Then grab world coordinates via a matrix multiply times the inverse of the view projection matrix,

Then back-project to the prior frame texture coordinates,

xy''''' = xy'''' * (1.0/w'''')

xy'''''' = {x''''' * 0.5 + 0.5, y''''' * (-0.5) + 0.5}

Then diff to generate the motion vector,

The example code below is public domain as it is just a simple equation in HLSL and C form which was generated by Mathematica. If this isn't enough, this generic code is covered in the unlicense,

Anyone is free to copy, modify, publish, use, compile, sell, or distribute this software, either in source code form or as a compiled binary, for any purpose, commercial or non-commercial, and by any means.

In jurisdictions that recognize copyright laws, the author or authors of this software dedicate any and all copyright interest in the software to the public domain. We make this dedication for the benefit of the public at large and to the detriment of our heirs and successors. We intend this dedication to be an overt act of relinquishment in perpetuity of all present and future rights to this software under copyright law.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Next grab your favorite symbolic solver like Mathematica and simply this. This step is required for good precision. The

float3 xyd, // {x,y} := {0,0}=topLeft, {1,1}=bottomRight screen, d = fetched depth

// constants generated by CameraMotionConstants() function

float4 c0, float4 c1, float4 c2, float4 c3, float4 c4) {

float2 mv;

float rcpScale = 1.0/(dot(xyd, c0.xyz) + c0.w);

mv.x = ((xyd.x * ((c1.x * xyd.y) + (c1.y * xyd.z) + c1.z)) + (c1.w * xyd.y) + (c2.x * xyd.x * xyd.x) + (c2.y * xyd.z) + c2.z) * rcpScale;

mv.y = ((xyd.y * ((c3.x * xyd.x) + (c3.y * xyd.z) + c3.z)) + (c3.w * xyd.x) + (c4.x * xyd.y * xyd.y) + (c4.y * xyd.z) + c4.z) * rcpScale;

return mv; }

And the code which can generate the constants,

float* restrict const cb, // constant buffer output

const double* restrict const c, // current view projection matrix inverse

const float* restrict const p // prior view projection matrix

) {

const double cxx = c[ 0]; const double cxy = c[ 1]; const double cxz = c[ 2]; const double cxw = c[ 3];

const double cyx = c[ 4]; const double cyy = c[ 5]; const double cyz = c[ 6]; const double cyw = c[ 7];

const double czx = c[ 8]; const double czy = c[ 9]; const double czz = c[10]; const double czw = c[11];

const double cwx = c[12]; const double cwy = c[13]; const double cwz = c[14]; const double cww = c[15];

const double pxx = (double)(p[ 0]); const double pxy = (double)(p[ 1]); const double pxz = (double)(p[ 2]); const double pxw = (double)(p[ 3]);

const double pyx = (double)(p[ 4]); const double pyy = (double)(p[ 5]); const double pyz = (double)(p[ 6]); const double pyw = (double)(p[ 7]);

const double pwx = (double)(p[12]); const double pwy = (double)(p[13]); const double pwz = (double)(p[14]); const double pww = (double)(p[15]);

// c0

cb[0] = (float)(4.0*(cwx*pww + cxx*pwx + cyx*pwy + czx*pwz));

cb[1] = (float)((-4.0)*(cwy*pww + cxy*pwx + cyy*pwy + czy*pwz));

cb[2] = (float)(2.0*(cwz*pww + cxz*pwx + cyz*pwy + czz*pwz));

cb[3] = (float)(2.0*(cww*pww - cwx*pww + cwy*pww + (cxw - cxx + cxy)*pwx + (cyw - cyx + cyy)*pwy + (czw - czx + czy)*pwz));

// c1

cb[4] = (float)(( 4.0)*(cwy*pww + cxy*pwx + cyy*pwy + czy*pwz));

cb[5] = (float)((-2.0)*(cwz*pww + cxz*pwx + cyz*pwy + czz*pwz));

cb[6] = (float)((-2.0)*(cww*pww + cwy*pww + cxw*pwx - 2.0*cxx*pwx + cxy*pwx + cyw*pwy - 2.0*cyx*pwy + cyy*pwy + czw*pwz - 2.0*czx*pwz + czy*pwz - cwx*(2.0*pww + pxw) - cxx*pxx - cyx*pxy - czx*pxz));

cb[7] = (float)(-2.0*(cyy*pwy + czy*pwz + cwy*(pww + pxw) + cxy*(pwx + pxx) + cyy*pxy + czy*pxz));

// c2

cb[ 8] = (float)((-4.0)*(cwx*pww + cxx*pwx + cyx*pwy + czx*pwz));

cb[ 9] = (float)(cyz*pwy + czz*pwz + cwz*(pww + pxw) + cxz*(pwx + pxx) + cyz*pxy + czz*pxz);

cb[10] = (float)(cwy*pww + cwy*pxw + cww*(pww + pxw) - cwx*(pww + pxw) + (cxw - cxx + cxy)*(pwx + pxx) + (cyw - cyx + cyy)*(pwy + pxy) + (czw - czx + czy)*(pwz + pxz));

cb[11] = (float)(0);

// c3

cb[12] = (float)((-4.0)*(cwx*pww + cxx*pwx + cyx*pwy + czx*pwz));

cb[13] = (float)((-2.0)*(cwz*pww + cxz*pwx + cyz*pwy + czz*pwz));

cb[14] = (float)(2.0*((-cww)*pww + cwx*pww - 2.0*cwy*pww - cxw*pwx + cxx*pwx - 2.0*cxy*pwx - cyw*pwy + cyx*pwy - 2.0*cyy*pwy - czw*pwz + czx*pwz - 2.0*czy*pwz + cwy*pyw + cxy*pyx + cyy*pyy + czy*pyz));

cb[15] = (float)(2.0*(cyx*pwy + czx*pwz + cwx*(pww - pyw) + cxx*(pwx - pyx) - cyx*pyy - czx*pyz));

// c4

cb[16] = (float)(4.0*(cwy*pww + cxy*pwx + cyy*pwy + czy*pwz));

cb[17] = (float)(cyz*pwy + czz*pwz + cwz*(pww - pyw) + cxz*(pwx - pyx) - cyz*pyy - czz*pyz);

cb[18] = (float)(cwy*pww + cww*(pww - pyw) - cwy*pyw + cwx*((-pww) + pyw) + (cxw - cxx + cxy)*(pwx - pyx) + (cyw - cyx + cyy)*(pwy - pyy) + (czw - czx + czy)*(pwz - pyz));

cb[19] = (float)(0); }