diff --git a/README.md b/README.md
index c3e6d62..dfe91f7 100644
--- a/README.md
+++ b/README.md
@@ -111,6 +111,20 @@ Buckminster Fuller’s Airocean projection (also known as “Dymaxion”), based
 
 The Cahill-Keyes projection, designed by Gene Keyes (1975), is built on Bernard J. S. Cahill’s 1909 octant design. Implementation by Mary Jo Graça (2011), ported to D3 by Enrico Spinielli (2013).
 
+<a href="#geoImago" name="geoImago">#</a> d3.<b>geoImago</b>() [<>](https://github.com/d3/d3-geo-polygon/blob/master/src/imago.js "Source")
+
+[<img src="https://raw.githubusercontent.com/d3/d3-geo-polygon/master/img/imago.png" width="480" height="250">](https://kunimune.home.blog/2017/11/23/the-secrets-of-the-authagraph-revealed/)
+
+The Imago projection, engineered by Justin Kunimune (2017), is inspired by Hajime Narukawa’s AuthaGraph design (1999).
+
+<a href="#imago_k" name="imago_k">#</a> <i>imago</i>.<b>k</b>([<i>k</i>])
+
+Exponent. Useful values include 0.59 (default, minimizes angular distortion of the continents), 0.68 (gives the best approximation of the Authagraph) and 0.72 (prevents kinks in the graticule).
+
+<a href="#imago_shift" name="imago_cut">#</a> <i>imago</i>.<b>shift</b>([<i>shift</i>])
+
+Horizontal shift. Defaults to 1.16.
+
 <a href="#geoTetrahedralLee" name="geoTetrahedralLee">#</a> d3.<b>geoTetrahedralLee</b>() [<>](https://github.com/d3/d3-geo-polygon/blob/master/src/tetrahedralLee.js "Source")
 <br><a href="#geoLeeRaw" name="geoLeeRaw">#</a> d3.<b>geoLeeRaw</b>
 
diff --git a/img/imago.png b/img/imago.png
new file mode 100644
index 0000000..55affbc
Binary files /dev/null and b/img/imago.png differ
diff --git a/src/imago.js b/src/imago.js
new file mode 100644
index 0000000..ca621f0
--- /dev/null
+++ b/src/imago.js
@@ -0,0 +1,346 @@
+/*
+ * Imago projection, by Justin Kunimune
+ *
+ * Inspired by Hajime Narukawa’s AuthaGraph
+ *
+ */
+import {
+  abs,
+  acos,
+  asin,
+  atan,
+  atan2,
+  cos,
+  degrees,
+  epsilon,
+  floor,
+  halfPi,
+  pi,
+  pow,
+  sign,
+  sin,
+  sqrt,
+  tan
+} from "./math";
+import { geoProjectionMutator as projectionMutator } from "d3-geo";
+import { default as clipPolygon } from "./clip/polygon";
+import { solve } from "./newton.js";
+
+var hypot = Math.hypot;
+
+var ASIN_ONE_THD = asin(1 / 3),
+  centrums = [
+    [halfPi, 0, 0, -halfPi, 0, sqrt(3)],
+    [-ASIN_ONE_THD, 0, pi, halfPi, 0, -sqrt(3)],
+    [-ASIN_ONE_THD, (2 * pi) / 3, pi, (5 * pi) / 6, 3, 0],
+    [-ASIN_ONE_THD, (-2 * pi) / 3, pi, pi / 6, -3, 0]
+  ],
+  TETRAHEDRON_WIDE_VERTEX = {
+    sphereSym: 3,
+    planarSym: 6,
+    width: 6,
+    height: 2 * sqrt(3),
+    centrums,
+    rotateOOB: function(x, y, xCen, yCen) {
+      yCen * 0;
+      if (abs(x) > this.width / 2) return [2 * xCen - x, -y];
+      else return [-x, this.height * sign(y) - y];
+    },
+    inBounds: () => true
+  },
+  configuration = TETRAHEDRON_WIDE_VERTEX;
+
+export function imagoRaw(k) {
+  function faceProject(lon, lat) {
+    const tht = atan(((lon - asin(sin(lon) / sqrt(3))) / pi) * sqrt(12)),
+      p = (halfPi - lat) / atan(sqrt(2) / cos(lon));
+
+    return [(pow(p, k) * sqrt(3)) / cos(tht), tht];
+  }
+
+  function faceInverse(r, th) {
+    const l = solve(
+        l => atan(((l - asin(sin(l) / sqrt(3))) / pi) * sqrt(12)),
+        th,
+        th / 2
+      ),
+      R = r / (sqrt(3) / cos(th));
+    return [halfPi - pow(R, 1 / k) * atan(sqrt(2) / cos(l)), l];
+  }
+
+  function obliquifySphc(latF, lonF, pole) {
+    if (pole == null)
+      // null pole indicates that this procedure should be bypassed
+      return [latF, lonF];
+
+    const lat0 = pole[0],
+      lon0 = pole[1],
+      tht0 = pole[2];
+
+    let lat1, lon1;
+    if (lat0 == halfPi) lat1 = latF;
+    else
+      lat1 = asin(
+        sin(lat0) * sin(latF) + cos(lat0) * cos(latF) * cos(lon0 - lonF)
+      ); // relative latitude
+
+    if (lat0 == halfPi)
+      // accounts for all the 0/0 errors at the poles
+      lon1 = lonF - lon0;
+    else if (lat0 == -halfPi) lon1 = lon0 - lonF - pi;
+    else {
+      lon1 =
+        acos(
+          (cos(lat0) * sin(latF) - sin(lat0) * cos(latF) * cos(lon0 - lonF)) /
+            cos(lat1)
+        ) - pi; // relative longitude
+      if (isNaN(lon1)) {
+        if (
+          (cos(lon0 - lonF) >= 0 && latF < lat0) ||
+          (cos(lon0 - lonF) < 0 && latF < -lat0)
+        )
+          lon1 = 0;
+        else lon1 = -pi;
+      } else if (sin(lonF - lon0) > 0)
+        // it's a plus-or-minus arccos.
+        lon1 = -lon1;
+    }
+    lon1 = lon1 - tht0;
+
+    return [lat1, lon1];
+  }
+
+  function obliquifyPlnr(coords, pole) {
+    if (pole == null)
+      //this indicates that you just shouldn't do this calculation
+      return coords;
+
+    let lat1 = coords[0],
+      lon1 = coords[1];
+    const lat0 = pole[0],
+      lon0 = pole[1],
+      tht0 = pole[2];
+
+    lon1 += tht0;
+    let latf = asin(sin(lat0) * sin(lat1) - cos(lat0) * cos(lon1) * cos(lat1)),
+      lonf,
+      innerFunc = sin(lat1) / cos(lat0) / cos(latf) - tan(lat0) * tan(latf);
+    if (lat0 == halfPi)
+      // accounts for special case when lat0 = pi/2
+      lonf = lon1 + lon0;
+    else if (lat0 == -halfPi)
+      // accounts for special case when lat0 = -pi/2
+      lonf = -lon1 + lon0 + pi;
+    else if (abs(innerFunc) > 1) {
+      // accounts for special case when cos(lat1) -> 0
+      if ((lon1 == 0 && lat1 < -lat0) || (lon1 != 0 && lat1 < lat0))
+        lonf = lon0 + pi;
+      else lonf = lon0;
+    } else if (sin(lon1) > 0) lonf = lon0 + acos(innerFunc);
+    else lonf = lon0 - acos(innerFunc);
+
+    let thtf = pole[2];
+
+    return [latf, lonf, thtf];
+  }
+
+  function forward(lon, lat) {
+    const width = configuration.width,
+      height = configuration.height;
+    const numSym = configuration.sphereSym; //we're about to be using this variable a lot
+    let latR = -Infinity;
+    let lonR = -Infinity;
+    let centrum = null;
+    for (const testCentrum of centrums) {
+      //iterate through the centrums to see which goes here
+      const relCoords = obliquifySphc(lat, lon, testCentrum);
+      if (relCoords[0] > latR) {
+        latR = relCoords[0];
+        lonR = relCoords[1];
+        centrum = testCentrum;
+      }
+    }
+
+    const lonR0 =
+      floor((lonR + pi / numSym) / ((2 * pi) / numSym)) * ((2 * pi) / numSym);
+
+    const rth = faceProject(lonR - lonR0, latR);
+    const r = rth[0];
+    const th = rth[1] + centrum[3] + (lonR0 * numSym) / configuration.planarSym;
+    const x0 = centrum[4];
+    const y0 = centrum[5];
+
+    let output = [r * cos(th) + x0, r * sin(th) + y0];
+    if (abs(output[0]) > width / 2 || abs(output[1]) > height / 2) {
+      output = configuration.rotateOOB(output[0], output[1], x0, y0);
+    }
+    return output;
+  }
+
+  function invert(x, y) {
+    if (isNaN(x) || isNaN(y)) return null;
+
+    if (!configuration.inBounds(x, y)) return null;
+
+    const numSym = configuration.planarSym;
+
+    let rM = +Infinity;
+    let centrum = null; //iterate to see which centrum we get
+    for (const testCentrum of centrums) {
+      const rR = hypot(x - testCentrum[4], y - testCentrum[5]);
+      if (rR < rM) {
+        //pick the centrum that minimises r
+        rM = rR;
+        centrum = testCentrum;
+      }
+    }
+    const th0 = centrum[3],
+      x0 = centrum[4],
+      y0 = centrum[5],
+      r = hypot(x - x0, y - y0),
+      th = atan2(y - y0, x - x0) - th0,
+      thBase =
+        floor((th + pi / numSym) / ((2 * pi) / numSym)) * ((2 * pi) / numSym);
+
+    let relCoords = faceInverse(r, th - thBase);
+
+    if (relCoords == null) return null;
+
+    relCoords[1] = (thBase * numSym) / configuration.sphereSym + relCoords[1];
+    let absCoords = obliquifyPlnr(relCoords, centrum);
+    return [absCoords[1], absCoords[0]];
+  }
+
+  forward.invert = invert;
+
+  return forward;
+}
+
+export function imagoBlock() {
+  var k = 0.68,
+    m = projectionMutator(imagoRaw),
+    p = m(k);
+
+  p.k = function(_) {
+    return arguments.length ? m((k = +_)) : k;
+  };
+
+  var a = -atan(1 / sqrt(2)) * degrees,
+    border = [
+      [-180 + epsilon, a + epsilon],
+      [0, 90],
+      [180 - epsilon, a + epsilon],
+      [180 - epsilon, a - epsilon],
+      [-180 + epsilon, a - epsilon],
+      [-180 + epsilon, a + epsilon]
+    ];
+
+  return p
+    .preclip(clipPolygon({
+      type: "Polygon",
+      coordinates: [border]
+      }))
+    .scale(144.04)
+    .rotate([18, -12.5, 3.5])
+    .center([0, 35.2644]);
+}
+
+function imagoWideRaw(k, shift) {
+  var imago = imagoRaw(k);
+  const height = configuration.height;
+
+  function forward(lon, lat) {
+    const p = imago(lon, lat),
+      q = [p[1], -p[0]];
+
+    if (q[1] > 0) {
+      q[0] = height - q[0];
+      q[1] *= -1;
+    }
+
+    q[0] += shift;
+    if (q[0] < 0) q[0] += height * 2;
+
+    return q;
+  }
+
+  function invert(x, y) {
+    x = (x - shift) / height;
+
+    if (x > 1.5) {
+      x -= 2;
+    }
+
+    if (x > 0.5) {
+      x = 1 - x;
+      y *= -1;
+    }
+
+    return imago.invert(-y, x * height);
+  }
+
+  forward.invert = invert;
+  return forward;
+}
+
+export default function() {
+  var k = 0.59,
+    shift = 1.16,
+    m = projectionMutator(imagoWideRaw),
+    p = m(k, shift);
+
+  p.shift = function(_) {
+    return arguments.length ? clipped(m(k, (shift = +_))) : shift;
+  };
+  p.k = function(_) {
+    return arguments.length ? clipped(m((k = +_), shift)) : k;
+  };
+
+  function clipped(p) {
+    const N = 100 + 2 * epsilon,
+      border = [],
+      e = 3e-3;
+
+    const scale = p.scale(),
+      center = p.center(),
+      translate = p.translate(),
+      rotate = p.rotate();
+    p.scale(1)
+      .center([0, 90])
+      .rotate([0, 0])
+      .translate([shift, 0]);
+    for (let i = N - epsilon; i > 0; i--) {
+      border.unshift(
+        p.invert([
+          1.5 * configuration.height - e,
+          ((configuration.width / 2) * i) / N
+        ])
+      );
+      border.push(
+        p.invert([
+          -0.5 * configuration.height + e,
+          ((configuration.width / 2) * i) / N
+        ])
+      );
+    }
+    border.push(border[0]);
+
+    return p
+      .scale(scale)
+      .center(center)
+      .translate(translate)
+      .rotate(rotate)
+      .preclip(
+        clipPolygon({
+          type: "Polygon",
+          coordinates: [border]
+        })
+      );
+  }
+
+  return clipped(p)
+    .rotate([18, -12.5, 3.5])
+    .scale(138.42)
+    .translate([480, 250])
+    .center([-139.405, 40.5844]);
+}
diff --git a/src/index.js b/src/index.js
index 5ef6aa6..b0833df 100644
--- a/src/index.js
+++ b/src/index.js
@@ -10,5 +10,6 @@ export {default as geoTetrahedralLee, leeRaw as geoLeeRaw} from "./tetrahedralLe
 export {default as geoGrayFullerRaw} from "./grayfuller";
 export {default as geoAirocean} from "./airocean";
 export {default as geoIcosahedral} from "./icosahedral";
+export {default as geoImago, imagoBlock as geoImagoBlock, imagoRaw as geoImagoRaw} from "./imago";
 export {default as geoCubic} from "./cubic";
 export {default as geoCahillKeyes, cahillKeyesRaw as geoCahillKeyesRaw} from "./cahillKeyes";
diff --git a/src/newton.js b/src/newton.js
new file mode 100644
index 0000000..8d53099
--- /dev/null
+++ b/src/newton.js
@@ -0,0 +1,16 @@
+import {abs, epsilon} from "./math";
+
+// Newton-Raphson
+// Solve f(x) = y, start from x
+export function solve(f, y, x) {
+  var steps = 100, delta, f0, f1;
+  x = x === undefined ? 0 : +x;
+  y = +y;
+  do {
+    f0 = f(x);
+    f1 = f(x + epsilon);
+    if (f0 === f1) f1 = f0 + epsilon;
+    x -= delta = (-1 * epsilon * (f0 - y)) / (f0 - f1);
+  } while (steps-- > 0 && abs(delta) > epsilon);
+  return steps < 0 ? NaN : x;
+}
diff --git a/test/compare-images b/test/compare-images
index bf0ff0b..67d64bf 100755
--- a/test/compare-images
+++ b/test/compare-images
@@ -7,6 +7,7 @@ for i in \
     cubic \
     dodecahedral \
     icosahedral \
+    imago \
     polyhedralButterfly \
     polyhedralCollignon \
     polyhedralWaterman \
diff --git a/test/invert-test.js b/test/invert-test.js
index c6b4617..c07d7d9 100644
--- a/test/invert-test.js
+++ b/test/invert-test.js
@@ -17,3 +17,16 @@ tape("inverse polyhedrals", function(test) {
   test.end();
 });
 
+tape("inverse Imago", function(test) {
+
+  [ d3.geoImago(), d3.geoImagoBlock() ]
+  .forEach(function(projection) {
+    [ [-23, 12], [10,10], [100,-45] ]
+    .forEach(function(location) {
+      projection.angle(Math.random()*360);
+      test.projectionEqual(projection, location, projection(location), 1e-5);
+    });
+  });
+  test.end();
+});
+