diff --git a/.vscode/settings.json b/.vscode/settings.json
index 03aa0229..df34dbec 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -31,4 +31,21 @@
     "editor.codeActionsOnSave": {
         "source.organizeImports": "explicit",
     },
+    "files.associations": {
+        "vector": "cpp",
+        "utility": "cpp",
+        "algorithm": "cpp",
+        "complex": "cpp",
+        "memory": "cpp",
+        "optional": "cpp",
+        "ostream": "cpp",
+        "sstream": "cpp",
+        "system_error": "cpp",
+        "type_traits": "cpp",
+        "xlocmon": "cpp",
+        "xmemory": "cpp",
+        "xstring": "cpp",
+        "xtr1common": "cpp",
+        "xutility": "cpp"
+    },
 }
\ No newline at end of file
diff --git a/setup.py b/setup.py
index ea0d0a6f..e3c5bb80 100644
--- a/setup.py
+++ b/setup.py
@@ -50,7 +50,7 @@ def compile_cython_extensions():
 # -------------------------------------------------------------------------------------
 
 MACROS_CPP = [("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")]
-CFLAGS_CPP = ["/EHsc"] if is_msvc() else ["-Wall"]
+CFLAGS_CPP = ["/EHsc", "/fp:fast"] if is_msvc() else ["-Wall"]
 if platform.system() == "Darwin":
     CFLAGS_CPP += ["-std=c++11"]
 
diff --git a/src/arrays.hpp b/src/arrays.hpp
index aea3a130..89f8f012 100644
--- a/src/arrays.hpp
+++ b/src/arrays.hpp
@@ -23,15 +23,15 @@ class PixelIterator
         check_img_ptr("pixeliter:", cur, out, img);
         return *cur;
     }
-    value_type &operator()(int dx, int dy)
+    value_type &operator()(npy_intp dx, npy_intp dy)
     {
         return *(cur + dy * img.si + dx * img.sj);
     }
-    void move(int dx, int dy)
+    void move(npy_intp dx, npy_intp dy)
     {
         cur += dy * img.si + dx * img.sj;
     }
-    void moveto(int x, int y)
+    void moveto(npy_intp x, npy_intp y)
     {
         cur = img.base + y * img.si + x * img.sj;
     }
@@ -52,19 +52,19 @@ class Array1D
     public:
         iterator() : pos(0L), stride(0) {}
         iterator(const Array1D &arr) : pos(arr.base), stride(arr.si) {}
-        iterator(const iterator &it, int n = 0) : pos(it.pos), stride(it.stride) { *this += n; }
+        iterator(const iterator &it, npy_intp n = 0) : pos(it.pos), stride(it.stride) { *this += n; }
         T &operator*() { return *pos; }
         const T &operator*() const { return *pos; }
-        T &operator[](int n) { return *(pos + n * stride); }
-        const T &operator[](int n) const { return *(pos + n * stride); }
-        iterator &operator+=(int n)
+        T &operator[](npy_intp n) { return *(pos + n * stride); }
+        const T &operator[](npy_intp n) const { return *(pos + n * stride); }
+        iterator &operator+=(npy_intp n)
         {
             pos += stride * n;
             return *this;
         }
-        int operator-(const iterator &it) { return (pos - it.pos) / stride; }
-        iterator operator+(int n) { return iterator(*this, n); }
-        iterator operator-(int n) { return iterator(*this, -n); }
+        npy_intp operator-(const iterator &it) { return (pos - it.pos) / stride; }
+        iterator operator+(npy_intp n) { return iterator(*this, n); }
+        iterator operator-(npy_intp n) { return iterator(*this, -n); }
         iterator &operator=(const iterator &it)
         {
             pos = it.pos;
@@ -99,7 +99,7 @@ class Array1D
 
     protected:
         T *pos;
-        int stride;
+        npy_intp stride;
     };
     Array1D() {}
     Array1D(PyArrayObject *arr)
@@ -109,8 +109,8 @@ class Array1D
         si = PyArray_STRIDE(arr, 0) / sizeof(value_type);
     }
 
