diff --git a/common/mathlib.cc b/common/mathlib.cc
index 045aff65..12b3ee04 100644
--- a/common/mathlib.cc
+++ b/common/mathlib.cc
@@ -244,6 +244,27 @@ float Filter_Gaussian(float width, float height, float x, float y)
         * Gaussian1D(height, y, alpha);
 }
 
+// from https://en.wikipedia.org/wiki/Lanczos_resampling
+static float Lanczos1D(float x, float a)
+{
+    if (x == 0)
+        return 1;
+    
+    if (x < -a || x >= a)
+        return 0;
+    
+    float lanczos = (a * sinf(M_PI * x) * sinf(M_PI * x / a)) / (M_PI * M_PI * x * x);
+    return lanczos;
+}
+
+// from https://en.wikipedia.org/wiki/Lanczos_resampling#Multidimensional_interpolation
+float Lanczos2D(float x, float y, float a)
+{
+    float dist = sqrtf((x*x) + (y*y));
+    float lanczos = Lanczos1D(dist, a);
+    return lanczos;
+}
+
 using namespace glm;
 using namespace std;
 
diff --git a/include/common/mathlib.hh b/include/common/mathlib.hh
index 43bafa05..fe3e5fc8 100644
--- a/include/common/mathlib.hh
+++ b/include/common/mathlib.hh
@@ -302,6 +302,9 @@ int SampleCDF(const std::vector<float> &cdf, float sample);
 // width (height) are the filter "radius" (not "diameter")
 float Filter_Gaussian(float width, float height, float x, float y);
 
+// sqrt(x^2 + y^2) should be <= a, returns 0 outside that range.
+float Lanczos2D(float x, float y, float a);
+
 // glm geometry
 
 glm::vec3 GLM_FaceNormal(std::vector<glm::vec3> points);