diff --git a/pyproject.toml b/pyproject.toml
index da2537b2137..54a4bcaec0c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -146,7 +146,6 @@ exclude = [
   '^src/PIL/DdsImagePlugin.py$',
   '^src/PIL/FpxImagePlugin.py$',
   '^src/PIL/Image.py$',
-  '^src/PIL/ImageMath.py$',
   '^src/PIL/ImageMorph.py$',
   '^src/PIL/ImageQt.py$',
   '^src/PIL/ImageShow.py$',
diff --git a/src/PIL/Image.py b/src/PIL/Image.py
index c56da545803..ac13c6c0c25 100644
--- a/src/PIL/Image.py
+++ b/src/PIL/Image.py
@@ -873,7 +873,7 @@ def verify(self):
 
     def convert(
         self, mode=None, matrix=None, dither=None, palette=Palette.WEB, colors=256
-    ):
+    ) -> Image:
         """
         Returns a converted copy of this image. For the "P" mode, this
         method translates pixels through the palette.  If mode is
@@ -1305,7 +1305,7 @@ def getbands(self):
         """
         return ImageMode.getmode(self.mode).bands
 
-    def getbbox(self, *, alpha_only=True):
+    def getbbox(self, *, alpha_only: bool = True) -> tuple[int, int, int, int]:
         """
         Calculates the bounding box of the non-zero regions in the
         image.
diff --git a/src/PIL/ImageMath.py b/src/PIL/ImageMath.py
index b77f4bce567..bc3318c04f4 100644
--- a/src/PIL/ImageMath.py
+++ b/src/PIL/ImageMath.py
@@ -17,6 +17,8 @@
 from __future__ import annotations
 
 import builtins
+from types import CodeType
+from typing import Any
 
 from . import Image, _imagingmath
 
@@ -24,10 +26,10 @@
 class _Operand:
     """Wraps an image operand, providing standard operators"""
 
-    def __init__(self, im):
+    def __init__(self, im: Image.Image):
         self.im = im
 
-    def __fixup(self, im1):
+    def __fixup(self, im1: _Operand | float) -> Image.Image:
         # convert image to suitable mode
         if isinstance(im1, _Operand):
             # argument was an image.
@@ -45,122 +47,131 @@ def __fixup(self, im1):
             else:
                 return Image.new("F", self.im.size, im1)
 
-    def apply(self, op, im1, im2=None, mode=None):
-        im1 = self.__fixup(im1)
+    def apply(
+        self,
+        op: str,
+        im1: _Operand | float,
+        im2: _Operand | float | None = None,
+        mode: str | None = None,
+    ) -> _Operand:
+        im_1 = self.__fixup(im1)
         if im2 is None:
             # unary operation
-            out = Image.new(mode or im1.mode, im1.size, None)
-            im1.load()
+            out = Image.new(mode or im_1.mode, im_1.size, None)
+            im_1.load()
             try:
-                op = getattr(_imagingmath, op + "_" + im1.mode)
+                op = getattr(_imagingmath, op + "_" + im_1.mode)
             except AttributeError as e:
                 msg = f"bad operand type for '{op}'"
                 raise TypeError(msg) from e
-            _imagingmath.unop(op, out.im.id, im1.im.id)
+            _imagingmath.unop(op, out.im.id, im_1.im.id)
         else:
             # binary operation
-            im2 = self.__fixup(im2)
-            if im1.mode != im2.mode:
+            im_2 = self.__fixup(im2)
+            if im_1.mode != im_2.mode:
                 # convert both arguments to floating point
-                if im1.mode != "F":
-                    im1 = im1.convert("F")
-                if im2.mode != "F":
-                    im2 = im2.convert("F")
-            if im1.size != im2.size:
+                if im_1.mode != "F":
+                    im_1 = im_1.convert("F")
+                if im_2.mode != "F":
+                    im_2 = im_2.convert("F")
+            if im_1.size != im_2.size:
                 # crop both arguments to a common size
-                size = (min(im1.size[0], im2.size[0]), min(im1.size[1], im2.size[1]))
-                if im1.size != size:
-                    im1 = im1.crop((0, 0) + size)
-                if im2.size != size:
-                    im2 = im2.crop((0, 0) + size)
-            out = Image.new(mode or im1.mode, im1.size, None)
-            im1.load()
-            im2.load()
+                size = (
+                    min(im_1.size[0], im_2.size[0]),
+                    min(im_1.size[1], im_2.size[1]),
+                )
+                if im_1.size != size:
+                    im_1 = im_1.crop((0, 0) + size)
+                if im_2.size != size:
+                    im_2 = im_2.crop((0, 0) + size)
+            out = Image.new(mode or im_1.mode, im_1.size, None)
+            im_1.load()
+            im_2.load()
             try:
-                op = getattr(_imagingmath, op + "_" + im1.mode)
+                op = getattr(_imagingmath, op + "_" + im_1.mode)
             except AttributeError as e:
                 msg = f"bad operand type for '{op}'"
                 raise TypeError(msg) from e
-            _imagingmath.binop(op, out.im.id, im1.im.id, im2.im.id)
+            _imagingmath.binop(op, out.im.id, im_1.im.id, im_2.im.id)
         return _Operand(out)
 
     # unary operators
-    def __bool__(self):
+    def __bool__(self) -> bool:
         # an image is "true" if it contains at least one non-zero pixel
         return self.im.getbbox() is not None
 
-    def __abs__(self):
+    def __abs__(self) -> _Operand:
         return self.apply("abs", self)
 
-    def __pos__(self):
+    def __pos__(self) -> _Operand:
         return self
 
-    def __neg__(self):
+    def __neg__(self) -> _Operand:
         return self.apply("neg", self)
 
     # binary operators
-    def __add__(self, other):
+    def __add__(self, other: _Operand | float) -> _Operand:
         return self.apply("add", self, other)
 
-    def __radd__(self, other):
+    def __radd__(self, other: _Operand | float) -> _Operand:
         return self.apply("add", other, self)
 
-    def __sub__(self, other):
+    def __sub__(self, other: _Operand | float) -> _Operand:
         return self.apply("sub", self, other)
 
-    def __rsub__(self, other):
+    def __rsub__(self, other: _Operand | float) -> _Operand:
         return self.apply("sub", other, self)
 
-    def __mul__(self, other):
+    def __mul__(self, other: _Operand | float) -> _Operand:
         return self.apply("mul", self, other)
 
-    def __rmul__(self, other):
+    def __rmul__(self, other: _Operand | float) -> _Operand:
         return self.apply("mul", other, self)
 
-    def __truediv__(self, other):
+    def __truediv__(self, other: _Operand | float) -> _Operand:
         return self.apply("div", self, other)
 
-    def __rtruediv__(self, other):
+    def __rtruediv__(self, other: _Operand | float) -> _Operand:
         return self.apply("div", other, self)
 
-    def __mod__(self, other):
+    def __mod__(self, other: _Operand | float) -> _Operand:
         return self.apply("mod", self, other)
 
-    def __rmod__(self, other):
+    def __rmod__(self, other: _Operand | float) -> _Operand:
         return self.apply("mod", other, self)
 
-    def __pow__(self, other):
+    def __pow__(self, other: _Operand | float) -> _Operand:
         return self.apply("pow", self, other)
 
-    def __rpow__(self, other):
+    def __rpow__(self, other: _Operand | float) -> _Operand:
         return self.apply("pow", other, self)
 
     # bitwise
-    def __invert__(self):
+    def __invert__(self) -> _Operand:
         return self.apply("invert", self)
 
-    def __and__(self, other):
+    def __and__(self, other: _Operand | float) -> _Operand:
         return self.apply("and", self, other)
 
-    def __rand__(self, other):
+    def __rand__(self, other: _Operand | float) -> _Operand:
         return self.apply("and", other, self)
 
-    def __or__(self, other):
+    def __or__(self, other: _Operand | float) -> _Operand:
         return self.apply("or", self, other)
 
-    def __ror__(self, other):
+    def __ror__(self, other: _Operand | float) -> _Operand:
         return self.apply("or", other, self)
 
-    def __xor__(self, other):
+    def __xor__(self, other: _Operand | float) -> _Operand:
         return self.apply("xor", self, other)
 
-    def __rxor__(self, other):
+    def __rxor__(self, other: _Operand | float) -> _Operand:
         return self.apply("xor", other, self)
 
-    def __lshift__(self, other):
+    def __lshift__(self, other: _Operand | float) -> _Operand:
         return self.apply("lshift", self, other)
 
-    def __rshift__(self, other):
+    def __rshift__(self, other: _Operand | float) -> _Operand:
         return self.apply("rshift", self, other)
 
     # logical
@@ -170,46 +181,46 @@ def __eq__(self, other):
     def __ne__(self, other):
         return self.apply("ne", self, other)
 
-    def __lt__(self, other):
+    def __lt__(self, other: _Operand | float) -> _Operand:
         return self.apply("lt", self, other)
 
-    def __le__(self, other):
+    def __le__(self, other: _Operand | float) -> _Operand:
         return self.apply("le", self, other)
 
-    def __gt__(self, other):
+    def __gt__(self, other: _Operand | float) -> _Operand:
         return self.apply("gt", self, other)
 
-    def __ge__(self, other):
+    def __ge__(self, other: _Operand | float) -> _Operand:
         return self.apply("ge", self, other)
 
 
 # conversions
-def imagemath_int(self):
+def imagemath_int(self: _Operand) -> _Operand:
     return _Operand(self.im.convert("I"))
 
 
-def imagemath_float(self):
+def imagemath_float(self: _Operand) -> _Operand:
     return _Operand(self.im.convert("F"))
 
 
 # logical
-def imagemath_equal(self, other):
+def imagemath_equal(self: _Operand, other: _Operand | float | None) -> _Operand:
     return self.apply("eq", self, other, mode="I")
 
 
-def imagemath_notequal(self, other):
+def imagemath_notequal(self: _Operand, other: _Operand | float | None) -> _Operand:
     return self.apply("ne", self, other, mode="I")
 
 
-def imagemath_min(self, other):
+def imagemath_min(self: _Operand, other: _Operand | float | None) -> _Operand:
     return self.apply("min", self, other)
 
 
-def imagemath_max(self, other):
+def imagemath_max(self: _Operand, other: _Operand | float | None) -> _Operand:
     return self.apply("max", self, other)
 
 
-def imagemath_convert(self, mode):
+def imagemath_convert(self: _Operand, mode: str) -> _Operand:
     return _Operand(self.im.convert(mode))
 
 
@@ -219,7 +230,7 @@ def imagemath_convert(self, mode):
         ops[k[10:]] = v
 
 
-def eval(expression, _dict={}, **kw):
+def eval(expression: str, _dict: dict[str, Any] = {}, **kw: Any) -> Any:
     """
     Evaluates an image expression.
 
@@ -247,7 +258,7 @@ def eval(expression, _dict={}, **kw):
 
     compiled_code = compile(expression, "<string>", "eval")
 
-    def scan(code):
+    def scan(code: CodeType) -> None:
         for const in code.co_consts:
             if type(const) is type(compiled_code):
                 scan(const)
diff --git a/src/PIL/_imagingmath.pyi b/src/PIL/_imagingmath.pyi
new file mode 100644
index 00000000000..b0235555dc5
--- /dev/null
+++ b/src/PIL/_imagingmath.pyi
@@ -0,0 +1,5 @@
+from __future__ import annotations
+
+from typing import Any
+
+def __getattr__(name: str) -> Any: ...