[base] Implement fast vector normalization. The function uses Newton's iterations instead of dividing vector components by its length, which needs a square root. This is, literally, a bit less accurate but a lot faster. * src/base/ftcalc.c (FT_Vector_NormLen): New function.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
diff --git a/ChangeLog b/ChangeLog
index a20be8f..af25864 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+2015-06-29 Alexei Podtelezhnikov <apodtele@gmail.com>
+
+ [base] Implement fast vector normalization.
+
+ The function uses Newton's iterations instead of dividing vector
+ components by its length, which needs a square root. This is,
+ literally, a bit less accurate but a lot faster.
+
+ * src/base/ftcalc.c (FT_Vector_NormLen): New function.
+
2015-06-28 Werner Lemberg <wl@gnu.org>
* CMakeLists.txt: Always create `ftconfig.h'.
diff --git a/include/freetype/internal/ftcalc.h b/include/freetype/internal/ftcalc.h
index 75752c3..d4d4337 100644
--- a/include/freetype/internal/ftcalc.h
+++ b/include/freetype/internal/ftcalc.h
@@ -300,6 +300,18 @@ FT_BEGIN_HEADER
/*
+ * This function normalizes a vector and returns its original length.
+ * The normalized vector is a 16.16 fixed-point unit vector with length
+ * close to 0x10000. The accuracy of the returned length is limited to
+ * 16 bits also. The function utilizes quick inverse square root
+ * aproximation without divisions and square roots relying on Newton's
+ * iterations instead.
+ */
+ FT_BASE( FT_UInt32 )
+ FT_Vector_NormLen( FT_Vector* vector );
+
+
+ /*
* Return -1, 0, or +1, depending on the orientation of a given corner.
* We use the Cartesian coordinate system, with positive vertical values
* going upwards. The function returns +1 if the corner turns to the
diff --git a/src/base/ftcalc.c b/src/base/ftcalc.c
index dca0e1d..f3595e7 100644
--- a/src/base/ftcalc.c
+++ b/src/base/ftcalc.c
@@ -785,6 +785,88 @@
}
+ /* documentation is in ftcalc.h */
+
+ FT_BASE_DEF( FT_UInt32 )
+ FT_Vector_NormLen( FT_Vector* vector )
+ {
+ FT_Int32 x = vector->x;
+ FT_Int32 y = vector->y;
+ FT_Int32 b, z;
+ FT_UInt32 u, v, l;
+ FT_Int sx = 1, sy = 1, shift;
+
+
+ FT_MOVE_SIGN( x, sx );
+ FT_MOVE_SIGN( y, sy );
+
+ /* trivial cases */
+ if ( x == 0 )
+ {
+ if ( y > 0 )
+ vector->y = sy * 0x10000;
+ return y;
+ }
+ else if ( y == 0 )
+ {
+ if ( x > 0 )
+ vector->x = sx * 0x10000;
+ return x;
+ }
+
+ /* estimate length and prenormalize */
+ l = x > y ? (FT_UInt32)x + ( y >> 1 )
+ : (FT_UInt32)y + ( x >> 1 );
+
+ shift = 31 - FT_MSB( l );
+ shift -= 15 + ( l >= 0xAAAAAAAAUL >> shift );
+
+ if ( shift > 0 )
+ {
+ x <<= shift;
+ y <<= shift;
+
+ /* reestimate length for tiny vectors */
+ l = x > y ? (FT_UInt32)x + ( y >> 1 )
+ : (FT_UInt32)y + ( x >> 1 );
+ }
+ else
+ {
+ x >>= -shift;
+ y >>= -shift;
+ l >>= -shift;
+ }
+
+ /* lower linear approximation for reciprocal length minus one */
+ b = 0x10000 - (FT_Int32)l;
+
+ /* Newton's iterations */
+ do
+ {
+ u = (FT_UInt32)( x + ( x * b >> 16 ) );
+ v = (FT_UInt32)( y + ( y * b >> 16 ) );
+
+ /* converting to signed gives difference with 2^32 */
+ z = -(FT_Int32)( u * u + v * v ) / 0x200;
+ z = z * ( ( 0x10000 + b ) >> 8 ) / 0x10000;
+
+ b += z;
+ } while ( z > 0 );
+
+ vector->x = sx < 0 ? -(FT_Pos)u : (FT_Pos)u;
+ vector->y = sy < 0 ? -(FT_Pos)v : (FT_Pos)v;
+
+ /* true length, again taking advantage of signed difference with 2^32 */
+ l = 0x10000 + (FT_Int32)( u * x + v * y ) / 0x10000;
+ if ( shift > 0 )
+ l = ( l + ( 1 << ( shift - 1 ) ) ) >> shift;
+ else
+ l <<= -shift;
+
+ return l;
+ }
+
+
#if 0
/* documentation is in ftcalc.h */