Wrong porosity value on large images
The function `ps.metrics.porosity()` returns an incorrect value if applied to a large image. This happens on Windows if the number of voxels in any phase exceeds the maximum allowable value for `np.int32` (2_147_483_647).

Possible fix:
```python
def porespy_porosity(im):
    im = np.array(im, dtype=int)
    Vp = np.sum(im == 1)
    Vs = np.sum(im == 0)
    e = Vp / (Vs + Vp)
    return e

def fixed_porosity(im):
    im = np.array(im, dtype=int)
    Vp = np.sum(im == 1, dtype = np.int64)
    Vs = np.sum(im == 0, dtype = np.int64)
    e = Vp / (Vs + Vp)
    return e

img = np.zeros((1500,1500,2000), dtype=int)
img[:,:,:1500] = 1
print(porespy_porosity(img))
print(fixed_porosity(img))
``` 
diff --git a/porespy/filters/_size_seq_satn.py b/porespy/filters/_size_seq_satn.py
index ed268366..6b83bc26 100644
--- a/porespy/filters/_size_seq_satn.py
+++ b/porespy/filters/_size_seq_satn.py
@@ -209,11 +209,12 @@ def seq_to_satn(seq, im=None, mode='drainage'):
         seq[uninvaded_mask] = 0
     seq = rankdata(seq, method='dense') - 1
     b = np.bincount(seq)
-    if (solid_mask.sum() > 0) or (uninvaded_mask.sum() > 0):
+    if (solid_mask.sum(dtype=np.int64) > 0) or \\
+            (uninvaded_mask.sum(dtype=np.int64) > 0):
         b[0] = 0
     c = np.cumsum(b)
     seq = np.reshape(seq, solid_mask.shape)
-    satn = c[seq]/(~solid_mask).sum()
+    satn = c[seq]/(seq.size - solid_mask.sum(dtype=np.int64))
     satn[solid_mask] = 0
     satn[uninvaded_mask] = -1
     return satn
diff --git a/porespy/filters/_snows.py b/porespy/filters/_snows.py
index 1da8995d..192c6082 100644
--- a/porespy/filters/_snows.py
+++ b/porespy/filters/_snows.py
@@ -422,7 +422,7 @@ def trim_saddle_points(peaks, dt, maxiter=20):
             if np.all(peaks_extended == peaks_i):
                 new_peaks[sx] += peaks_i
                 break  # Found a true peak
-            elif np.sum(peaks_extended * peaks_i) == 0:
+            elif np.sum(peaks_extended * peaks_i, dtype=np.int64) == 0:
                 break  # Found a saddle point
             peaks_i = peaks_extended
         if iters >= maxiter:
@@ -495,7 +495,7 @@ def trim_saddle_points_legacy(peaks, dt, maxiter=10):
             peaks_extended = (peaks_max == dt_i) * im_i
             if np.all(peaks_extended == peaks_i):
                 break  # Found a true peak
-            elif np.sum(peaks_extended * peaks_i) == 0:
+            elif np.sum(peaks_extended * peaks_i, dtype=np.int64) == 0:
                 peaks_i = False
                 break  # Found a saddle point
             # The following line was also missing from the original.  It has
