diff --git a/docs/harfbuzz-sections.txt b/docs/harfbuzz-sections.txt
index b475bbae3..4cf1aeb32 100644
--- a/docs/harfbuzz-sections.txt
+++ b/docs/harfbuzz-sections.txt
@@ -780,6 +780,7 @@ hb_set_t
 <FILE>hb-shape</FILE>
 hb_shape
 hb_shape_full
+hb_shape_justify
 hb_shape_list_shapers
 </SECTION>
 
diff --git a/src/gen-def.py b/src/gen-def.py
index 9b4efd4f3..6011817bc 100755
--- a/src/gen-def.py
+++ b/src/gen-def.py
@@ -19,7 +19,8 @@ symbols = sorted (re.findall (r"^hb_\w+(?= \()", "\n".join (headers_content), re
 if '--experimental-api' not in sys.argv:
 	# Move these to harfbuzz-sections.txt when got stable
 	experimental_symbols = \
-"""hb_subset_repack_or_fail
+"""hb_shape_justify
+hb_subset_repack_or_fail
 hb_subset_input_override_name_table
 """.splitlines ()
 	symbols = [x for x in symbols if x not in experimental_symbols]
diff --git a/src/hb-algs.hh b/src/hb-algs.hh
index 8c7ad434d..13587eac0 100644
--- a/src/hb-algs.hh
+++ b/src/hb-algs.hh
@@ -1344,4 +1344,62 @@ struct
 HB_FUNCOBJ (hb_dec);
 
 
+/* Adapted from kurbo implementation with extra parameters added,
+ * and finding for a particular range instead of 0.
+ *
+ * For documentation and implementation see:
+ *
+ * [ITP method]: https://en.wikipedia.org/wiki/ITP_Method
+ * [An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality]: https://dl.acm.org/doi/10.1145/3423597
+ * https://docs.rs/kurbo/0.8.1/kurbo/common/fn.solve_itp.html
+ * https://github.com/linebender/kurbo/blob/fd839c25ea0c98576c7ce5789305822675a89938/src/common.rs#L162-L248
+ */
+template <typename func_t>
+double solve_itp (func_t f,
+		  double a, double b,
+		  double epsilon,
+		  double min_y, double max_y,
+		  double &ya, double &yb, double &y)
+{
+  unsigned n1_2 = (unsigned) (hb_max (ceil (log2 ((b - a) / epsilon)) - 1.0, 0.0));
+  const unsigned n0 = 1; // Hardwired
+  const double k1 = 0.2 / (b - a); // Hardwired.
+  unsigned nmax = n0 + n1_2;
+  double scaled_epsilon = epsilon * double (1llu << nmax);
+  double _2_epsilon = 2.0 * epsilon;
+  while (b - a > _2_epsilon)
+  {
+    double x1_2 = 0.5 * (a + b);
+    double r = scaled_epsilon - 0.5 * (b - a);
+    double xf = (yb * a - ya * b) / (yb - ya);
+    double sigma = x1_2 - xf;
+    double b_a = b - a;
+    // This has k2 = 2 hardwired for efficiency.
+    double b_a_k2 = b_a * b_a;
+    double delta = k1 * b_a_k2;
+    int sigma_sign = sigma >= 0 ? +1 : -1;
+    double xt = delta <= fabs (x1_2 - xf) ? xf + delta * sigma_sign : x1_2;
+    double xitp = fabs (xt - x1_2) <= r ? xt : x1_2 - r * sigma_sign;
+    double yitp = f (xitp);
+    if (yitp > max_y)
+    {
+      b = xitp;
+      yb = yitp;
+    }
+    else if (yitp < min_y)
+    {
+      a = xitp;
+      ya = yitp;
+    }
+    else
+    {
+      y = yitp;
+      return xitp;
+    }
+    scaled_epsilon *= 0.5;
+  }
+  return 0.5 * (a + b);
+}
+
+
 #endif /* HB_ALGS_HH */
diff --git a/src/hb-debug.hh b/src/hb-debug.hh
index efab37464..0ac4515fa 100644
--- a/src/hb-debug.hh
+++ b/src/hb-debug.hh
@@ -373,6 +373,10 @@ struct hb_no_trace_t {
 #define HB_DEBUG_FT (HB_DEBUG+0)
 #endif
 
+#ifndef HB_DEBUG_JUSTIFY
+#define HB_DEBUG_JUSTIFY (HB_DEBUG+0)
+#endif
+
 #ifndef HB_DEBUG_OBJECT
 #define HB_DEBUG_OBJECT (HB_DEBUG+0)
 #endif
diff --git a/src/hb-shape.cc b/src/hb-shape.cc
index 7b5bf2c5e..13bda36f3 100644
--- a/src/hb-shape.cc
+++ b/src/hb-shape.cc
@@ -26,6 +26,8 @@
  * Google Author(s): Behdad Esfahbod
  */
 
+#define HB_DEBUG_JUSTIFY 1
+
 #include "hb.hh"
 
 #include "hb-shaper.hh"
@@ -196,4 +198,212 @@ hb_shape (hb_font_t           *font,
 }
 
 
+#ifdef HB_EXPERIMENTAL_API
+
+static float
+buffer_advance (hb_buffer_t *buffer)
+{
+  float a = 0;
+  auto *pos = buffer->pos;
+  unsigned count = buffer->len;
+  if (HB_DIRECTION_IS_HORIZONTAL (buffer->props.direction))
+    for (unsigned i = 0; i < count; i++)
+      a += pos[i].x_advance;
+  else
+    for (unsigned i = 0; i < count; i++)
+      a += pos[i].y_advance;
+  return a;
+}
+
+static void
+reset_buffer (hb_buffer_t *buffer,
+	      hb_array_t<const hb_glyph_info_t> text)
+{
+  assert (buffer->ensure (text.length));
+  buffer->have_positions = false;
+  buffer->len = text.length;
+  memcpy (buffer->info, text.arrayZ, text.length * sizeof (buffer->info[0]));
+  hb_buffer_set_content_type (buffer, HB_BUFFER_CONTENT_TYPE_UNICODE);
+}
+
+/**
+ * hb_shape_justify:
+ * @font: a mutable #hb_font_t to use for shaping
+ * @buffer: an #hb_buffer_t to shape
+ * @features: (array length=num_features) (nullable): an array of user
+ *    specified #hb_feature_t or `NULL`
+ * @num_features: the length of @features array
+ * @shaper_list: (array zero-terminated=1) (nullable): a `NULL`-terminated
+ *    array of shapers to use or `NULL`
+ * @min_target_advance: Minimum advance width/height to aim for.
+ * @max_target_advance: Maximum advance width/height to aim for.
+ * @advance: (inout): Input/output advance width/height of the buffer.
+ * @var_tag: (out): Variation-axis tag used for justification.
+ * @var_value: (out): Variation-axis value used to reach target justification.
+ *
+ * See hb_shape_full() for basic details. If @shaper_list is not `NULL`, the specified
+ * shapers will be used in the given order, otherwise the default shapers list
+ * will be used.
+ *
+ * In addition, justify the shaping results such that the shaping results reach
+ * the target advance width/height, depending on the buffer direction.
+ *
+ * This API is currently experimental and will probably change in the future.
+ *
+ * Return value: false if all shapers failed, true otherwise
+ *
+ * XSince: EXPERIMENTAL
+ **/
+hb_bool_t
+hb_shape_justify (hb_font_t          *font,
+		  hb_buffer_t        *buffer,
+		  const hb_feature_t *features,
+		  unsigned int        num_features,
+		  const char * const *shaper_list,
+		  float               min_target_advance,
+		  float               max_target_advance,
+		  float              *advance, /* IN/OUT */
+		  hb_tag_t           *var_tag, /* OUT */
+		  float              *var_value /* OUT */)
+{
+  // TODO Negative font scales?
+
+  if (min_target_advance <= *advance && *advance <= max_target_advance)
+    return hb_shape_full (font, buffer,
+			  features, num_features,
+			  shaper_list);
+
+  hb_face_t *face = font->face;
+
+  hb_tag_t tag = HB_TAG_NONE;
+  hb_ot_var_axis_info_t axis_info;
+
+  hb_tag_t tags[] =
+  {
+    HB_TAG ('j','s','t','f'),
+    HB_TAG ('w','d','t','h'),
+  };
+  for (unsigned i = 0; i < ARRAY_LENGTH (tags); i++)
+    if (hb_ot_var_find_axis_info (face, tags[i], &axis_info))
+    {
+      tag = *var_tag = tags[i];
+      break;
+    }
+
+  if (!tag)
+  {
+    if (hb_shape_full (font, buffer,
+		       features, num_features,
+		       shaper_list))
+    {
+      *advance = buffer_advance (buffer);
+      return true;
+    }
+    else
+      return false;
+  }
+
+  unsigned text_len = buffer->len;
+  auto *text_info = (hb_glyph_info_t *) hb_malloc (text_len * sizeof (buffer->info[0]));
+  if (unlikely (text_len && !text_info))
+    return false;
+  hb_memcpy (text_info, buffer->info, text_len * sizeof (buffer->info[0]));
+  auto text = hb_array<const hb_glyph_info_t> (text_info, text_len);
+
+  if (!*advance)
+  {
+    hb_font_set_variation (font, tag, axis_info.default_value);
+    if (!hb_shape_full (font, buffer,
+			features, num_features,
+			shaper_list))
+      return false;
+    *advance = buffer_advance (buffer);
+  }
+
+  if (min_target_advance <= *advance && *advance <= max_target_advance)
+    return true;
+
+  double a, b, ya, yb;
+  if (*advance < min_target_advance)
+  {
+    ya = *advance;
+    a = axis_info.default_value;
+    b = axis_info.max_value;
+
+    hb_font_set_variation (font, tag, (float) b);
+    reset_buffer (buffer, text);
+    if (!hb_shape_full (font, buffer,
+			features, num_features,
+			shaper_list))
+      return false;
+    yb = buffer_advance (buffer);
+    if (yb <= 0)
+    {
+      *advance = (float) yb;
+      return true;
+    }
+  }
+  else
+  {
+    yb = *advance;
+    a = axis_info.min_value;
+    b = axis_info.default_value;
+
+    hb_font_set_variation (font, tag, (float) a);
+    reset_buffer (buffer, text);
+    if (!hb_shape_full (font, buffer,
+			features, num_features,
+			shaper_list))
+      return false;
+    ya = buffer_advance (buffer);
+    if (ya >= 0)
+    {
+      *advance = (float) ya;
+      return true;
+    }
+  }
+
+  double epsilon = (b - a) / (1<<14);
+  bool failed = false;
+
+  auto f = [&] (double x)
+  {
+    hb_font_set_variation (font, tag, (float) x);
+    reset_buffer (buffer, text);
+    if (unlikely (!hb_shape_full (font, buffer,
+				  features, num_features,
+				  shaper_list)))
+    {
+      failed = true;
+      return (double) min_target_advance;
+    }
+
+    double w = buffer_advance (buffer);
+    DEBUG_MSG (JUSTIFY, nullptr, "Trying '%c%c%c%c' axis parameter %f. Advance %g. Target: min %g max %g",
+	       HB_UNTAG (tag), x, w,
+	       (double) min_target_advance, (double) max_target_advance);
+    return w;
+  };
+
+  double y = 0;
+  double itp = solve_itp (f,
+			  a, b,
+			  epsilon,
+			  min_target_advance, max_target_advance,
+			  ya, yb, y);
+
+  hb_free (text_info);
+
+  if (failed)
+    return false;
+
+  *var_value = (float) itp;
+  *advance = (float) y;
+
+  return true;
+}
+
+#endif
+
+
 #endif
diff --git a/src/hb-shape.h b/src/hb-shape.h
index 922f8c011..d4d4fdfd2 100644
--- a/src/hb-shape.h
+++ b/src/hb-shape.h
@@ -53,6 +53,18 @@ hb_shape_full (hb_font_t          *font,
 	       unsigned int        num_features,
 	       const char * const *shaper_list);
 
+HB_EXTERN hb_bool_t
+hb_shape_justify (hb_font_t          *font,
+		  hb_buffer_t        *buffer,
+		  const hb_feature_t *features,
+		  unsigned int        num_features,
+		  const char * const *shaper_list,
+		  float               min_target_advance,
+		  float               max_target_advance,
+		  float              *advance, /* IN/OUT */
+		  hb_tag_t           *var_tag, /* OUT */
+		  float              *var_value /* OUT */);
+
 HB_EXTERN const char **
 hb_shape_list_shapers (void);
 
diff --git a/util/shape-options.hh b/util/shape-options.hh
index b4d9b2c00..b22b95883 100644
--- a/util/shape-options.hh
+++ b/util/shape-options.hh
@@ -161,12 +161,52 @@ struct shape_options_t
     }
     else
     {
-      if (!hb_shape_full (font, buffer, features, num_features, shapers))
+      if (advance <= 0)
       {
-	if (error)
-	  *error = "Shaping failed.";
-	goto fail;
+	if (!hb_shape_full (font, buffer, features, num_features, shapers))
+	{
+	  if (error)
+	    *error = "Shaping failed.";
+	  goto fail;
+	}
+
+	if (advance < 0)
+	{
+	  float unit = (1 << SUBPIXEL_BITS);
+
+	  /* Calculate buffer advance */
+	  float w = 0;
+	  unsigned count = 0;
+	  hb_glyph_position_t *pos = hb_buffer_get_glyph_positions (buffer, &count);
+	  if (HB_DIRECTION_IS_HORIZONTAL (hb_buffer_get_direction (buffer)))
+	    for (unsigned i = 0; i < count; i++)
+	      w += pos[i].x_advance;
+	  else
+	    for (unsigned i = 0; i < count; i++)
+	      w += pos[i].y_advance;
+
+	  printf ("Default size: %u\n", (unsigned) roundf (w / unit));
+	  exit (0);
+	}
       }
+#ifdef HB_EXPERIMENTAL_API
+      else
+      {
+        float unit = (1 << SUBPIXEL_BITS);
+        float target_advance = advance * unit;
+	float w = 0;
+	hb_tag_t var_tag;
+	float var_value;
+	if (!hb_shape_justify (font, buffer, features, num_features, shapers,
+			       target_advance - unit * 0.5f, target_advance + unit * 0.5f,
+			       &w, &var_tag, &var_value))
+	{
+	  if (error)
+	    *error = "Shaping failed.";
+	  goto fail;
+	}
+      }
+#endif
     }
 
     if (normalize_glyphs)
@@ -202,6 +242,7 @@ struct shape_options_t
   hb_feature_t *features = nullptr;
   unsigned int num_features = 0;
   char **shapers = nullptr;
+  signed advance = 0;
   hb_bool_t utf8_clusters = false;
   hb_codepoint_t invisible_glyph = 0;
   hb_codepoint_t not_found_glyph = 0;
@@ -328,6 +369,10 @@ shape_options_t::add_options (option_parser_t *parser)
     {"script",		0, 0, G_OPTION_ARG_STRING,	&this->script,			"Set text script (default: auto)",	"ISO-15924 tag"},
     {"bot",		0, 0, G_OPTION_ARG_NONE,	&this->bot,			"Treat text as beginning-of-paragraph",	nullptr},
     {"eot",		0, 0, G_OPTION_ARG_NONE,	&this->eot,			"Treat text as end-of-paragraph",	nullptr},
+#ifdef HB_EXPERIMENTAL_API
+    {"justify-to",	0, 0,
+			      G_OPTION_ARG_INT,		&this->advance,			"Target size to justify to",		"SIZE, or -1"},
+#endif
     {"preserve-default-ignorables",0, 0, G_OPTION_ARG_NONE,	&this->preserve_default_ignorables,	"Preserve Default-Ignorable characters",	nullptr},
     {"remove-default-ignorables",0, 0, G_OPTION_ARG_NONE,	&this->remove_default_ignorables,	"Remove Default-Ignorable characters",	nullptr},
     {"invisible-glyph",	0, 0, G_OPTION_ARG_INT,		&this->invisible_glyph,		"Glyph value to replace Default-Ignorables with",	nullptr},
@@ -339,7 +384,7 @@ shape_options_t::add_options (option_parser_t *parser)
     {"safe-to-insert-tatweel",0, 0, G_OPTION_ARG_NONE,	&this->safe_to_insert_tatweel,	"Produce safe-to-insert-tatweel glyph flag",	nullptr},
     {"glyphs",		0, 0, G_OPTION_ARG_NONE,	&this->glyphs,			"Interpret input as glyph string",	nullptr},
     {"verify",		0, 0, G_OPTION_ARG_NONE,	&this->verify,			"Perform sanity checks on shaping results",	nullptr},
-    {"num-iterations", 'n', G_OPTION_FLAG_IN_MAIN,
+    {"num-iterations",	'n',G_OPTION_FLAG_IN_MAIN,
 			      G_OPTION_ARG_INT,		&this->num_iterations,		"Run shaper N times (default: 1)",	"N"},
     {nullptr}
   };