diff --git a/examples/test_metrics.py b/examples/test_metrics.py
index c9c5bd1..fe4eaa3 100644
--- a/examples/test_metrics.py
+++ b/examples/test_metrics.py
@@ -18,6 +18,7 @@
 SM = py_sod_metrics.Smeasure()
 EM = py_sod_metrics.Emeasure()
 MAE = py_sod_metrics.MAE()
+MSIOU = py_sod_metrics.MSIoU()
 
 sample_gray = dict(with_adaptive=True, with_dynamic=True)
 sample_bin = dict(with_adaptive=False, with_dynamic=False, with_binary=True, sample_based=True)
@@ -78,6 +79,7 @@
     SM.step(pred=pred, gt=mask)
     EM.step(pred=pred, gt=mask)
     MAE.step(pred=pred, gt=mask)
+    MSIOU.step(pred=pred, gt=mask)
     FMv2.step(pred=pred, gt=mask)
 
 fm = FM.get_results()["fm"]
@@ -85,12 +87,14 @@
 sm = SM.get_results()["sm"]
 em = EM.get_results()["em"]
 mae = MAE.get_results()["mae"]
+msiou = MSIOU.get_results()["msiou"]
 fmv2 = FMv2.get_results()
 
 curr_results = {
     "MAE": mae,
     "Smeasure": sm,
     "wFmeasure": wfm,
+    "MSIOU": msiou,
     # E-measure for sod
     "adpEm": em["adp"],
     "meanEm": em["curve"].mean(),
@@ -253,6 +257,7 @@
     },
     "v1_4_1": {
         "MAE": 0.03705558476661653,
+        "MSIOU": 0.8228002109838289,
         "Smeasure": 0.9029761578759272,
         "adpEm": 0.9408760066970617,
         "adpFm": 0.5816750824038355,
@@ -336,6 +341,9 @@ def test_wfm(self):
     def test_mae(self):
         self.assertEqual(curr_results["MAE"], self.default_results["MAE"])
 
+    def test_msiou(self):
+        self.assertEqual(curr_results["MSIOU"], self.default_results["MSIOU"])
+
     def test_fm(self):
         self.assertEqual(curr_results["adpFm"], self.default_results["adpFm"])
         self.assertEqual(curr_results["meanFm"], self.default_results["meanFm"])
diff --git a/py_sod_metrics/__init__.py b/py_sod_metrics/__init__.py
index d68ec80..1eb846f 100755
--- a/py_sod_metrics/__init__.py
+++ b/py_sod_metrics/__init__.py
@@ -15,6 +15,7 @@
     TNRHandler,
     TPRHandler,
 )
+from py_sod_metrics.multiscale_iou import MSIoU
 from py_sod_metrics.sod_metrics import (
     MAE,
     Emeasure,
diff --git a/py_sod_metrics/multiscale_iou.py b/py_sod_metrics/multiscale_iou.py
new file mode 100644
index 0000000..374bf4a
--- /dev/null
+++ b/py_sod_metrics/multiscale_iou.py
@@ -0,0 +1,105 @@
+import numpy as np
+from scipy import ndimage
+
+from .utils import TYPE
+
+
+class MSIoU:
+    def __init__(self):
+        """
+        Multiscale Intersection over Union (MS-IoU) metric.
+
+        ::
+
+            @inproceedings{ahmadzadehMultiscaleIOUMetric2021,
+                title = {Multiscale IOU: A Metric for Evaluation of Salient Object Detection with Fine Structures},
+                booktitle = {2021 IEEE International Conference on Image Processing (ICIP)},
+                author = {Ahmadzadeh, Azim and Kempton, Dustin J. and Chen, Yang and Angryk, Rafal A.},
+                year = {2021},
+                doi = {10.1109/ICIP42928.2021.9506337}
+            }
+        """
+        # The values of this collection determines the resolutions based on which MIoU is computed.
+        # It is set as the original implementation
+        self.cell_sizes = np.power(2, np.linspace(0, 9, num=10, dtype=int))
+        self.msious = []
+
+    def get_edge(self, mask: np.ndarray):
+        """Edge detection based on the `scipy.ndimage.sobel` function.
+
+        :param mask: a binary mask of an object whose edges are of interest.
+        :return: a binary mask of 1's as edges and 0's as background.
+        """
+        sx = ndimage.sobel(mask, axis=0, mode="constant")
+        sy = ndimage.sobel(mask, axis=1, mode="constant")
+        sob = np.hypot(sx, sy)
+        sob[sob > 0] = 1
+        sob[sob <= 0] = 0
+        return sob
+
+    def shrink_by_grid(self, image: np.ndarray, cell_size: int) -> np.ndarray:
+        """Box-counting after the zero padding if needed."""
+        h, w = image.shape[:2]
+
+        pad_h = h % cell_size
+        if pad_h != 0:
+            pad_h = cell_size - pad_h
+        pad_w = w % cell_size
+        if pad_w != 0:
+            pad_w = cell_size - pad_w
+        if pad_h != 0 or pad_w != 0:
+            image = np.pad(image, ((pad_h, 0), (pad_w, 0)), mode="constant", constant_values=0)
+
+        h = image.shape[0]
+        w = image.shape[1]
+        image = image.reshape(h // cell_size, cell_size, w // cell_size, cell_size)
+        image = image.sum(axis=(1, 3))
+        image[image > 0] = 1
+        return image
+
+    def cal_msiou(self, pred: np.ndarray, gt: np.ndarray) -> float:
+        """
+        Args:
+            pred (np.ndarray[bool]):
+            gt (np.ndarray[bool]):
+
+        Returns:
+            float:
+        """
+        pred = self.get_edge(pred)
+        gt = self.get_edge(gt)
+
+        ratios = []
+        for cell_size in self.cell_sizes:
+            s_pred = self.shrink_by_grid(pred, cell_size=cell_size)
+            s_gt = self.shrink_by_grid(gt, cell_size=cell_size)
+            numerator = np.logical_and(s_pred, s_gt).sum() + 1
+            denominator = s_gt.sum() + 1
+            ratios.append(numerator / denominator)
+
+        # Calculates area under the curves using Trapezoids.
+        msiou = np.trapz(y=ratios, dx=1 / (len(self.cell_sizes) - 1))
+        return msiou
+
+    def step(self, pred: np.ndarray, gt: np.ndarray):
+        """
+        computes the metric for the two given regions. Use 'miou/utils/mask_loader.py'
+        to properly load masks as binary arrays for `pred` and `gt`.
+
+        :param pred: mask segmentation of the reference (ground-truth) object.
+        :param gt: mask segmentation of the target (detected) object.
+        """
+        gt = gt > 128
+        pred = pred > 128
+
+        msiou = self.cal_msiou(pred, gt)
+        self.msious.append(msiou)
+        return msiou
+
+    def get_results(self) -> dict:
+        """Return the results about MS-IoU.
+
+        :return: dict(msiou=msiou)
+        """
+        msiou = np.mean(np.array(self.msious, TYPE))
+        return dict(msiou=msiou)