diff --git a/src/base/cg.rs b/src/base/cg.rs
index 6a7db4bd..935388af 100644
--- a/src/base/cg.rs
+++ b/src/base/cg.rs
@@ -11,11 +11,12 @@ use crate::base::allocator::Allocator;
 use crate::base::dimension::{DimName, DimNameDiff, DimNameSub, U1};
 use crate::base::storage::{Storage, StorageMut};
 use crate::base::{
-    DefaultAllocator, Matrix3, Matrix4, MatrixN, Scalar, SquareMatrix, Unit, Vector, Vector3,
-    VectorN,
+    DefaultAllocator, Matrix3, Matrix4, MatrixN, Scalar, SquareMatrix, Unit, Vector, Vector2,
+    Vector3, VectorN,
 };
 use crate::geometry::{
-    Isometry, IsometryMatrix3, Orthographic3, Perspective3, Point, Point3, Rotation2, Rotation3,
+    Isometry, IsometryMatrix3, Orthographic3, Perspective3, Point, Point2, Point3, Rotation2,
+    Rotation3,
 };
 
 use simba::scalar::{ClosedAdd, ClosedMul, RealField};
@@ -70,6 +71,26 @@ impl<N: RealField> Matrix3<N> {
     pub fn new_rotation(angle: N) -> Self {
         Rotation2::new(angle).to_homogeneous()
     }
+
+    /// Creates a new homogeneous matrix that applies a scaling factor for each dimension with respect to point.
+    ///
+    /// Can be used to implement "zoom_to" functionality.
+    #[inline]
+    pub fn new_nonuniform_scaling_wrt_point(scaling: &Vector2<N>, pt: &Point2<N>) -> Self {
+        let _0 = N::zero();
+        let _1 = N::one();
+        Matrix3::new(
+            scaling.x,
+            _0,
+            pt.x - pt.x * scaling.x,
+            _0,
+            scaling.y,
+            pt.y - pt.y * scaling.y,
+            _0,
+            _0,
+            _1,
+        )
+    }
 }
 
 impl<N: RealField> Matrix4<N> {
@@ -90,6 +111,33 @@ impl<N: RealField> Matrix4<N> {
         Isometry::rotation_wrt_point(rot, pt).to_homogeneous()
     }
 
+    /// Creates a new homogeneous matrix that applies a scaling factor for each dimension with respect to point.
+    ///
+    /// Can be used to implement "zoom_to" functionality.
+    #[inline]
+    pub fn new_nonuniform_scaling_wrt_point(scaling: &Vector3<N>, pt: &Point3<N>) -> Self {
+        let _0 = N::zero();
+        let _1 = N::one();
+        Matrix4::new(
+            scaling.x,
+            _0,
+            _0,
+            pt.x - pt.x * scaling.x,
+            _0,
+            scaling.y,
+            _0,
+            pt.y - pt.y * scaling.y,
+            _0,
+            _0,
+            scaling.z,
+            pt.z - pt.z * scaling.z,
+            _0,
+            _0,
+            _0,
+            _1,
+        )
+    }
+
     /// Builds a 3D homogeneous rotation matrix from an axis and an angle (multiplied together).
     ///
     /// Returns the identity matrix if the given argument is zero.
diff --git a/tests/core/cg.rs b/tests/core/cg.rs
new file mode 100644
index 00000000..b061efbb
--- /dev/null
+++ b/tests/core/cg.rs
@@ -0,0 +1,59 @@
+use na::{Matrix3, Matrix4, Point2, Point3, Vector2, Vector3};
+
+/// See Example 3.4 of "Graphics and Visualization: Principles & Algorithms"
+/// by Theoharis, Papaioannou, Platis, Patrikalakis.
+#[test]
+fn test_scaling_wrt_point_1() {
+    let a = Point2::new(0.0, 0.0);
+    let b = Point2::new(1.0, 1.0);
+    let c = Point2::new(5.0, 2.0);
+
+    let scaling = Vector2::new(2.0, 2.0);
+    let scale_about = Matrix3::new_nonuniform_scaling_wrt_point(&scaling, &c);
+
+    let expected_a = Point2::new(-5.0, -2.0);
+    let expected_b = Point2::new(-3.0, 0.0);
+    let result_a = scale_about.transform_point(&a);
+    let result_b = scale_about.transform_point(&b);
+    let result_c = scale_about.transform_point(&c);
+
+    assert!(expected_a == result_a);
+    assert!(expected_b == result_b);
+    assert!(c == result_c);
+}
+
+/// Based on the same example as the test above.
+#[test]
+fn test_scaling_wrt_point_2() {
+    let a = Point3::new(0.0, 0.0, 1.0);
+    let b = Point3::new(1.0, 1.0, 1.0);
+    let c = Point3::new(5.0, 2.0, 1.0);
+
+    let scaling = Vector3::new(2.0, 2.0, 1.0);
+    let scale_about = Matrix4::new_nonuniform_scaling_wrt_point(&scaling, &c);
+
+    let expected_a = Point3::new(-5.0, -2.0, 1.0);
+    let expected_b = Point3::new(-3.0, 0.0, 1.0);
+
+    let result_a = scale_about.transform_point(&a);
+    let result_b = scale_about.transform_point(&b);
+    let result_c = scale_about.transform_point(&c);
+
+    assert!(expected_a == result_a);
+    assert!(expected_b == result_b);
+    assert!(c == result_c);
+}
+
+/// Based on https://github.com/emlowry/AiE/blob/50bae4068edb686cf8ffacdf6fab8e7cb22e7eb1/Year%201%20Classwork/MathTest/Matrix4x4TestGroup.cpp#L145
+#[test]
+fn test_scaling_wrt_point_3() {
+    let about = Point3::new(2.0, 1.0, -2.0);
+    let scale = Vector3::new(2.0, 0.5, -1.0);
+    let pt = Point3::new(1.0, 2.0, 3.0);
+    let scale_about = Matrix4::new_nonuniform_scaling_wrt_point(&scale, &about);
+
+    let expected = Point3::new(0.0, 1.5, -7.0);
+    let result = scale_about.transform_point(&pt);
+
+    assert!(result == expected);
+}
diff --git a/tests/core/mod.rs b/tests/core/mod.rs
index ec1c4e3e..7c738368 100644
--- a/tests/core/mod.rs
+++ b/tests/core/mod.rs
@@ -1,6 +1,7 @@
 #[cfg(feature = "abomonation-serialize")]
 mod abomonation;
 mod blas;
+mod cg;
 mod conversion;
 mod edition;
 mod empty;