diff --git a/porespy/generators/_imgen.py b/porespy/generators/_imgen.py
index 40907908..0ac6cfa8 100644
--- a/porespy/generators/_imgen.py
+++ b/porespy/generators/_imgen.py
@@ -409,8 +409,8 @@ def rsa(
         n_max = np.inf
     # Compute volume fraction info
     vf_final = volume_fraction
-    vf_start = im.sum() / im.size
-    vf_template = ps_round(r, ndim=im.ndim, smooth=smooth).sum() / im.size
+    vf_start = im.sum(dtype=np.int64) / im.size
+    vf_template = ps_round(r, ndim=im.ndim, smooth=smooth).sum(dtype=np.int64) / im.size
     logger.debug(f"Initial volume fraction: {vf_start}")
     # Dilate existing objects by strel to remove pixels near them
     # from consideration for sphere placement
@@ -592,7 +592,7 @@ def bundle_of_tubes(
                           offset=0.5*spacing,
                           spacing=spacing,
                           lattice='sc')
-    N = im.sum()
+    N = im.sum(dtype=np.int64)
     if distribution is None:
         # +1 below is because randint 4.X gives a max of 3
         distribution = spst.randint(low=3, high=int(spacing/2 + 1))
@@ -672,7 +672,7 @@ def polydisperse_spheres(shape: List[int],
     phi_desired = 1 - (1 - porosity) / (len(Rs))
     im = np.ones(shape, dtype=bool)
     for r in Rs:
-        phi_im = im.sum() / np.prod(shape)
+        phi_im = im.sum(dtype=np.int64) / np.prod(shape)
         phi_corrected = 1 - (1 - phi_desired) / phi_im
         temp = overlapping_spheres(shape=shape, r=r, porosity=phi_corrected)
         im = im * temp
@@ -979,8 +979,9 @@ def overlapping_spheres(
     shape = np.array(shape)
     if np.size(shape) == 1:
         shape = np.full((3, ), int(shape))
-    ndim = (shape != 1).sum()
-    s_vol = ps_disk(r).sum() if ndim == 2 else ps_ball(r).sum()
+    ndim = (shape != 1).sum(dtype=np.int64)
+    s_vol = ps_disk(r).sum(dtype=np.int64) if ndim == 2 \\
+        else ps_ball(r).sum(dtype=np.int64)
 
     bulk_vol = np.prod(shape)
     N = int(np.ceil((1 - porosity) * bulk_vol / s_vol))
@@ -992,7 +993,7 @@ def overlapping_spheres(
 
     def g(im):
         r"""Returns fraction of 0s, given a binary image"""
-        return 1 - im.sum() / np.prod(shape)
+        return 1 - im.sum(dtype=np.int64) / np.prod(shape)
 
     # # Newton's method for getting image porosity match the given
     # w = 1.0                         # Damping factor
@@ -1177,7 +1178,7 @@ def _cylinders(
     elif np.size(shape) == 2:
         raise Exception("2D cylinders don't make sense")
     # Find hypotenuse of domain from [0,0,0] to [Nx,Ny,Nz]
-    H = np.sqrt(np.sum(np.square(shape))).astype(int)
+    H = np.sqrt(np.sum(np.square(shape), dtype=np.int64)).astype(int)
     if length is None:  # Assume cylinders span domain if length not given
         length = 2 * H
     R = min(int(length / 2), 2 * H)  # Trim given length to 2H if too long
@@ -1349,7 +1350,7 @@ def cylinders(
     # Rough estimate of n_fibers
     n_fibers_added = 0
     # Calculate fraction of fibers to be added in each iteration.
-    subdif = 0.8 / np.sum(np.arange(1, maxiter) ** 2)
+    subdif = 0.8 / np.sum(np.arange(1, maxiter) ** 2, dtype=np.int64)
     fractions = [0.2]
     for i in range(1, maxiter):
         fractions.append(fractions[i - 1] + (maxiter - i) ** 2 * subdif)
diff --git a/porespy/generators/_pseudo_packings.py b/porespy/generators/_pseudo_packings.py
index 1510e0f1..1f42b131 100644
--- a/porespy/generators/_pseudo_packings.py
+++ b/porespy/generators/_pseudo_packings.py
@@ -221,7 +221,7 @@ def pseudo_electrostatic_packing(
     options = np.where(dt == 1)
     for _ in tqdm(range(maxiter), **settings.tqdm):
         hits = dt[options] < dtmax
-        if hits.sum() == 0:
+        if hits.sum(dtype=np.int64) == 0:
             if dt.min() == dtmax:
                 break
             options = np.where(dt == dt.min())
diff --git a/porespy/metrics/_funcs.py b/porespy/metrics/_funcs.py
index 4ceb5222..d1001a35 100644
--- a/porespy/metrics/_funcs.py
+++ b/porespy/metrics/_funcs.py
@@ -172,7 +172,7 @@ def representative_elementary_volume(im, npoints=1000):
         p = pads[i]
         new_s = extend_slice(s, shape=im.shape, pad=p)
         temp = im[new_s]
-        Vp = np.sum(temp)
+        Vp = np.sum(temp, dtype=np.int64)
         Vt = np.size(temp)
         porosity[i] = Vp / Vt
         volume[i] = Vt
@@ -225,9 +225,9 @@ def porosity(im):
     to view online example.
 
     """
-    im = np.array(im, dtype=int)
-    Vp = np.sum(im == 1)
-    Vs = np.sum(im == 0)
+    im = np.array(im, dtype=np.int64)
+    Vp = np.sum(im == 1, dtype=np.int64)
+    Vs = np.sum(im == 0, dtype=np.int64)
     e = Vp / (Vs + Vp)
     return e
 
@@ -263,7 +263,8 @@ def porosity_profile(im, axis=0):
     im = np.atleast_3d(im)
     a = set(range(im.ndim)).difference(set([axis]))
     a1, a2 = a
-    prof = np.sum(np.sum(im == 1, axis=a2), axis=a1) / (im.shape[a2] * im.shape[a1])
+    tmp = np.sum(np.sum(im == 1, axis=a2, dtype=np.int64), axis=a1, dtype=np.int64)
+    prof = tmp / (im.shape[a2] * im.shape[a1])
     return prof
 
 
@@ -533,8 +534,10 @@ def chord_length_distribution(im, bins=10, log=False, voxel_size=1,
         x = np.log10(x)
     if normalization == 'length':
         h = list(np.histogram(x, bins=bins, density=False))
-        h[0] = h[0] * (h[1][1:] + h[1][:-1]) / 2  # Scale bin heigths by length
-        h[0] = h[0] / h[0].sum() / (h[1][1:] - h[1][:-1])  # Normalize h[0] manually
+        # Scale bin heigths by length
+        h[0] = h[0] * (h[1][1:] + h[1][:-1]) / 2
+        # Normalize h[0] manually
+        h[0] = h[0] / h[0].sum(dtype=np.int64) / (h[1][1:] - h[1][:-1])
     elif normalization in ['number', 'count']:
         h = np.histogram(x, bins=bins, density=True)
     else:
@@ -765,12 +768,13 @@ def _radial_profile(autocorr, bins, pf=None, voxel_size=1):
     return tpcf
 
 
-@njit(parallel=True)
+@njit(parallel=False)
 def _get_radial_sum(dt, bins, bin_size, autocorr):
     radial_sum = np.zeros_like(bins[:-1])
     for i, r in enumerate(bins[:-1]):
         mask = (dt <= r) * (dt > (r - bin_size[i]))
-        radial_sum[i] = np.sum(np.ravel(autocorr)[np.ravel(mask)]) / np.sum(mask)
+        radial_sum[i] = np.sum(np.ravel(autocorr)[np.ravel(mask)], dtype=np.int64) \\
+            / np.sum(mask)
     return radial_sum
 
 
@@ -948,7 +952,8 @@ def phase_fraction(im, normed=True):
     labels = np.unique(im)
     results = {}
     for label in labels:
-        results[label] = np.sum(im == label) * (1 / im.size if normed else 1)
+        results[label] = np.sum(im == label, dtype=np.int64) * \\
+            (1 / im.size if normed else 1)
     return results
 
 
@@ -1054,7 +1059,8 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
                 r = sizes[mask][0]*voxel_size
                 pc = -2*sigma*np.cos(np.deg2rad(theta))/r
                 x.append(pc)
-                snwp = ((seq <= n)*(seq > 0)*(im == 1)).sum()/im.sum()
+                snwp = ((seq <= n)*(seq > 0) *
+                        (im == 1)).sum(dtype=np.int64)/im.sum(dtype=np.int64)
                 y.append(snwp)
         pc_curve = Results()
         pc_curve.pc = x
@@ -1072,7 +1078,7 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
                 r = n*voxel_size
                 pc = -2*sigma*np.cos(np.deg2rad(theta))/r
                 x.append(pc)
-                snwp = ((sizes >= n)*(im == 1)).sum()/im.sum()
+                snwp = ((sizes >= n)*(im == 1)).sum(dtype=np.int64)/im.sum(dtype=np.int64)
                 y.append(snwp)
         pc_curve = Results()
         pc_curve.pc = x
@@ -1088,10 +1094,10 @@ def pc_curve(im, sizes=None, pc=None, seq=None,
             # Add a point at begining to ensure curve starts a 0, if no residual
             Ps = np.hstack((Ps[0] - np.abs(Ps[0]/2), Ps))
         y = []
-        Vp = im.sum()
+        Vp = im.sum(dtype=np.int64)
         temp = pc[im]
         for p in tqdm(Ps, **settings.tqdm):
-            y.append((temp <= p).sum()/Vp)
+            y.append((temp <= p).sum(dtype=np.int64)/Vp)
         pc_curve = Results()
         pc_curve.pc = Ps
         pc_curve.snwp = y
@@ -1158,7 +1164,7 @@ def satn_profile(satn, s, axis=0, span=10, mode='tile'):
                 void = satn[i*span:(i+1)*span, ...] != 0
                 nwp = (satn[i*span:(i+1)*span, ...] < s) \\
                     *(satn[i*span:(i+1)*span, ...] > 0)
-                y[i] = nwp.sum()/void.sum()
+                y[i] = nwp.sum(dtype=np.int64)/void.sum(dtype=np.int64)
                 z[i] = i*span + (span-1)/2
         if mode == 'slide':
             y = np.zeros(int(satn.shape[0]-span))
@@ -1166,7 +1172,7 @@ def satn_profile(satn, s, axis=0, span=10, mode='tile'):
             for i in range(int(satn.shape[0]-span)):
                 void = satn[i:i+span, ...] != 0
                 nwp = (satn[i:i+span, ...] < s)*(satn[i:i+span, ...] > 0)
-                y[i] = nwp.sum()/void.sum()
+                y[i] = nwp.sum(dtype=np.int64)/void.sum(dtype=np.int64)
                 z[i] = i + (span-1)/2
         return z, y
 
diff --git a/porespy/metrics/_meshtools.py b/porespy/metrics/_meshtools.py
index dcd70bca..ad571763 100644
--- a/porespy/metrics/_meshtools.py
+++ b/porespy/metrics/_meshtools.py
@@ -63,7 +63,7 @@ def region_volumes(regions, mode='marching_cubes'):
         if mode == 'marching_cubes':
             vols[i] = mesh_volume(region)
         elif mode.startswith('voxel'):
-            vols[i] = region.sum()
+            vols[i] = region.sum(dtype=np.int64)
     return vols
 
 
diff --git a/porespy/networks/_funcs.py b/porespy/networks/_funcs.py
index d5dc52cc..484bb148 100644
--- a/porespy/networks/_funcs.py
+++ b/porespy/networks/_funcs.py
@@ -269,7 +269,7 @@ def generate_voxel_image(network, pore_shape="sphere", throat_shape="cylinder",
         logger.debug(f"Maximum dimension: {max_dim} voxels")
         im = _generate_voxel_image(
             network, pore_shape, throat_shape, max_dim=max_dim)
-        eps = im.astype(bool).sum() / np.prod(im.shape)
+        eps = im.astype(bool).sum(dtype=np.int64) / np.prod(im.shape)
         err = abs(1 - eps / eps_old)
         eps_old = eps
         max_dim = int(max_dim * 1.25)
diff --git a/porespy/networks/_getnet.py b/porespy/networks/_getnet.py
index 5ab27d35..244c9dfc 100644
--- a/porespy/networks/_getnet.py
+++ b/porespy/networks/_getnet.py
@@ -178,11 +178,11 @@ def regions_to_network(regions, phases=None, voxel_size=1, accuracy='standard'):
         p_phase[pore] = (phases[s]*pore_im).max()
         temp = np.vstack(np.where(sub_dt == sub_dt.max()))[:, 0]
         p_coords_dt_global[pore, :] = temp + s_offset
-        p_volume[pore] = np.sum(pore_im)
+        p_volume[pore] = np.sum(pore_im, dtype=np.int64)
         p_dia_local[pore] = 2*np.amax(pore_dt)
         p_dia_global[pore] = 2*np.amax(sub_dt)
         # The following is overwritten if accuracy is set to 'high'
-        p_area_surf[pore] = np.sum(pore_dt == 1)
+        p_area_surf[pore] = np.sum(pore_dt == 1, dtype=np.int64)
         im_w_throats = spim.binary_dilation(input=pore_im, structure=struc_elem(1))
         im_w_throats = im_w_throats*sub_im
         Pn = np.unique(im_w_throats)[1:] - 1
@@ -192,7 +192,7 @@ def regions_to_network(regions, phases=None, voxel_size=1, accuracy='standard'):
                 vx = np.where(im_w_throats == (j + 1))
                 t_dia_inscribed.append(2*np.amax(sub_dt[vx]))
                 # The following is overwritten if accuracy is set to 'high'
-                t_perimeter.append(np.sum(sub_dt[vx] < 2))
+                t_perimeter.append(np.sum(sub_dt[vx] < 2, dtype=np.int64))
                 # The following is overwritten if accuracy is set to 'high'
                 t_area.append(np.size(vx[0]))
                 p_area_surf[pore] -= np.size(vx[0])
@@ -238,7 +238,7 @@ def regions_to_network(regions, phases=None, voxel_size=1, accuracy='standard'):
     PT1 = PT1-p_dia_local[P12[:, 0]]/2*voxel_size
     PT2 = PT2-p_dia_local[P12[:, 1]]/2*voxel_size
     dist = (p_coords[P12[:, 0]] - p_coords[P12[:, 1]])*voxel_size
-    net['throat.direct_length'] = np.sqrt(np.sum(dist**2, axis=1))
+    net['throat.direct_length'] = np.sqrt(np.sum(dist**2, axis=1, dtype=np.int64))
     net['throat.perimeter'] = np.array(t_perimeter)*voxel_size
     if (accuracy == 'high') and (im.ndim == 2):
         msg = "accuracy='high' only available in 3D, reverting to 'standard'"
diff --git a/porespy/simulations/_dns.py b/porespy/simulations/_dns.py
index eea5a6b0..352599d7 100644
--- a/porespy/simulations/_dns.py
+++ b/porespy/simulations/_dns.py
@@ -59,18 +59,18 @@ def tortuosity_fd(im, axis, solver=None):
     openpnm_v3 = op.__version__.startswith('3')
 
     # Obtain original porosity
-    eps0 = im.sum() / im.size
+    eps0 = im.sum(dtype=np.int64) / im.size
 
     # Remove floating pores
     inlets = faces(im.shape, inlet=axis)
     outlets = faces(im.shape, outlet=axis)
     im = trim_nonpercolating_paths(im, inlets=inlets, outlets=outlets)
     # Check if porosity is changed after trimmimg floating pores
-    eps = im.sum() / im.size
+    eps = im.sum(dtype=np.int64) / im.size
     if not eps:
         raise Exception('No pores remain after trimming floating pores')
     if eps < eps0:  # pragma: no cover
-        logger.warning(f'Found non-percolating regions, were filled to percolate')
+        logger.warning('Found non-percolating regions, were filled to percolate')
 
     # Generate a Cubic network to be used as an orthogonal grid
     net = op.network.CubicTemplate(template=im, spacing=1.0)
diff --git a/porespy/tools/_funcs.py b/porespy/tools/_funcs.py
index a786642c..5ebf3c4a 100644
--- a/porespy/tools/_funcs.py
+++ b/porespy/tools/_funcs.py
@@ -1342,7 +1342,7 @@ def insert_cylinder(im, xyz0, xyz1, r):
     template = np.zeros(shape=shape_template)
 
     # Shortcut for orthogonal cylinders
-    if (xyz0 == xyz1).sum() == 2:
+    if (xyz0 == xyz1).sum(dtype=np.int64) == 2:
         unique_dim = [xyz0[i] != xyz1[i] for i in range(3)].index(True)
         shape_template[unique_dim] = 1
         template_2D = disk(radius=r).reshape(shape_template)
diff --git a/porespy/visualization/_funcs.py b/porespy/visualization/_funcs.py
index 74116331..4eecc771 100644
--- a/porespy/visualization/_funcs.py
+++ b/porespy/visualization/_funcs.py
@@ -201,8 +201,9 @@ def satn_to_panels(satn, im, bins=None, axis=0, slice=None, **kwargs):
         im_data.pop('vmax')
         [im_data.pop(i) for i in kwargs]
         ax[i // n][i % n].imshow(**im_data, vmax=2, **kwargs)
-        ax[i // n][i % n].set_title(str(np.around(temp.sum()/im.sum(),
-                                                  decimals=5)))
+        ax[i // n][i % n].set_title(
+            str(np.around(temp.sum(dtype=np.int64)/im.sum(dtype=np.int64),
+                          decimals=5)))
         temp_old = np.copy(temp)
     return fig, ax
 
diff --git a/porespy/visualization/_views.py b/porespy/visualization/_views.py
index 6c43882b..3c61324d 100644
--- a/porespy/visualization/_views.py
+++ b/porespy/visualization/_views.py
@@ -183,5 +183,5 @@ def xray(im, axis=0):  # pragma: no cover
         im = np.transpose(im, axes=[1, 0, 2])
     if axis == 2:
         im = np.transpose(im, axes=[2, 1, 0])
-    im = np.sum(im, axis=0)
+    im = np.sum(im, axis=0, dtype=np.int64)
     return im
['porespy/filters/_snows.py', 'porespy/generators/_pseudo_packings.py', 'porespy/metrics/_funcs.py', 'porespy/visualization/_views.py', 'porespy/tools/_funcs.py', 'porespy/networks/_funcs.py', 'porespy/filters/_size_seq_satn.py', 'porespy/metrics/_meshtools.py', 'porespy/networks/_getnet.py', 'porespy/generators/_imgen.py', 'porespy/visualization/_funcs.py', 'porespy/simulations/_dns.py']