diff --git a/invesalius_cy/cy_mesh.pyx b/invesalius_cy/cy_mesh.pyx
index 4c8b605..3983254 100644
--- a/invesalius_cy/cy_mesh.pyx
+++ b/invesalius_cy/cy_mesh.pyx
@@ -6,6 +6,7 @@
 #cython: cdivision=True
 #cython: nonecheck=False
 
+import os
 import sys
 import time
 cimport numpy as np
@@ -68,7 +69,7 @@ cdef class Mesh:
             self.faces = _faces
             self.normals = _normals
 
-            for i in xrange(_faces.shape[0]):
+            for i in range(_faces.shape[0]):
                 self.map_vface[self.faces[i, 1]].push_back(i)
                 self.map_vface[self.faces[i, 2]].push_back(i)
                 self.map_vface[self.faces[i, 3]].push_back(i)
@@ -201,12 +202,12 @@ cdef class Mesh:
             idfaces = self.get_faces_by_vertex(v_id)
             nf = idfaces.size()
 
-            for nid in xrange(nf):
+            for nid in range(nf):
                 f_id = deref(idfaces)[nid]
                 if status_f.find(f_id) == status_f.end():
                     status_f[f_id] = True
 
-                    for j in xrange(3):
+                    for j in range(3):
                         vj = self.faces[f_id, j+1]
                         if status_v.find(vj) == status_v.end():
                             status_v[vj] = True
