@@ -713,6 +713,7 @@ void resample(
     using span_alloc_t = agg::span_allocator<color_type>;
     using span_conv_alpha_t = span_conv_alpha<color_type>;
 
+    using nn_affine_interpolator_t = accurate_interpolator_affine_nn<>;
     using affine_interpolator_t = agg::span_interpolator_linear<>;
     using arbitrary_interpolator_t =
         agg::span_interpolator_adaptor<agg::span_interpolator_linear<>, lookup_distortion>;
@@ -774,10 +775,10 @@ void resample(
 
     if (params.interpolation == NEAREST) {
         if (params.is_affine) {
-            using span_gen_t = typename type_mapping_t::template span_gen_nn_type<image_accessor_t, affine_interpolator_t>;
+            using span_gen_t = typename type_mapping_t::template span_gen_nn_type<image_accessor_t, nn_affine_interpolator_t>;
             using span_conv_t = agg::span_converter<span_gen_t, span_conv_alpha_t>;
             using nn_renderer_t = agg::renderer_scanline_aa<renderer_t, span_alloc_t, span_conv_t>;
-            affine_interpolator_t interpolator(inverted);
+            nn_affine_interpolator_t interpolator(inverted);
             span_gen_t span_gen(input_accessor, interpolator);
             span_conv_t span_conv(span_gen, conv_alpha);
             nn_renderer_t nn_renderer(renderer, span_alloc, span_conv);
 
@@ -167,17 +167,12 @@ image_resample(py::array input_array,
 
         if (is_affine) {
             convert_trans_affine(transform, params.affine);
-            // If affine parameters will make subpixels visible, treat as nonaffine instead
-            if (params.affine.sx >= agg::image_subpixel_scale / 2 || params.affine.sy >= agg::image_subpixel_scale / 2) {
-                is_affine = false;
-                params.affine = agg::trans_affine();  // reset to identity affine parameters
-            }
-        }
-        if (!is_affine) {
+            params.is_affine = is_affine;
+        } else {
             transform_mesh = _get_transform_mesh(transform, output_array.shape());
             params.transform_mesh = transform_mesh.data();
+            params.is_affine = false;
         }
-        params.is_affine = is_affine;
     }
 
     if (auto resampler =
 
@@ -2,6 +2,7 @@
 #define MPL_AGG_WORKAROUND_H
 
 #include "agg_pixfmt_rgba.h"
+#include "agg_trans_affine.h"
 
 /**********************************************************************
  WORKAROUND: This class is to workaround a bug in Agg SVN where the
@@ -42,4 +43,87 @@ struct fixed_blender_rgba_plain : agg::conv_rgba_plain<ColorT, Order>
     }
 };
 
+
+/**********************************************************************
+ This class provides higher-accuracy nearest-neighbor interpolation for
+ affine transforms than span_interpolator_linear by using
+ floating-point-based interpolation instead of integer-based
+*/
+
+template<class Transformer = agg::trans_affine, unsigned SubpixelShift = 8>
+class accurate_interpolator_affine_nn
+{
+public:
+    typedef Transformer trans_type;
+
+    enum subpixel_scale_e
+    {
+        subpixel_shift = SubpixelShift,
+        subpixel_scale  = 1 << subpixel_shift
+    };
+
+    //--------------------------------------------------------------------
+    accurate_interpolator_affine_nn() {}
+    accurate_interpolator_affine_nn(trans_type& trans) : m_trans(&trans) {}
+    accurate_interpolator_affine_nn(trans_type& trans,
+                                    double x, double y, unsigned len) :
+        m_trans(&trans)
+    {
+        begin(x, y, len);
+    }
+
+    //----------------------------------------------------------------
+    const trans_type& transformer() const { return *m_trans; }
+    void transformer(trans_type& trans) { m_trans = &trans; }
+
+    //----------------------------------------------------------------
+    void begin(double x, double y, unsigned len)
+    {
+        m_len = len - 1;
+
+        m_stx1 = x;
+        m_sty1 = y;
+        m_trans->transform(&m_stx1, &m_sty1);
+        m_stx1 *= subpixel_scale;
+        m_sty1 *= subpixel_scale;
+
+        m_stx2 = x + m_len;
+        m_sty2 = y;
+        m_trans->transform(&m_stx2, &m_sty2);
+        m_stx2 *= subpixel_scale;
+        m_sty2 *= subpixel_scale;
+    }
+
+    //----------------------------------------------------------------
+    void resynchronize(double xe, double ye, unsigned len)
+    {
+        m_len = len - 1;
+
+        m_trans->transform(&xe, &ye);
+        m_stx2 = xe * subpixel_scale;
+        m_sty2 = ye * subpixel_scale;
+    }
+
+    //----------------------------------------------------------------
+    void operator++()
+    {
+        m_stx1 += (m_stx2 - m_stx1) / m_len;
+        m_sty1 += (m_sty2 - m_sty1) / m_len;
+        m_len--;
+    }
+
+    //----------------------------------------------------------------
+    void coordinates(int* x, int* y) const
+    {
+        // Truncate instead of round because this interpolator needs to
+        // match the definitions for nearest-neighbor interpolation
+        *x = int(m_stx1);
+        *y = int(m_sty1);
+    }
+
+private:
+    trans_type* m_trans;
+    unsigned m_len;
+    double m_stx1, m_sty1, m_stx2, m_sty2;
+};
 #endif