-    Array1D(value_type *_base, int _ni, int _si) : base(_base), ni(_ni),
-                                                   si(_si / sizeof(value_type))
+    Array1D(value_type *_base, npy_intp _ni, npy_intp _si) : base(_base), ni(_ni),
+                                                             si(_si / sizeof(value_type))
     {
     }
     iterator begin() { return iterator(*this); }
@@ -120,7 +120,7 @@ class Array1D
         it += ni;
         return it;
     }
-    void init(value_type *_base, int _ni, int _si)
+    void init(value_type *_base, npy_intp _ni, npy_intp _si)
     {
         base = _base;
         ni = _ni;
@@ -128,12 +128,12 @@ class Array1D
     }
 
     // Pixel accessors
-    value_type &value(int x)
+    value_type &value(npy_intp x)
     {
         check("array1d:", x, ni, outside);
         return *(base + x * si);
     }
-    const value_type &value(int x) const
+    const value_type &value(npy_intp x) const
     {
         check("array1d:", x, ni, outside);
         return *(base + x * si);
@@ -142,8 +142,8 @@ class Array1D
 public:
     value_type outside;
     value_type *base;
-    int ni; // dimensions
-    int si; // strides in sizeof(value_type)
+    npy_intp ni; // dimensions
+    npy_intp si; // strides in sizeof(value_type)
 };
 
 template <class T>
@@ -162,11 +162,11 @@ class Array2D
         sj = PyArray_STRIDE(arr, 1) / sizeof(value_type);
     }
 
-    Array2D(value_type *_base, int _ni, int _nj, int _si, int _sj) : base(_base), ni(_ni), nj(_nj),
-                                                                     si(_si / sizeof(value_type)), sj(_sj / sizeof(value_type))
+    Array2D(value_type *_base, npy_intp _ni, npy_intp _nj, npy_intp _si, npy_intp _sj) : base(_base), ni(_ni), nj(_nj),
+                                                                                         si(_si / sizeof(value_type)), sj(_sj / sizeof(value_type))
     {
     }
-    void init(value_type *_base, int _ni, int _nj, int _si, int _sj)
+    void init(value_type *_base, npy_intp _ni, npy_intp _nj, npy_intp _si, npy_intp _sj)
     {
         base = _base;
         ni = _ni;
@@ -176,13 +176,13 @@ class Array2D
     }
 
     // Pixel accessors
-    value_type &value(int x, int y)
+    value_type &value(npy_intp x, npy_intp y)
     {
         check("array2d x:", x, nj, outside);
         check("array2d y:", y, ni, outside);
         return *(base + x * sj + y * si);
     }
-    const value_type &value(int x, int y) const
+    const value_type &value(npy_intp x, npy_intp y) const
     {
         check("array2d x:", x, nj, outside);
         check("array2d y:", y, ni, outside);
@@ -192,8 +192,8 @@ class Array2D
 public:
     value_type outside;
     value_type *base;
-    int ni, nj; // dimensions
-    int si, sj; // strides in sizeof(value_type)
+    npy_intp ni, nj; // dimensions
+    npy_intp si, sj; // strides in sizeof(value_type)
 };
 
 template <class Image>
diff --git a/src/pcolor.cpp b/src/pcolor.cpp
index 1db14a4d..40d63612 100644
--- a/src/pcolor.cpp
+++ b/src/pcolor.cpp
@@ -9,7 +9,6 @@
 #include <numpy/arrayobject.h>
 #ifdef _MSC_VER
 #include <float.h>
-#pragma fenv_access(on)
 #else
 #include <fenv.h>
 #endif
@@ -50,19 +49,19 @@ using std::vector;
 }
 #endif
 