@@ -255,7 +256,7 @@ cdef vector[weight_t]* calc_artifacts_weight(Mesh mesh, vector[vertex_id_t]& ver
         near_vertices = mesh.get_near_vertices_to_v(vi_id, tmax)
         nnv = near_vertices.size()
 
-        for j in xrange(nnv):
+        for j in range(nnv):
             vj_id = deref(near_vertices)[j]
             vj = &mesh.vertices[vj_id, 0]
 
@@ -269,13 +270,13 @@ cdef vector[weight_t]* calc_artifacts_weight(Mesh mesh, vector[vertex_id_t]& ver
 
         del near_vertices
 
-    #  for i in xrange(msize):
+    #  for i in range(msize):
         #  if mesh.is_border(i):
             #  deref(weights)[i] = 0.0
 
     #  cdef vertex_id_t v0, v1, v2
-    #  for i in xrange(mesh.faces.shape[0]):
-        #  for j in xrange(1, 4):
+    #  for i in range(mesh.faces.shape[0]):
+        #  for j in range(1, 4):
             #  v0 = mesh.faces[i, j]
             #  vi = &mesh.vertices[v0, 0]
             #  if mesh.is_border(v0):
@@ -361,7 +362,7 @@ cdef vector[vertex_id_t]* find_staircase_artifacts(Mesh mesh, double[3] stack_or
 
     nv = mesh.vertices.shape[0]
 
-    for v_id in xrange(nv):
+    for v_id in range(nv):
         max_z = -10000
         min_z = 10000
         max_y = -10000
@@ -372,7 +373,7 @@ cdef vector[vertex_id_t]* find_staircase_artifacts(Mesh mesh, double[3] stack_or
         f_ids = mesh.get_faces_by_vertex(v_id)
         nf = deref(f_ids).size()
 
-        for i in xrange(nf):
+        for i in range(nf):
             f_id = deref(f_ids)[i]
             normal = &mesh.normals[f_id, 0]
 
@@ -415,7 +416,7 @@ cdef void taubin_smooth(Mesh mesh, vector[weight_t]& weights, float l, float m, 
     nvertices = mesh.vertices.shape[0]
     cdef vector[Point] D = vector[Point](nvertices)
     cdef vertex_t* vi
-    for s in xrange(steps):
+    for s in range(steps):
         for i in prange(nvertices, nogil=True):
             D[i] = calc_d(mesh, i)
 
diff --git a/invesalius_cy/floodfill.pyx b/invesalius_cy/floodfill.pyx
index 4183459..91ed2e4 100644
--- a/invesalius_cy/floodfill.pyx
+++ b/invesalius_cy/floodfill.pyx
@@ -142,11 +142,11 @@ def floodfill_threshold(np.ndarray[image_t, ndim=3] data, list seeds, int t0, in
 
             out[z, y, x] = fill
 
-            for k in xrange(odz):
+            for k in range(odz):
                 zo = z + k - offset_z
-                for j in xrange(ody):
+                for j in range(ody):
                     yo = y + j - offset_y
-                    for i in xrange(odx):
+                    for i in range(odx):
                         if strct[k, j, i]:
                             xo = x + i - offset_x
                             if 0 <= xo < dx and 0 <= yo < dy and 0 <= zo < dz and out[zo, yo, xo] != fill and t0 <= data[zo, yo, xo] <= t1:
@@ -254,13 +254,13 @@ def fill_holes_automatically(np.ndarray[mask_t, ndim=3] mask, np.ndarray[np.uint
     dy = mask.shape[1]
     dx = mask.shape[2]
 
-    for z in xrange(dz):
-        for y in xrange(dy):
-            for x in xrange(dx):
+    for z in range(dz):
+        for y in range(dy):
+            for x in range(dx):
                 sizes[labels[z, y, x]] += 1
 
     #Checking if any hole will be filled
-    for i in xrange(nlabels + 1):
+    for i in range(nlabels + 1):
         if sizes[i] <= max_size:
             modified = True
 
@@ -268,8 +268,8 @@ def fill_holes_automatically(np.ndarray[mask_t, ndim=3] mask, np.ndarray[np.uint
         return 0
 
     for z in prange(dz, nogil=True):
-        for y in xrange(dy):
-            for x in xrange(dx):
+        for y in range(dy):
+            for x in range(dx):
                 if sizes[labels[z, y, x]] <= max_size:
                     mask[z, y, x] = 254
 
diff --git a/invesalius_cy/interpolation.pyx b/invesalius_cy/interpolation.pyx
index a314d32..cbcce91 100644
--- a/invesalius_cy/interpolation.pyx
+++ b/invesalius_cy/interpolation.pyx
@@ -172,28 +172,28 @@ cdef double lanczos3(image_t[:, :, :] V, double x, double y, double z) nogil:
     cdef int m, n, o
 
     m = 0
-    for k in xrange(zi, zf):
+    for k in range(zi, zf):
         n = 0
-        for j in xrange(yi, yf):
+        for j in range(yi, yf):
             lx = 0
-            for i in xrange(xi, xf):
+            for i in range(xi, xf):
                 lx += _G(V, i, j, k) * lanczos3_L(x - i, a)
             temp_x[m][n] = lx
             n += 1
         m += 1
 
     m = 0
-    for k in xrange(zi, zf):
+    for k in range(zi, zf):
         n = 0
         ly = 0
-        for j in xrange(yi, yf):
+        for j in range(yi, yf):
             ly += temp_x[m][n] * lanczos3_L(y - j, a)
             n += 1
         temp_y[m] = ly
         m += 1
 
     m = 0
-    for k in xrange(zi, zf):
+    for k in range(zi, zf):
         lz += temp_y[m] * lanczos3_L(z - k, a)
         m += 1
 
@@ -309,7 +309,7 @@ cdef void calc_coef_tricub(image_t[:, :, :] V, double x, double y, double z, dou
 
     for j in prange(64):
         coef[j] = 0.0
-        for i in xrange(64):
+        for i in range(64):
                 coef[j] += (temp[j][i] * _x[i])
 
 
@@ -328,9 +328,9 @@ cdef double tricub_interpolate(image_t[:, :, :] V, double x, double y, double z)
     cdef int yi = <int>floor(y)
     cdef int zi = <int>floor(z)
 
-    for i in xrange(4):
-        for j in xrange(4):
-            for k in xrange(4):
+    for i in range(4):
+        for j in range(4):
+            for k in range(4):
                 result += (coef[i+4*j+16*k] * ((x-xi)**i) * ((y-yi)**j) * ((z-zi)**k))
     # return V[<int>z, <int>y, <int>x]
     # with gil:
@@ -370,9 +370,9 @@ cdef double tricubicInterpolate(image_t[:, :, :] V, double x, double y, double z
 
     cdef int i, j, k
 
-    for i in xrange(4):
-        for j in xrange(4):
-            for k in xrange(4):
+    for i in range(4):
+        for j in range(4):
+            for k in range(4):
                 p[i][j][k] = _G(V, xi + i -1, yi + j -1, zi + k - 1)
 
     cdef double arr[4]
diff --git a/invesalius_cy/mips.pyx b/invesalius_cy/mips.pyx
index 93b2641..a0954c4 100644
--- a/invesalius_cy/mips.pyx
+++ b/invesalius_cy/mips.pyx
@@ -28,14 +28,14 @@ def lmip(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t tmin,
 
     # AXIAL
     if axis == 0:
-        for x in xrange(sx):
-            for y in xrange(sy):
+        for x in range(sx):
+            for y in range(sy):
                 max = image[0, y, x]
                 if max >= tmin and max <= tmax:
                     start = 1
                 else:
                     start = 0
-                for z in xrange(sz):
+                for z in range(sz):
                     if image[z, y, x] > max:
                         max = image[z, y, x]
 
@@ -49,14 +49,14 @@ def lmip(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t tmin,
 
     #CORONAL
     elif axis == 1:
-        for z in xrange(sz):
-            for x in xrange(sx):
+        for z in range(sz):
+            for x in range(sx):
                 max = image[z, 0, x]
                 if max >= tmin and max <= tmax:
                     start = 1
                 else:
                     start = 0
-                for y in xrange(sy):
+                for y in range(sy):
                     if image[z, y, x] > max:
                         max = image[z, y, x]
 
@@ -70,14 +70,14 @@ def lmip(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t tmin,
 
     #CORONAL
     elif axis == 2:
-        for z in xrange(sz):
-            for y in xrange(sy):
+        for z in range(sz):
+            for y in range(sy):
                 max = image[z, y, 0]
                 if max >= tmin and max <= tmax:
                     start = 1
                 else:
                     start = 0
-                for x in xrange(sx):
+                for x in range(sx):
                     if image[z, y, x] > max:
                         max = image[z, y, x]
 
@@ -164,11 +164,11 @@ def mida(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t wl,
     # AXIAL
     if axis == 0:
         for x in prange(sx, nogil=True):
-            for y in xrange(sy):
+            for y in range(sy):
                 fmax = 0.0
                 alpha_p = 0.0
                 colour_p = 0.0
-                for z in xrange(sz):
+                for z in range(sz):
                     vl = image[z, y, x]
                     fpi = 1.0/(max - min) * (vl - min)
                     if fpi > fmax:
@@ -198,11 +198,11 @@ def mida(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t wl,
     #CORONAL
     elif axis == 1:
         for z in prange(sz, nogil=True):
-            for x in xrange(sx):
+            for x in range(sx):
                 fmax = 0.0
                 alpha_p = 0.0
                 colour_p = 0.0
-                for y in xrange(sy):
+                for y in range(sy):
                     vl = image[z, y, x]
                     fpi = 1.0/(max - min) * (vl - min)
                     if fpi > fmax:
@@ -229,11 +229,11 @@ def mida(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t wl,
     #AXIAL
     elif axis == 2:
         for z in prange(sz, nogil=True):
-            for y in xrange(sy):
+            for y in range(sy):
                 fmax = 0.0
                 alpha_p = 0.0
                 colour_p = 0.0
-                for x in xrange(sx):
+                for x in range(sx):
                     vl = image[z, y, x]
                     fpi = 1.0/(max - min) * (vl - min)
                     if fpi > fmax:
diff --git a/invesalius_cy/transforms.pyx b/invesalius_cy/transforms.pyx
index 8c2c10a..3b43b80 100644
--- a/invesalius_cy/transforms.pyx
+++ b/invesalius_cy/transforms.pyx
@@ -116,23 +116,23 @@ def apply_view_matrix_transform(image_t[:, :, :] volume,
         f_interp = lanczos3
 
     if orientation == 'AXIAL':
-        for z in xrange(n, n+odz):
+        for z in range(n, n+odz):
             for y in prange(dy, nogil=True):
-                for x in xrange(dx):
+                for x in range(dx):
                     out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
             count += 1
 
     elif orientation == 'CORONAL':
-        for y in xrange(n, n+ody):
+        for y in range(n, n+ody):
             for z in prange(dz, nogil=True):
-                for x in xrange(dx):
+                for x in range(dx):
                     out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
             count += 1
 
     elif orientation == 'SAGITAL':
-        for x in xrange(n, n+odx):
+        for x in range(n, n+odx):
             for z in prange(dz, nogil=True):
-                for y in xrange(dy):
+                for y in range(dy):
                     out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
             count += 1
 
@@ -157,14 +157,14 @@ def convolve_non_zero(image_t[:, :, :] volume,
     skx = kernel.shape[2]
 
     for z in prange(sz, nogil=True):
-        for y in xrange(sy):
-            for x in xrange(sx):
+        for y in range(sy):
+            for x in range(sx):
                 if volume[z, y, x] != 0:
-                    for k in xrange(skz):
+                    for k in range(skz):
                         kz = z - skz // 2 + k
-                        for j in xrange(sky):
+                        for j in range(sky):
                             ky = y - sky // 2 + j
-                            for i in xrange(skx):
+                            for i in range(skx):
                                 kx = x - skx // 2 + i
 
                                 if 0 <= kz < sz and 0 <= ky < sy and 0 <= kx < sx:
--
libgit2 0.21.2