Reland: Optimize the transformRect and transformPoint methods in matrix_utils. (#36396) (#37275)

Primarily these methods no longer allocate any objects other than their
return values.

Additionally, the math in the methods is reduced compared to the general
case math based on known input conditions.

Modified to no longer generate infinite values in some finite cases.

packages/flutter/lib/src/painting/matrix_utils.dart

@@ -2,7 +2,6 @@
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.
 
-import 'dart:math' as math;
 import 'dart:typed_data';
 
 import 'package:flutter/foundation.dart';
@@ -124,9 +123,70 @@ class MatrixUtils {
   /// This function assumes the given point has a z-coordinate of 0.0. The
   /// z-coordinate of the result is ignored.
   static Offset transformPoint(Matrix4 transform, Offset point) {
-    final Vector3 position3 = Vector3(point.dx, point.dy, 0.0);
-    final Vector3 transformed3 = transform.perspectiveTransform(position3);
-    return Offset(transformed3.x, transformed3.y);
+    final Float64List storage = transform.storage;
+    final double x = point.dx;
+    final double y = point.dy;
+
+    // Directly simulate the transform of the vector (x, y, 0, 1),
+    // dropping the resulting Z coordinate, and normalizing only
+    // if needed.
+
+    final double rx = storage[0] * x + storage[4] * y + storage[12];
+    final double ry = storage[1] * x + storage[5] * y + storage[13];
+    final double rw = storage[3] * x + storage[7] * y + storage[15];
+    if (rw == 1.0) {
+      return Offset(rx, ry);
+    } else {
+      return Offset(rx / rw, ry / rw);
+    }
+  }
+
+  /// Returns a rect that bounds the result of applying the given matrix as a
+  /// perspective transform to the given rect.
+  ///
+  /// This version of the operation is slower than the regular transformRect
+  /// method, but it avoids creating infinite values from large finite values
+  /// if it can.
+  static Rect _safeTransformRect(Matrix4 transform, Rect rect) {
+    final Float64List storage = transform.storage;
+    final bool isAffine = storage[3] == 0.0 &&
+        storage[7] == 0.0 &&
+        storage[15] == 1.0;
+
+    _minMax ??= Float64List(4);
+
+    _accumulate(storage, rect.left,  rect.top,    true,  isAffine);
+    _accumulate(storage, rect.right, rect.top,    false, isAffine);
+    _accumulate(storage, rect.left,  rect.bottom, false, isAffine);
+    _accumulate(storage, rect.right, rect.bottom, false, isAffine);
+
+    return Rect.fromLTRB(_minMax[0], _minMax[1], _minMax[2], _minMax[3]);
+  }
+
+  static Float64List _minMax;
+  static void _accumulate(Float64List m, double x, double y,
+      bool first, bool isAffine)
+  {
+    final double w = isAffine ? 1.0 : 1.0 / (m[3] * x + m[7] * y + m[15]);
+    final double tx = (m[0] * x + m[4] * y + m[12]) * w;
+    final double ty = (m[1] * x + m[5] * y + m[13]) * w;
+    if (first) {
+      _minMax[0] = _minMax[2] = tx;
+      _minMax[1] = _minMax[3] = ty;
+    } else {
+      if (tx < _minMax[0]) {
+        _minMax[0] = tx;
+      }
+      if (ty < _minMax[1]) {
+        _minMax[1] = ty;
+      }
+      if (tx > _minMax[2]) {
+        _minMax[2] = tx;
+      }
+      if (ty > _minMax[3]) {
+        _minMax[3] = ty;
+      }
+    }
   }
 
   /// Returns a rect that bounds the result of applying the given matrix as a
@@ -136,23 +196,232 @@ class MatrixUtils {
   /// The transformed rect is then projected back into the plane with z equals
   /// 0.0 before computing its bounding rect.
   static Rect transformRect(Matrix4 transform, Rect rect) {
-    final Offset point1 = transformPoint(transform, rect.topLeft);
-    final Offset point2 = transformPoint(transform, rect.topRight);
-    final Offset point3 = transformPoint(transform, rect.bottomLeft);
-    final Offset point4 = transformPoint(transform, rect.bottomRight);
+    final Float64List storage = transform.storage;
+    final double x = rect.left;
+    final double y = rect.top;
+    final double w = rect.right - x;
+    final double h = rect.bottom - y;
+
+    // We want to avoid turning a finite rect into an infinite one if we can.
+    if (!w.isFinite || !h.isFinite) {
+      return _safeTransformRect(transform, rect);
+    }
+
+    // Transforming the 4 corners of a rectangle the straightforward way
+    // incurs the cost of transforming 4 points using vector math which
+    // involves 48 multiplications and 48 adds and then normalizing
+    // the points using 4 inversions of the homogeneous weight factor
+    // and then 12 multiplies. Once we have transformed all of the points
+    // we then need to turn them into a bounding box using 4 min/max
+    // operations each on 4 values yielding 12 total comparisons.
+    //
+    // On top of all of those operations, using the vector_math package to
+    // do the work for us involves allocating several objects in order to
+    // communicate the values back and forth - 4 allocating getters to extract
+    // the [Offset] objects for the corners of the [Rect], 4 conversions to
+    // a [Vector3] to use [Matrix4.perspectiveTransform()], and then 4 new
+    // [Offset] objects allocated to hold those results, yielding 8 [Offset]
+    // and 4 [Vector3] object allocations per rectangle transformed.
+    //
+    // But the math we really need to get our answer is actually much less
+    // than that.
+    //
+    // First, consider that a full point transform using the vector math
+    // package involves expanding it out into a vector3 with a Z coordinate
+    // of 0.0 and then performing 3 multiplies and 3 adds per coordinate:
+    // ```
+    // xt = x*m00 + y*m10 + z*m20 + m30;
+    // yt = x*m01 + y*m11 + z*m21 + m31;
+    // zt = x*m02 + y*m12 + z*m22 + m32;
+    // wt = x*m03 + y*m13 + z*m23 + m33;
+    // ```
+    // Immediately we see that we can get rid of the 3rd column of multiplies
+    // since we know that Z=0.0. We can also get rid of the 3rd row because
+    // we ignore the resulting Z coordinate. Finally we can get rid of the
+    // last row if we don't have a perspective transform since we can verify
+    // that the results are 1.0 for all points.  This gets us down to 16
+    // multiplies and 16 adds in the non-perspective case and 24 of each for
+    // the perspective case. (Plus the 12 comparisons to turn them back into
+    // a bounding box.)
+    //
+    // But we can do even better than that.
+    //
+    // Under optimal conditions of no perspective transformation,
+    // which is actually a very common condition, we can transform
+    // a rectangle in as little as 3 operations:
+    //
+    // (rx,ry) = transform of upper left corner of rectangle
+    // (wx,wy) = delta transform of the (w, 0) width relative vector
+    // (hx,hy) = delta transform of the (0, h) height relative vector
+    //
+    // A delta transform is a transform of all elements of the matrix except
+    // for the translation components. The translation components are added
+    // in at the end of each transform computation so they represent a
+    // constant offset for each point transformed. A delta transform of
+    // a horizontal or vertical vector involves a single multiplication due
+    // to the fact that it only has one non-zero coordinate and no addition
+    // of the translation component.
+    //
+    // In the absence of a perspective transform, the transformed
+    // rectangle will be mapped into a parallelogram with corners at:
+    // corner1 = (rx, ry)
+    // corner2 = corner1 + dTransformed width vector = (rx+wx, ry+wy)
+    // corner3 = corner1 + dTransformed height vector = (rx+hx, ry+hy)
+    // corner4 = corner1 + both dTransformed vectors = (rx+wx+hx, ry+wy+hy)
+    // In all, this method of transforming the rectangle requires only
+    // 8 multiplies and 12 additions (which we can reduce to 8 additions if
+    // we only need a bounding box, see below).
+    //
+    // In the presence of a perspective transform, the above conditions
+    // continue to hold with respect to the non-normalized coordinates so
+    // we can still save a lot of multiplications by computing the 4
+    // non-normalized coordinates using relative additions before we normalize
+    // them and they lose their "pseudo-parallelogram" relationships.  We still
+    // have to do the normalization divisions and min/max all 4 points to
+    // get the resulting transformed bounding box, but we save a lot of
+    // calculations over blindly transforming all 4 coordinates independently.
+    // In all, we need 12 multiplies and 22 additions to construct the
+    // non-normalized vectors and then 8 divisions (or 4 inversions and 8
+    // multiplies) for normalization (plus the standard set of 12 comparisons
+    // for the min/max bounds operations).
+    //
+    // Back to the non-perspective case, the optimization that lets us get
+    // away with fewer additions if we only need a bounding box comes from
+    // analyzing the impact of the relative vectors on expanding the
+    // bounding box of the parallelogram. First, the bounding box always
+    // contains the transformed upper-left corner of the rectangle. Next,
+    // each relative vector either pushes on the left or right side of the
+    // bounding box and also either the top or bottom side, depending on
+    // whether it is positive or negative. Finally, you can consider the
+    // impact of each vector on the bounding box independently. If, say,
+    // wx and hx have the same sign, then the limiting point in the bounding
+    // box will be the one that involves adding both of them to the origin
+    // point. If they have opposite signs, then one will push one wall one
+    // way and the other will push the opposite wall the other way and when
+    // you combine both of them, the resulting "opposite corner" will
+    // actually be between the limits they established by pushing the walls
+    // away from each other, as below:
+    // ```
+    //         +---------(originx,originy)--------------+
+    //         |            -----^----                  |
+    //         |       -----          ----              |
+    //         |  -----                   ----          |
+    // (+hx,+hy)<                             ----      |
+    //         |  ----                            ----  |
+    //         |      ----                             >(+wx,+wy)
+    //         |          ----                   -----  |
+    //         |              ----          -----       |
+    //         |                  ---- -----            |
+    //         |                      v                 |
+    //         +---------------(+wx+hx,+wy+hy)----------+
+    // ```
+    // In this diagram, consider that:
+    // ```
+    // wx would be a positive number
+    // hx would be a negative number
+    // wy and hy would both be positive numbers
+    // ```
+    // As a result, wx pushes out the right wall, hx pushes out the left wall,
+    // and both wy and hy push down the bottom wall of the bounding box. The
+    // wx,hx pair (of opposite signs) worked on opposite walls and the final
+    // opposite corner had an X coordinate between the limits they established.
+    // The wy,hy pair (of the same sign) both worked together to push the
+    // bottom wall down by their sum.
+    //
+    // This relationship allows us to simply start with the point computed by
+    // transforming the upper left corner of the rectangle, and then
+    // conditionally adding wx, wy, hx, and hy to either the left or top
+    // or right or bottom of the bounding box independently depending on sign.
+    // In that case we only need 4 comparisons and 4 additions total to
+    // compute the bounding box, combined with the 8 multiplications and
+    // 4 additions to compute the transformed point and relative vectors
+    // for a total of 8 multiplies, 8 adds, and 4 comparisons.
+    //
+    // An astute observer will note that we do need to do 2 subtractions at
+    // the top of the method to compute the width and height.  Add those to
+    // all of the relative solutions listed above.  The test for perspective
+    // also adds 3 compares to the affine case and up to 3 compares to the
+    // perspective case (depending on which test fails, the rest are omitted).
+    //
+    // The final tally:
+    // basic method          = 60 mul + 48 add + 12 compare
+    // optimized perspective = 12 mul + 22 add + 15 compare + 2 sub
+    // optimized affine      =  8 mul +  8 add +  7 compare + 2 sub
+    //
+    // Since compares are essentially subtractions and subtractions are
+    // the same cost as adds, we end up with:
+    // basic method          = 60 mul + 60 add/sub/compare
+    // optimized perspective = 12 mul + 39 add/sub/compare
+    // optimized affine      =  8 mul + 17 add/sub/compare
+
+    final double wx = storage[0] * w;
+    final double hx = storage[4] * h;
+    final double rx = storage[0] * x + storage[4] * y + storage[12];
+
+    final double wy = storage[1] * w;
+    final double hy = storage[5] * h;
+    final double ry = storage[1] * x + storage[5] * y + storage[13];
+
+    if (storage[3] == 0.0 && storage[7] == 0.0 && storage[15] == 1.0) {
+      double left  = rx;
+      double right = rx;
+      if (wx < 0) {
+        left  += wx;
+      } else {
+        right += wx;
+      }
+      if (hx < 0) {
+        left  += hx;
+      } else {
+        right += hx;
+      }
+
+      double top    = ry;
+      double bottom = ry;
+      if (wy < 0) {
+        top    += wy;
+      } else {
+        bottom += wy;
+      }
+      if (hy < 0) {
+        top    += hy;
+      } else {
+        bottom += hy;
+      }
+
+      return Rect.fromLTRB(left, top, right, bottom);
+    } else {
+      final double ww = storage[3] * w;
+      final double hw = storage[7] * h;
+      final double rw = storage[3] * x + storage[7] * y + storage[15];
+
+      final double ulx =  rx            /  rw;
+      final double uly =  ry            /  rw;
+      final double urx = (rx + wx)      / (rw + ww);
+      final double ury = (ry + wy)      / (rw + ww);
+      final double llx = (rx      + hx) / (rw      + hw);
+      final double lly = (ry      + hy) / (rw      + hw);
+      final double lrx = (rx + wx + hx) / (rw + ww + hw);
+      final double lry = (ry + wy + hy) / (rw + ww + hw);
+
       return Rect.fromLTRB(
-        _min4(point1.dx, point2.dx, point3.dx, point4.dx),
-        _min4(point1.dy, point2.dy, point3.dy, point4.dy),
-        _max4(point1.dx, point2.dx, point3.dx, point4.dx),
-        _max4(point1.dy, point2.dy, point3.dy, point4.dy),
+        _min4(ulx, urx, llx, lrx),
+        _min4(uly, ury, lly, lry),
+        _max4(ulx, urx, llx, lrx),
+        _max4(uly, ury, lly, lry),
       );
     }
+  }
 
   static double _min4(double a, double b, double c, double d) {
-    return math.min(a, math.min(b, math.min(c, d)));
+    final double e = (a < b) ? a : b;
+    final double f = (c < d) ? c : d;
+    return (e < f) ? e : f;
   }
   static double _max4(double a, double b, double c, double d) {
-    return math.max(a, math.max(b, math.max(c, d)));
+    final double e = (a > b) ? a : b;
+    final double f = (c > d) ? c : d;
+    return (e > f) ? e : f;
   }
 
   /// Returns a rect that bounds the result of applying the inverse of the given

packages/flutter/test/painting/matrix_utils_test.dart

@@ -9,7 +9,7 @@ import 'package:flutter/painting.dart';
 import 'package:vector_math/vector_math_64.dart';
 
 void main() {
-  test('MatrixUtils.transformRect handles very small values', () {
+  test('MatrixUtils.transformRect handles very large finite values', () {
     const Rect evilRect = Rect.fromLTRB(0.0, -1.7976931348623157e+308, 800.0, 1.7976931348623157e+308);
     final Matrix4 transform = Matrix4.identity()..translate(10.0, 0.0);
     final Rect transformedRect = MatrixUtils.transformRect(transform, evilRect);
@@ -128,4 +128,78 @@ void main() {
       forcedOffset,
     );
   });
+
+  test('transformRect with no perspective (w = 1)', () {
+    const Rect rectangle20x20 = Rect.fromLTRB(10, 20, 30, 40);
+
+    // Identity
+    expect(
+      MatrixUtils.transformRect(Matrix4.identity(), rectangle20x20),
+      rectangle20x20,
+    );
+
+    // 2D Scaling
+    expect(
+      MatrixUtils.transformRect(Matrix4.diagonal3Values(2, 2, 2), rectangle20x20),
+      const Rect.fromLTRB(20, 40, 60, 80),
+    );
+
+    // Rotation
+    expect(
+      MatrixUtils.transformRect(Matrix4.rotationZ(pi / 2.0), rectangle20x20),
+      within<Rect>(distance: 0.00001, from: const Rect.fromLTRB(-40.0, 10.0, -20.0, 30.0)),
+    );
+  });
+
+  test('transformRect with perspective (w != 1)', () {
+    final Matrix4 transform = MatrixUtils.createCylindricalProjectionTransform(
+      radius: 10.0,
+      angle: pi / 8.0,
+      perspective: 0.3,
+    );
+
+    for (int i = 1; i < 10000; i++) {
+      final Rect rect = Rect.fromLTRB(11.0 * i, 12.0 * i, 15.0 * i, 18.0 * i);
+      final Rect golden = _vectorWiseTransformRect(transform, rect);
+      expect(
+        MatrixUtils.transformRect(transform, rect),
+        within<Rect>(distance: 0.00001, from: golden),
+      );
+    }
+  });
+}
+
+// Produces the same computation as `MatrixUtils.transformPoint` but it uses
+// the built-in perspective transform methods in the Matrix4 class as a
+// golden implementation of the optimized `MatrixUtils.transformPoint`
+// to make sure optimizations do not contain bugs.
+Offset _transformPoint(Matrix4 transform, Offset point) {
+  final Vector3 position3 = Vector3(point.dx, point.dy, 0.0);
+  final Vector3 transformed3 = transform.perspectiveTransform(position3);
+  return Offset(transformed3.x, transformed3.y);
+}
+
+// Produces the same computation as `MatrixUtils.transformRect` but it does this
+// one point at a time. This function is used as the golden implementation of the
+// optimized `MatrixUtils.transformRect` to make sure optimizations do not contain
+// bugs.
+Rect _vectorWiseTransformRect(Matrix4 transform, Rect rect) {
+  final Offset point1 = _transformPoint(transform, rect.topLeft);
+  final Offset point2 = _transformPoint(transform, rect.topRight);
+  final Offset point3 = _transformPoint(transform, rect.bottomLeft);
+  final Offset point4 = _transformPoint(transform, rect.bottomRight);
+  return Rect.fromLTRB(
+    _min4(point1.dx, point2.dx, point3.dx, point4.dx),
+    _min4(point1.dy, point2.dy, point3.dy, point4.dy),
+    _max4(point1.dx, point2.dx, point3.dx, point4.dx),
+    _max4(point1.dy, point2.dy, point3.dy, point4.dy),
+  );
+}
+
+double _min4(double a, double b, double c, double d) {
+  return min(a, min(b, min(c, d)));
+}
+
+double _max4(double a, double b, double c, double d) {
+  return max(a, max(b, max(c, d)));
 }