diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7508c06351deb3b0fefebec84fee5b8329c6f550..2e00f872f31f964c7370260ada1b55df868931f3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,6 +32,10 @@ When preparing for a bug fix release create a new 2nd heading above the Fixed
 heading to indicate that only the bug fixes and security fixes are in the bug fix
 release.
 -->
+## [Unreleased]
+
+### Added
+- Added a mutual information matcher [#559](https://github.com/USGS-Astrogeology/autocnet/pull/559)
 
 ## [0.6.0]
 
diff --git a/autocnet/matcher/mutual_information.py b/autocnet/matcher/mutual_information.py
new file mode 100644
index 0000000000000000000000000000000000000000..51b044f69590924dd1a048739e6e2263be290b7a
--- /dev/null
+++ b/autocnet/matcher/mutual_information.py
@@ -0,0 +1,132 @@
+from math import floor
+
+import numpy as np
+
+from scipy.ndimage.measurements import center_of_mass
+
+def mutual_information(t1, t2, **kwargs):
+    """
+    Computes the correlation coefficient between two images using a histogram
+    comparison (Mutual information for joint histograms). The corr_map coefficient
+    will be between 0 and 4
+
+    Parameters
+    ----------
+
+    t1 : ndarray
+         First image to use in the histogram comparison
+
+    t2 : ndarray
+         Second image to use in the histogram comparison
+
+    Returns
+    -------
+
+    : float
+      Correlation coefficient computed between the two images being compared
+      between 0 and 4
+
+    See Also
+    --------
+    numpy.histogram2d : for the kwargs that can be passed to the comparison
+    """
+    hgram, x_edges, y_edges = np.histogram2d(t1.ravel(),t2.ravel(), **kwargs)
+
+    # Convert bins counts to probability values
+    pxy = hgram / float(np.sum(hgram))
+    px = np.sum(pxy, axis=1) # marginal for x over y
+    py = np.sum(pxy, axis=0) # marginal for y over x
+    px_py = px[:, None] * py[None, :] # Broadcast to multiply marginals
+    # Now we can do the calculation using the pxy, px_py 2D arrays
+    nzs = pxy > 0 # Only non-zero pxy values contribute to the sum
+    return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
+
+def mutual_information_match(d_template, s_image, subpixel_size=3,
+                             func=None, **kwargs):
+    """
+    Applys the mutual information matcher function over a search image using a
+    defined template
+
+
+    Parameters
+    ----------
+    d_template : ndarray
+                 The input search template used to 'query' the destination
+                 image
+
+    s_image : ndarray
+              The image or sub-image to be searched
+
+    subpixel_size : int
+                    Subpixel area size to search for the center of mass
+                    calculation
+
+    func : function
+           Function object to be used to compute the histogram comparison
+
+    Returns
+    -------
+    x : float
+        The x offset
+
+    y : float
+        The y offset
+
+    max_corr : float
+               The strength of the correlation in the range [0, 4].
+
+    corr_map : ndarray
+               Map of corrilation coefficients when comparing the template to
+               locations within the search area
+    """
+
+    if func == None:
+        func = mutual_information
+
+    image_size = s_image.shape
+    template_size = d_template.shape
+
+    y_diff = image_size[0] - template_size[0]
+    x_diff = image_size[1] - template_size[1]
+
+    max_corr = -np.inf
+    corr_map = np.zeros((y_diff+1, x_diff+1))
+    max_i = -1  # y
+    max_j = -1  # x
+    for i in range(y_diff+1):
+        for j in range(x_diff+1):
+            sub_image = s_image[i:i+template_size[1],  # y
+                                j:j+template_size[0]]  # x
+            corr = func(sub_image, d_template, **kwargs)
+            if corr > max_corr:
+                max_corr = corr
+                max_i = i
+                max_j = j
+            corr_map[i, j] = corr
+
+    y, x = np.unravel_index(np.argmax(corr_map, axis=None), corr_map.shape)
+
+    upper = int(2 + floor(subpixel_size / 2))
+    lower = upper - 1
+    # x, y are the location of the upper left hand corner of the template in the image
+    area = corr_map[y-lower:y+upper,
+                    x-lower:x+upper]
+
+    # Compute the y, x shift (subpixel) using scipys center_of_mass function
+    cmass  = center_of_mass(area)
+
+    if area.shape != (subpixel_size+2, subpixel_size+2):
+        print("Max correlation is too close to the boundary.")
+        return None, None, 0, None
+
+    subpixel_y_shift = subpixel_size - 1 - cmass[0]
+    subpixel_x_shift = subpixel_size - 1 - cmass[1]
+
+    y += subpixel_y_shift
+    x += subpixel_x_shift
+
+    # Compute the idealized shift (image center)
+    y -= (s_image.shape[0] / 2) - (d_template.shape[0] / 2)
+    x -= (s_image.shape[1] / 2) - (d_template.shape[1] / 2)
+
+    return float(x), float(y), float(max_corr), corr_map
diff --git a/autocnet/matcher/tests/test_mutual_information.py b/autocnet/matcher/tests/test_mutual_information.py
new file mode 100644
index 0000000000000000000000000000000000000000..1b0a00a8877f941abbd6988d97d745465e83bd73
--- /dev/null
+++ b/autocnet/matcher/tests/test_mutual_information.py
@@ -0,0 +1,34 @@
+import math
+import os
+import sys
+import unittest
+from unittest.mock import patch
+
+import pytest
+import numpy as np
+
+from .. import mutual_information
+
+def test_good_mi():
+    test_image1 = np.array([[i for i in range(50)] for j in range(50)])
+    corrilation = mutual_information.mutual_information(test_image1, test_image1)
+    assert corrilation == pytest.approx(2.30258509299404)
+
+def test_bad_mi():
+    test_image1 = np.array([[i for i in range(50)] for j in range(50)])
+    test_image2 = np.ones((50, 50))
+    corrilation = mutual_information.mutual_information(test_image1, test_image2)
+    assert corrilation == pytest.approx(0)
+
+def test_mutual_information():
+    d_template = np.array([[i for i in range(50, 100)] for j in range(50)])
+    s_image = np.ones((100, 100))
+
+    s_image[25:75, 25:75] = d_template
+
+    x_offset, y_offset, max_corr, corr_map = mutual_information.mutual_information_match(d_template, s_image, bins=20)
+    assert x_offset == 0.01711861257171421
+    assert y_offset == 0.0
+    assert max_corr == 2.9755967600033015
+    assert corr_map.shape == (51, 51)
+    assert np.min(corr_map) >= 0.0
diff --git a/autocnet/matcher/tests/test_subpixel.py b/autocnet/matcher/tests/test_subpixel.py
index 6b2cf1b0250afb2cd6f054c32812cc2b362b3d48..9d386006330ad6796b85136c404b60ca632128a5 100644
--- a/autocnet/matcher/tests/test_subpixel.py
+++ b/autocnet/matcher/tests/test_subpixel.py
@@ -32,14 +32,6 @@ def iris_pair():
     rts_image = rescaled[:sizer, :sizec] 
     return image, rts_image
 
-
-@pytest.fixture
-def apollo_subsets():
-    # These need to be geodata sets or just use mocks...
-    arr1 = imread(get_path('AS15-M-0295_SML(1).png'))[100:201, 123:224]
-    arr2 = imread(get_path('AS15-M-0295_SML(2).png'))[235:336, 95:196]
-    return arr1, arr2
-
 @pytest.fixture
 def apollo_subsets():
     # These need to be geodata sets or just use mocks...