-static bool vert_line(double _x0, double _y0, double _x1, double _y1, int NX,
-                      vector<int> &imin, vector<int> &imax,
-                      bool draw, npy_uint32 col, Array2D<npy_uint32> &D)
+static bool vert_line(double _x0, double _y0, double _x1, double _y1, npy_intp NX,
+                      vector<npy_intp> &imin, vector<npy_intp> &imax,
+                      bool draw, npy_intp col, Array2D<npy_intp> &D)
 {
-    int x0 = lrint(_x0);
-    int y0 = lrint(_y0);
-    int x1 = lrint(_x1);
-    int y1 = lrint(_y1);
-    int dx = abs(x1 - x0);
-    int dy = abs(y1 - y0);
-    int sx, sy;
-    int NY = imin.size() - 1;
-    int err, e2;
+    npy_intp x0 = lrint(_x0);
+    npy_intp y0 = lrint(_y0);
+    npy_intp x1 = lrint(_x1);
+    npy_intp y1 = lrint(_y1);
+    npy_intp dx = abs(x1 - x0);
+    npy_intp dy = abs(y1 - y0);
+    npy_intp sx, sy;
+    npy_intp NY = imin.size() - 1;
+    npy_intp err, e2;
     bool visible = false;
     NX = NX - 1;
     if (x0 < x1)
@@ -79,8 +78,8 @@ static bool vert_line(double _x0, double _y0, double _x1, double _y1, int NX,
     {
         if (y0 >= 0 && y0 <= NY)
         {
-            int _min = min(imin[y0], x0);
-            int _max = max(imax[y0], x0);
+            npy_intp _min = min(imin[y0], x0);
+            npy_intp _max = max(imax[y0], x0);
             if (draw)
             {
                 if (x0 >= 0 && x0 <= NX)
@@ -88,8 +87,8 @@ static bool vert_line(double _x0, double _y0, double _x1, double _y1, int NX,
                     D.value(x0, y0) = col;
                 }
             }
-            imin[y0] = max(0, _min);
-            imax[y0] = min(NX, _max);
+            imin[y0] = max<npy_intp>(0, _min);
+            imax[y0] = min<npy_intp>(NX, _max);
             if (_min <= NX && _max >= 0)
             {
                 visible = true;
@@ -118,19 +117,19 @@ struct QuadHelper
     const Array2D<T> &X;
     const Array2D<T> &Y;
     const Array2D<T> &Z;
-    Array2D<npy_uint32> &D;
+    Array2D<npy_intp> &D;
     LutScale<T, npy_uint32> &scale;
     double x1, x2, y1, y2, m_dx, m_dy;
     npy_uint32 bgcolor;
     bool border;
     bool flat;
     double uflat, vflat;
-    int ixmin, ixmax, iymin, iymax;
+    npy_intp ixmin, ixmax, iymin, iymax;
 
     QuadHelper(const Array2D<T> &X_,
                const Array2D<T> &Y_,
                const Array2D<T> &Z_,
-               Array2D<npy_uint32> &D_,
+               Array2D<npy_intp> &D_,
                LutScale<T, npy_uint32> &scale_,
                double x1_, double x2_, double y1_, double y2_,
                bool _border, bool _flat,
@@ -146,8 +145,8 @@ struct QuadHelper
 
     void draw_triangles()
     {
-        int i, j;
-        vector<int> imin, imax;
+        npy_intp i, j;
+        vector<npy_intp> imin, imax;
         imin.resize(D.ni);
         imax.resize(D.ni);
         ixmin = D.nj;
@@ -163,10 +162,10 @@ struct QuadHelper
         }
     }
 
-    void draw_quad(int qi, int qj,
-                   vector<int> &imin, vector<int> &imax)
+    void draw_quad(npy_intp qi, npy_intp qj,
+                   vector<npy_intp> &imin, vector<npy_intp> &imax)
     {
-        int i, j;
+        npy_intp i, j;
         double u, v;
         double v0, v1, v2, v3, v4;
 
@@ -181,8 +180,8 @@ struct QuadHelper
         double ymin = min(ay, min(by, min(cy, dy)));
         double ymax = max(ay, max(by, max(cy, dy)));
 
-        int i0 = int(ymin + .5);
-        int i1 = int(ymax + .5);
+        npy_intp i0 = int(ymin + .5);
+        npy_intp i1 = int(ymax + .5);
         //	printf("Quads: i=%d->%d\n", i0, i1);
 
         if (i0 < 0)
@@ -237,7 +236,7 @@ struct QuadHelper
             // XXX Color = Alpha
             return;
         }
-        int dm = 0, dM = 0;
+        npy_intp dm = 0, dM = 0;
         if (border)
         {
             dm = 1;
@@ -251,8 +250,8 @@ struct QuadHelper
         {
             ixmin = min(ixmin, imin[i]);
             ixmax = max(ixmax, imax[i]);
-            int jmin = max(0, imin[i]) + dm;
-            int jmax = min(imax[i], D.nj - 1) + dM;
+            npy_intp jmin = max<npy_intp>(0, imin[i]) + dm;
+            npy_intp jmax = min<npy_intp>(imax[i], D.nj - 1) + dM;
             for (j = jmin; j <= jmax; ++j)
             {
                 if (!flat)
@@ -436,7 +435,7 @@ PyObject *py_scale_quads(PyObject *self, PyObject *args)
     Array2D<double> X(p_src_x), Y(p_src_y), Z(p_src_z);
     /* Destination is RGB */
     unsigned long bg = 0;
-    Array2D<npy_uint32> dest(p_dst);
+    Array2D<npy_intp> dest(p_dst);
     if (apply_bg)
     {
 #if PY_MAJOR_VERSION >= 3
@@ -483,9 +482,9 @@ PyObject *py_vert_line(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_TypeError, "imin, imax must be int ndarray");
         return NULL;
     }
-    Array1D<int> pmin(p_min), pmax(p_max);
-    vector<int> imin, imax;
-    int nx = int(max(y0, y1)) + 1;
+    Array1D<npy_intp> pmin(p_min), pmax(p_max);
+    vector<npy_intp> imin, imax;
+    npy_intp nx = int(max(y0, y1)) + 1;
     if (pmin.ni < nx || pmax.ni < nx)
     {
         PyErr_SetString(PyExc_TypeError, "imin, imax not large enough");
@@ -497,14 +496,14 @@ PyObject *py_vert_line(PyObject *self, PyObject *args)
     }
     imin.resize(nx);
     imax.resize(nx);
-    for (int i = 0; i < nx; ++i)
+    for (npy_intp i = 0; i < nx; ++i)
     {
         imin[i] = pmin.value(i);
         imax[i] = pmax.value(i);
     }
-    Array2D<npy_uint32> dummy;
+    Array2D<npy_intp> dummy;
     vert_line(x0, y0, x1, y1, xmax, imin, imax, false, 0, dummy);
-    for (int i = 0; i < nx; ++i)
+    for (npy_intp i = 0; i < nx; ++i)
     {
         pmin.value(i) = imin[i];
         pmax.value(i) = imax[i];
diff --git a/src/points.hpp b/src/points.hpp
index 6b8ca9d8..d88c3d6c 100644
--- a/src/points.hpp
+++ b/src/points.hpp
@@ -21,8 +21,8 @@ class Point
     Point() : _ix(0), _iy(0), _x(0.0), _y(0.0) {}
     real x() const { return _x; }
     real y() const { return _y; }
-    int ix() const { return _ix; }
-    int iy() const { return _iy; }
+    npy_intp ix() const { return _ix; }
+    npy_intp iy() const { return _iy; }
 
     void setx(real x)
     {
@@ -53,7 +53,7 @@ class Point
     }
 
 protected:
-    int _ix, _iy;
+    npy_intp _ix, _iy;
     real _x, _y;
 };
 
@@ -65,7 +65,7 @@ class Point2DRectilinear : public Point
 public:
     Point2DRectilinear() : _insidex(true), _insidey(true) {}
     bool inside() const { return _insidex && _insidey; }
-    void testx(int _m, int _M)
+    void testx(npy_intp _m, npy_intp _M)
     {
         if (_ix < _m || _ix >= _M)
         {
@@ -76,7 +76,7 @@ class Point2DRectilinear : public Point
             _insidex = true;
         }
     }
-    void testy(int _m, int _M)
+    void testy(npy_intp _m, npy_intp _M)
     {
         if (_iy < _m || _iy >= _M)
         {
@@ -104,7 +104,7 @@ class Point2DRectilinear : public Point
 };
 
 /* A special transformation operation that transforms
-   i,j (int) coordinates with a translation and scale
+   i,j (npy_intp) coordinates with a translation and scale
 */
 class ScaleTransform
 {
@@ -112,11 +112,11 @@ class ScaleTransform
     typedef Point2DRectilinear point;
     typedef point::real real;
 
-    ScaleTransform(int _nx, int _ny,
+    ScaleTransform(npy_intp _nx, npy_intp _ny,
                    real _x0, real _y0,
                    real _dx, real _dy) : nx(_nx), ny(_ny), x0(_x0), y0(_y0), dx(_dx), dy(_dy) {}
 
-    void set(point &p, int x, int y) const
+    void set(point &p, npy_intp x, npy_intp y) const
     {
         p.set(x0 + x * dx, y0 + y * dy);
         p.testx(0, nx);
@@ -134,7 +134,7 @@ class ScaleTransform
     }
 
 public:
-    int nx, ny;
+    npy_intp nx, ny;
     real x0, y0;
     real dx, dy;
 };
@@ -155,7 +155,7 @@ class Point2D : public Point
         return *this;
     }
 
-    void test(int _xm, int _xM, int _ym, int _yM)
+    void test(npy_intp _xm, npy_intp _xM, npy_intp _ym, npy_intp _yM)
     {
         if ((_ix < _xm) || (_ix >= _xM) || (_iy < _ym) || (_iy >= _yM))
         {
@@ -175,7 +175,7 @@ class LinearTransform
     typedef Point2D point;
     typedef point::real real;
 
-    LinearTransform(int _nx, int _ny,
+    LinearTransform(npy_intp _nx, npy_intp _ny,
                     real _x0, real _y0,
                     real _xx, real _xy,
                     real _yx, real _yy) : nx(_nx), ny(_ny),
@@ -184,7 +184,7 @@ class LinearTransform
                                           yx(_yx), yy(_yy)
     {
     }
-    void set(point &p, int x, int y) const
+    void set(point &p, npy_intp x, npy_intp y) const
     {
         p.set(x0 + x * xx + y * xy,
               y0 + x * yx + y * yy);
@@ -202,7 +202,7 @@ class LinearTransform
     }
 
 public:
-    int nx, ny;
+    npy_intp nx, ny;
     real x0, y0;
     real xx, xy, yx, yy;
 };
@@ -296,7 +296,7 @@ class XYTransform
     typedef Point2DAxis<axis_type> point;
     typedef typename point::real real;
 
-    XYTransform(int _nx, int _ny,
+    XYTransform(npy_intp _nx, npy_intp _ny,
                 const axis_type &_ax,
                 const axis_type &_ay,
                 real _x0, real _y0,
@@ -327,7 +327,7 @@ class XYTransform
             p._insidey = true;
         }
     }
-    void set(point &p, int x, int y) const
+    void set(point &p, npy_intp x, npy_intp y) const
     {
         p.set(ax, x0 + x * dx, ay, y0 + y * dy);
         testx(p);
@@ -345,7 +345,7 @@ class XYTransform
     }
 
 public:
-    int nx, ny;
+    npy_intp nx, ny;
     const axis_type &ax;
     const axis_type &ay;
     real x0, y0;
diff --git a/src/scaler.cpp b/src/scaler.cpp
index 10e56445..8342bc22 100644
--- a/src/scaler.cpp
+++ b/src/scaler.cpp
@@ -9,7 +9,6 @@
 #include <numpy/arrayobject.h>
 #ifdef _MSC_VER
 #include <float.h>
-#pragma fenv_access(on)
 #define FE_TOWARDZERO _RC_CHOP
 #define fegetround() (_controlfp(0, 0) & _MCW_RC)
 int fesetround(int r)
@@ -86,8 +85,8 @@ struct NearestInterpolation
 {
     T operator()(const Array2D<T> &src, const TR &tr, const typename TR::point &p)
     {
-        int nx = p.ix();
-        int ny = p.iy();
+        npy_intp nx = p.ix();
+        npy_intp ny = p.iy();
         return src.value(nx, ny);
     }
 };
@@ -97,8 +96,8 @@ struct LinearInterpolation
 {
     T operator()(const Array2D<T> &src, const TR &tr, const typename TR::point &p)
     {
-        int nx = p.ix();
-        int ny = p.iy();
+        npy_intp nx = p.ix();
+        npy_intp ny = p.iy();
         double v = src.value(nx, ny);
         double a = 0;
 
@@ -132,8 +131,8 @@ struct LinearInterpolation<npy_uint32, TR>
     npy_uint32 operator()(const Array2D<npy_uint32> &src, const TR &tr, const typename TR::point &p)
     {
         int k;
-        int nx = p.ix();
-        int ny = p.iy();
+        npy_intp nx = p.ix();
+        npy_intp ny = p.iy();
         rgba_t p1, p2, p3, p4, r;
         p1.v = src.value(nx, ny);
         float v[4], v2[4];
@@ -224,8 +223,8 @@ struct LinearInterpolation<T, XYScale>
 {
     T operator()(const Array2D<T> &src, const XYScale &tr, const typename XYScale::point &p)
     {
-        int nx = p.ix();
-        int ny = p.iy();
+        npy_intp nx = p.ix();
+        npy_intp ny = p.iy();
         double v = src.value(nx, ny);
         double a = 0;
 
@@ -274,7 +273,7 @@ struct SubSampleInterpolation
     }
     T operator()(const Array2D<T> &src, const TR &tr, const typename TR::point &p0)
     {
-        int i, j;
+        npy_intp i, j;
         typename TR::point p, p1;
         typename num_trait<T>::large_type value = 0;
         typename num_trait<T>::large_type count = 0, msk, val;
@@ -300,9 +299,9 @@ struct SubSampleInterpolation
         }
         // printf("%d/%d\n", (int)value, (int)count);
         if (count)
-            return value / count;
+            return static_cast<T>(value / count);
         else
-            return value;
+            return static_cast<T>(value);
     }
     typename TR::real ki, kj;
     const Array2D<T> &mask;
@@ -314,7 +313,7 @@ void _scale_rgb(DEST &dest,
                 int dx1, int dy1, int dx2, int dy2,
                 Interpolation &interpolate)
 {
-    int i, j;
+    npy_intp i, j;
     ST val;
     int round = fegetround();
     PixelIterator<DEST> it(dest);
@@ -453,8 +452,8 @@ static bool check_transform(PyArrayObject *p_tr)
         PyErr_SetString(PyExc_TypeError, "transform data type must be float");
         return false;
     }
-    int ni = PyArray_DIM(p_tr, 0);
-    int nj = PyArray_DIM(p_tr, 1);
+    npy_intp ni = PyArray_DIM(p_tr, 0);
+    npy_intp nj = PyArray_DIM(p_tr, 1);
 
     if (ni != 3 || nj != 3)
     {
@@ -620,8 +619,8 @@ template <class Params>
 static PyObject *dispatch_source(Params &p)
 {
     bool ok;
-    int dni = PyArray_DIM(p.p_dst, 0);
-    int dnj = PyArray_DIM(p.p_dst, 1);
+    npy_intp dni = PyArray_DIM(p.p_dst, 0);
+    npy_intp dnj = PyArray_DIM(p.p_dst, 1);
 
     if (!PyArg_ParseTuple(p.p_dst_data, "iiii",
                           &p.dx1, &p.dy1, &p.dx2, &p.dy2))
@@ -633,8 +632,8 @@ static PyObject *dispatch_source(Params &p)
         swap(p.dx1, p.dx2);
     if (p.dy2 < p.dy1)
         swap(p.dy1, p.dy2);
-    check_image_bounds(dni, dnj, p.dx1, p.dy1);
-    check_image_bounds(dni, dnj, p.dx2, p.dy2);
+    check_image_bounds(static_cast<int>(dni), static_cast<int>(dnj), p.dx1, p.dy1);
+    check_image_bounds(static_cast<int>(dni), static_cast<int>(dnj), p.dx2, p.dy2);
 
     switch (PyArray_TYPE(p.p_src))
     {
@@ -717,10 +716,10 @@ static PyObject *py_scale_xy(PyObject *self, PyObject *args)
     {
         return NULL;
     }
-    int ni = PyArray_DIM(p_src, 0);
-    int nj = PyArray_DIM(p_src, 1);
-    int dni = PyArray_DIM(p_dst, 0);
-    int dnj = PyArray_DIM(p_dst, 1);
+    npy_intp ni = PyArray_DIM(p_src, 0);
+    npy_intp nj = PyArray_DIM(p_src, 1);
+    npy_intp dni = PyArray_DIM(p_dst, 0);
+    npy_intp dnj = PyArray_DIM(p_dst, 1);
     double dx = (x2 - x1) / dnj;
     double dy = (y2 - y1) / dni;
     Array1D<double> ax(p_ax), ay(p_ay);
@@ -754,8 +753,8 @@ static PyObject *py_scale_tr(PyObject *self, PyObject *args)
         return NULL;
     }
 
-    int ni = PyArray_DIM(p_src, 0);
-    int nj = PyArray_DIM(p_src, 1);
+    npy_intp ni = PyArray_DIM(p_src, 0);
+    npy_intp nj = PyArray_DIM(p_src, 1);
     Array2D<double> tr(p_tr);
     LinearTransform trans(nj, ni,
                           tr.value(2, 0), tr.value(2, 1), // x0, y0
@@ -793,10 +792,10 @@ static PyObject *py_scale_rect(PyObject *self, PyObject *args)
         return NULL;
     }
 
-    int ni = PyArray_DIM(p_src, 0);
-    int nj = PyArray_DIM(p_src, 1);
-    int dni = PyArray_DIM(p_dst, 0);
-    int dnj = PyArray_DIM(p_dst, 1);
+    npy_intp ni = PyArray_DIM(p_src, 0);
+    npy_intp nj = PyArray_DIM(p_src, 1);
+    npy_intp dni = PyArray_DIM(p_dst, 0);
+    npy_intp dnj = PyArray_DIM(p_dst, 1);
     double dx = (x2 - x1) / dnj;
     double dy = (y2 - y1) / dni;
     ScaleTransform trans(nj, ni, x1, y1, dx, dy);
diff --git a/src/scaler.hpp b/src/scaler.hpp
index 152bb864..c95878c0 100644
--- a/src/scaler.hpp
+++ b/src/scaler.hpp
@@ -23,7 +23,7 @@ struct Scaler
 template <class T>
 struct Scaler<T, true>
 {
-    typedef num_trait<fixed> trait;
+    typedef num_trait<T> trait;
 
     Scaler(double _a, double _b) : a(trait::fromdouble(_a)), b(trait::fromdouble(_b)) {}
     int scale(T x) const { return trait::toint(a * x + b); }
diff --git a/src/traits.hpp b/src/traits.hpp
index b9bf46c1..0eaf6d65 100644
--- a/src/traits.hpp
+++ b/src/traits.hpp
@@ -119,6 +119,28 @@ struct num_trait<fixed>
     static const bool is_integer = false;
 };
 
+template <>
+struct num_trait<npy_int64>
+{
+    typedef npy_int64 value_type;
+    typedef signed __int64 large_type;
+    static int toint(value_type v) { return (int)v; }
+    static value_type fromdouble(double v) { return (value_type)v; }
+
+    static const bool is_integer = true;
+};
+
+template <>
+struct num_trait<npy_uint64>
+{
+    typedef npy_uint64 value_type;
+    typedef unsigned __int64 large_type;
+    static int toint(value_type v) { return static_cast<int>(v); }
+    static value_type fromdouble(double v) { return static_cast<value_type>(v); }
+
+    static const bool is_integer = true;
+};
+
 template <class A>
 static void dispatch_array(int npy_type, A &algo)
 {