Vanishing Points and Lines

A vanishing point is the image point toward which the images of parallel scene lines converge: the rails of a track, the edges of a corridor, the horizontal courses of a brick wall. A vanishing line, also called a horizon line when the plane is the ground, is the image line that contains the vanishing points of all directions parallel to one scene plane. Because they are images of points and lines at infinity, vanishing points and lines do not depend on where the camera is, only on how it is turned and on its internal parameters. That makes them a direct link between a single photograph and the camera’s orientation, its focal length and principal point, and the relative sizes of things in the scene. This article covers their geometry, how they are detected among cluttered line segments, calibration from orthogonal vanishing points, the Manhattan-world assumption, and single-view metrology.

Definition

The vanishing point of a direction $\mathbf{d}$ in space is the image of the point at infinity $(\mathbf{d}^\top, 0)^\top$ of that direction, that is, the image point where the images of all lines with direction $\mathbf{d}$ meet. The vanishing line of a plane is the image of the plane’s line at infinity; it contains the vanishing points of every direction parallel to the plane. The horizon of a photograph, in this geometric sense, is the vanishing line of the ground plane. The pinhole camera model shows how parallel lines converge, and the projective geometry article places vanishing points among the objects at infinity; the sections below build on both.

Intuition

Walk along a straight road and look ahead. The two edges of the road are parallel, yet in the picture they draw together, because equal widths farther away project smaller. Follow them indefinitely and their images approach a single point, which they never reach. Every other line parallel to the road, a kerb, a painted line, a row of lamp tops, heads to the same point, wherever it starts. A different set of parallel lines, such as the vertical edges of the buildings, heads to a different point, often far above or below the picture.

Turning the head moves the vanishing points across the image, while walking along the road leaves them where they are: they behave like stars, fixed by direction alone. Since the road surface, the pavements and every other horizontal plane contain many directions, their vanishing points line up along one line, the horizon, which passes through the point straight ahead at eye level. Anything in the scene that the horizon passes through is at the height of the eye.

Formal Definition

Vanishing points

Write the camera as $\mathbf{x} \sim \mathbf{K}\,[\mathbf{R} \mid \mathbf{t}]\,\mathbf{X}$, with intrinsic matrix $\mathbf{K}$ and world-to-camera rotation $\mathbf{R}$ (see camera parameters). The point at infinity of direction $\mathbf{d}$ has fourth coordinate zero, so the translation drops out:

$$\mathbf{v} \sim \mathbf{K}\,[\mathbf{R} \mid \mathbf{t}] \begin{pmatrix} \mathbf{d} \\ 0 \end{pmatrix} = \mathbf{K}\mathbf{R}\,\mathbf{d}.$$

The vanishing point depends on the direction, the rotation and the intrinsics, and on nothing else. Conversely, a vanishing point fixes a direction in camera coordinates, $\mathbf{R}\mathbf{d} \sim \mathbf{K}^{-1}\mathbf{v}$, up to sign: a vanishing point cannot tell a direction from its opposite. When $\mathbf{R}\mathbf{d}$ has no component along the optical axis, the third coordinate of $\mathbf{v}$ is zero and the vanishing point is itself at infinity, which is why lines parallel to the image plane stay parallel in the image. Working in homogeneous coordinates, with $\mathbf{v}$ scaled to unit length, handles both cases alike; Szeliski (2022, Ex. 7.14) suggests this representation on a sphere for vanishing points near infinity.

Under an ideal pinhole camera the images of lines with direction $\mathbf{d}$ pass exactly through $\mathbf{v}$, so in noise-free data any two of them determine it: $\mathbf{v} \sim \mathbf{l}_1 \times \mathbf{l}_2$ for homogeneous image lines $\mathbf{l}_1, \mathbf{l}_2$.

Vanishing lines

A direction $\mathbf{d}$ lies in a plane with unit normal $\mathbf{n}$ when $\mathbf{n}^\top\mathbf{d} = 0$. Substituting $\mathbf{d} \sim \mathbf{R}^\top\mathbf{K}^{-1}\mathbf{v}$ gives the condition $\mathbf{n}^\top\mathbf{R}^\top\mathbf{K}^{-1}\mathbf{v} = 0$, so, by our derivation, the vanishing points of the plane’s directions lie on the line

$$\mathbf{l} \sim \mathbf{K}^{-\top}\mathbf{R}\,\mathbf{n}.$$

All parallel planes share this line, since it depends only on their common normal. In practice the vanishing line of a plane is usually found as the line through two vanishing points of directions in it, $\mathbf{l} \sim \mathbf{v}_1 \times \mathbf{v}_2$. Criminisi, Reid and Zisserman (1999, §2) note the property that makes the horizon useful for measurement: a scene point whose image lies on the vanishing line of a reference plane is as far from that plane as the camera center, points imaged on one side of the line are farther from the plane, and points on the other side are closer.

Angles between directions

Two vanishing points determine the angle between their scene directions once $\mathbf{K}$ is known. Since $\mathbf{R}\mathbf{d}_i \sim \mathbf{K}^{-1}\mathbf{v}_i$ and rotations preserve angles, it follows that

$$\cos\theta = \frac{\mathbf{v}_1^\top \boldsymbol{\omega}\, \mathbf{v}_2}{\sqrt{\mathbf{v}_1^\top \boldsymbol{\omega}\, \mathbf{v}_1}\,\sqrt{\mathbf{v}_2^\top \boldsymbol{\omega}\, \mathbf{v}_2}}, \qquad \boldsymbol{\omega} = \mathbf{K}^{-\top}\mathbf{K}^{-1}.$$

The matrix $\boldsymbol{\omega}$ is the image of the absolute conic, which Szeliski (2022, §11.1) describes as commonly used for camera calibration, pointing to Hartley and Zisserman (2004, §8.5) for the theory. Orthogonal directions satisfy $\mathbf{v}_1^\top \boldsymbol{\omega}\, \mathbf{v}_2 = 0$, a linear equation in the entries of $\boldsymbol{\omega}$; this is the constraint behind every calibration method in the sections below.

The Manhattan world

Buildings, rooms and streets are dominated by three mutually orthogonal directions, two horizontal and one vertical. Szeliski (2022, §11.4.8) writes that finding an orthogonal set of vanishing points to fix a global orientation is often called invoking the Manhattan world assumption, after Coughlan and Yuille (1999), and that the generalization in which horizontal directions meet at other angles was named the Atlanta world by Schindler and Dellaert (2004). Under the Manhattan assumption an image has three dominant vanishing points whose back-projected directions $\mathbf{K}^{-1}\mathbf{v}_i$ are orthogonal; they are the columns of $\mathbf{K}\mathbf{R}$ up to scale, as the code below uses.

Properties

• Invariance to translation. Moving the camera without turning it leaves every vanishing point in place. A sequence of images from a translating camera therefore has fixed vanishing points, and changes in them measure rotation alone, which is the idea behind using vanishing points as a visual compass (the title of Košecká and Zhang’s 2002 paper, which could not be read for this article).
• Finite and infinite. A vanishing point is finite unless its direction is parallel to the image plane. The vertical vanishing point of a level camera is at infinity, so vertical edges stay parallel in the image; tilting the camera up or down brings it to a finite position above or below the image. Under affine cameras, all vanishing points are at infinity (see pinhole camera model).
• The principal point is the orthocenter. For a camera with square pixels and no skew, the principal point is the orthocenter of the triangle formed by three finite vanishing points of orthogonal directions; Szeliski (2022, §11.1.1) attributes the result to Caprile and Torre (1990) and to Hartley and Zisserman (2004, §8.6). It follows from the orthogonality constraint in the next section: subtracting the equations for two pairs that share a vertex shows that the principal point lies on that vertex’s altitude.
• Sensitivity grows with distance. A vanishing point far outside the image comes from image lines that are nearly parallel, so a small error in their angles moves it a long way. Szeliski (2022, §11.1.1) notes that the focal length estimated from two vanishing points becomes more accurate as they move closer to the image center, and recommends tilting a calibration pattern a good deal about the 45° axis for that reason; the example below shows the effect.
• Every line at camera height lies on the horizon. By the property of vanishing lines given above, any horizontal segment at the height of the camera’s center projects onto the horizon, whatever its direction. Simon, Fond and Berger (2018, §3.2) use this to detect the horizon directly as an alignment of segments. In our reading, the same fact makes such segments ambiguous for detectors: they are consistent with every horizontal vanishing point at once.
• Independence from the calibration of measurements. Ratios of heights above a plane, and ratios of lengths and areas on planes parallel to it, can be measured from the plane’s vanishing line and one vanishing point of another direction, without knowing $\mathbf{K}$ (Criminisi, Reid and Zisserman, 1999). This is the basis of single-view metrology.

Detecting Vanishing Points

Detection usually runs in three stages: find straight line segments, group them by the vanishing point they pass through, and refine each vanishing point from its group (Lezama et al., 2014, §1; Szeliski, 2022, §7.4). Lens distortion must be removed first, since it curves the straight edges the method depends on (see lens distortion).

Line segments

Segments come from an edge detection step followed by grouping of edge pixels into lines. Szeliski (2022, §7.4.2) describes the main options: the Hough transform, in which each edge pixel votes for the lines through it and peaks of the accumulator are lines; RANSAC, which proposes a line from two edge pixels and counts the others that fall on it; and bottom-up grouping of pixels with similar gradient orientation into line-support regions, the approach of the LSD detector of Grompone von Gioi, Jakubowicz et al. (2008), which Szeliski describes as widely used in practice for its speed, its rejection of texture and its open-source code. Every later stage inherits the quality of these segments: short or fragmented segments have poorly determined directions.

Grouping and estimation

• Gaussian-sphere accumulation. Barnard (1983), described by Lezama et al. (2014) as the seminal work on the problem, mapped image lines onto the Gaussian sphere of directions around a calibrated camera’s center and found vanishing points with a Hough transform on the sphere. In our reading, representing directions on the sphere treats finite and infinite vanishing points alike; Lezama et al. (2014, §1) note that it needs the calibration.
• Voting in the image. In the first edition of his book, Szeliski had pairs of lines vote for candidate vanishing point locations and refined the winners by robust least squares (Szeliski, 2022, §7.4.3).
• Hypothesize and verify. RANSAC samples two segments, takes their intersection as a candidate, and scores it by the segments consistent with it; the best candidate is refined and its segments removed before the next search. Tardif (2009) instead built a set of candidates from pairs of segments and computed consensus sets with the J-linkage algorithm (as Simon et al., 2018, §2 describe it); Zhou et al. (2019) call LSD followed by J-linkage probably the most widely used line-based method.
• Expectation–maximization. Košecká and Zhang (2002) estimated vanishing points and the probability that each segment belongs to each one by alternating the two, as Simon et al. (2018, §2) describe; Simon et al. point out that EM depends on its initialization. Szeliski (2022, §7.4.3) notes that alternating between assigning lines and refitting vanishing points may give better results than a strict two-stage pipeline.
• Manhattan-world inference. Coughlan and Yuille (2000) skipped segment grouping altogether. Their Bayesian model explains the gradient at each pixel as an edge of an $x$, $y$ or $z$ line, an outlier edge, or no edge, with priors of 0.02, 0.02, 0.02, 0.04 and 0.9, sums over these explanations rather than deciding them, and searches over a quantized set of camera orientations (three Euler angles) for the maximum of the posterior. The direction of each pixel’s gradient is modeled as within 4° (indoors) or 6° (outdoors) of the predicted line normal with probability 0.9. Their earlier ICCV paper had fixed the elevation and twist at zero, which they report was invalid for many of the later images. The same model flags edges that fit none of the three directions, which they used to pick out people against a building background.
• Orthogonality as a constraint. Searching for pairs or triplets of mutually orthogonal vanishing points, rather than for each separately, makes detection more robust (Szeliski, 2022, §7.4.3). Wildenauer and Hanbury (2012) solve for three orthogonal vanishing points and the focal length from a minimal sample of four lines inside RANSAC (as Simon et al., 2018, §2 describe it).
• A contrario detection. Lezama et al. (2014) map each segment to a point in the PClines dual space, where lines through a common point become collinear points, and detect those alignments with a point-alignment detector whose threshold bounds the expected number of false detections under a model of random points. The same detector first joins aligned segments into longer lines. They refine each candidate from the segments within 2° of it, using the angle between a segment and the line from its midpoint to the candidate as the consistency measure, then choose an orthogonal triplet (with a calibrated Manhattan scene) or the horizon (otherwise).
• Horizon first. Zhai, Workman and Jacobs (2016) reverse the usual order: a convolutional network estimates a distribution over horizon lines from the whole image, candidate horizons are sampled from it, and each is scored by the horizontal vanishing points found on it, without a Manhattan assumption. Simon et al. (2018) also detect the horizon first, as an a-contrario alignment of horizontal segments at the camera’s height, and then search for horizontal vanishing points only on it, which they report gives far fewer spurious vanishing points.
• Learned detectors. NeurVPS (Zhou et al., 2019) samples candidate vanishing points and classifies each with “conic convolutions”, convolutions applied in a coordinate frame aligned with the lines through the candidate, so that features are gathered along the directions that converge on it. It is trained on labeled images; its authors name the need for large amounts of training data as a drawback, and did not evaluate it on YUD, ECD or HLW (below) because, they write, no suitable public training set exists for them.

Refinement

Given the segments of one vanishing point, a common refinement minimizes the squared distances from the vanishing point to the lines through the segments, weighting each segment by its length (Lezama et al., 2014, §2.4). With the homogeneous lines $\mathbf{l}_i$ scaled so that their first two entries form a unit normal, this is the least squares problem

$$\hat{\mathbf{v}} = \arg\min_{\|\mathbf{v}\| = 1} \sum_i w_i\, (\mathbf{l}_i^\top \mathbf{v})^2,$$

whose solution is the right singular vector of the stacked, weighted lines with the smallest singular value (see singular value decomposition). Strictly, $\mathbf{l}_i^\top\mathbf{v}$ is a point-to-line distance only when $\mathbf{v}$ is normalized to third coordinate 1; constraining $\|\mathbf{v}\| = 1$ instead keeps the problem well posed when the vanishing point is at or near infinity, at the price of weighting near and far candidates differently. Criminisi, Reid and Zisserman (1999, §4) go further and treat the measured points as uncertain, propagating covariances to the final measurements.

Calibration from Orthogonal Vanishing Points

Assume square pixels, zero skew, focal length $f$ in pixels and principal point $\mathbf{c} = (c_x, c_y)$. For finite vanishing points $\mathbf{v}_i, \mathbf{v}_j$ in pixel coordinates, the orthogonality constraint $\mathbf{v}_i^\top \boldsymbol{\omega}\, \mathbf{v}_j = 0$ becomes (Szeliski, 2022, eq. 11.3)

$$(\mathbf{v}_i - \mathbf{c})^\top(\mathbf{v}_j - \mathbf{c}) + f^2 = 0, \qquad i \neq j.$$

• Focal length. With the principal point assumed at the image center, which Szeliski (2022, §11.1.1) calls often safe for rough modeling, each pair of finite orthogonal vanishing points gives $f^2$ directly. The two vanishing points must lie on opposite sides of the principal point, as seen from it, for $f^2$ to come out positive: a pair whose vectors from $\mathbf{c}$ make an acute angle cannot be orthogonal for any focal length.
• Principal point. With three finite vanishing points, the principal point is the orthocenter of their triangle (Caprile and Torre, 1990, as cited by Szeliski), after which any pair gives $f$. Szeliski adds that it is more accurate in practice to re-estimate the unknown intrinsics by nonlinear least squares.
• Rotation. Once $\mathbf{K}$ is known, the normalized vectors $\mathbf{K}^{-1}\mathbf{v}_i$ are the columns of $\mathbf{R}$, up to sign and to the assignment of world axes. Szeliski (2022, §11.1.1) suggests taking cross products of two of them for the third, or orthogonal Procrustes, the nearest rotation to the noisy columns, computed with a singular value decomposition. The signs must be fixed by convention (for example, “up” is the vertical direction pointing toward the top of the image), because a vanishing point does not distinguish $\mathbf{d}$ from $-\mathbf{d}$.
• Degenerate configurations. A vanishing point at infinity gives no equation of this form, since the constraint then involves only its direction. A level camera, with its vertical vanishing point at infinity, gives $f$ from the two horizontal vanishing points when the principal point is known, but no orthocenter. A view in which a facade is nearly parallel to the image plane pushes one vanishing point far away, where its position is poorly determined.

This is the single-image counterpart of target-based camera calibration: instead of a printed pattern of known geometry, the scene’s own right angles provide the constraints. Szeliski (2022, §11.1.1) lists Caprile and Torre (1990), Liebowitz and Zisserman (1998), Cipolla, Drummond and Robertson (1999), Criminisi, Reid and Zisserman (2000) and Hartley and Zisserman (2004) among the works that calibrate from vanishing points, and calls this rudimentary calibration (§11.6).

Single-View Metrology

Criminisi, Reid and Zisserman (1999) showed that a single perspective image, with no knowledge of $\mathbf{K}$ or the camera’s pose, supports three kinds of measurement once the vanishing line $\mathbf{l}$ of a reference plane (such as the ground) and the vanishing point $\mathbf{v}$ of a reference direction not parallel to it (such as the vertical) are known: distances between planes parallel to the reference plane, ratios of lengths and areas on those planes, and the camera’s position relative to the plane. All of them are affine measurements: distances between planes and the camera’s position are found up to one common scale factor, which a single known reference length fixes. Their ICCV paper frames the work as reversing the rules of perspective drawing set out by Alberti in 1435, and applies it to Piero della Francesca’s Flagellation, reconstructing the scene from the painting’s floor pattern and vertical lines.

The algebra is compact. Choose world axes with $X$ and $Y$ in the reference plane and $Z$ along the reference direction. The columns of the projection matrix are then the vanishing points of the three axes and the image of the origin, and Criminisi et al. write it as

$$\mathbf{P} = \begin{bmatrix} \mathbf{v}_X & \mathbf{v}_Y & \alpha\,\mathbf{v} & \hat{\mathbf{l}} \end{bmatrix},$$

with $\mathbf{v}_X$ and $\mathbf{v}_Y$ on the vanishing line, $\hat{\mathbf{l}} = \mathbf{l}/\|\mathbf{l}\|$ standing in for the origin, and $\alpha$ an unknown scale (their eq. 1). For a point $\mathbf{t}$ at height $Z$ above a base point $\mathbf{b}$ on the plane, imaged so that $\mathbf{b}$, $\mathbf{t}$ and $\mathbf{v}$ are collinear, the height satisfies (their eq. 4)

$$\alpha Z = -\frac{\|\mathbf{b} \times \mathbf{t}\|}{(\hat{\mathbf{l}}^\top\mathbf{b})\,\|\mathbf{v} \times \mathbf{t}\|}.$$

One reference object of known height gives $\alpha$; every other height then follows from its top and base points. Mapping between two parallel planes is a planar homology, $\tilde{\mathbf{H}} = \mathbf{I} + \alpha Z\, \mathbf{v}\hat{\mathbf{l}}^\top$ (their eq. 5), whose axis is the vanishing line and whose vertex is the vanishing point; it transfers points and areas from one plane to another for comparison. The camera center follows from $\mathbf{P}\mathbf{C} = \mathbf{0}$ by Cramer’s rule (their eq. 7).

Criminisi et al. also propagated the uncertainty of the vanishing line, the vanishing point and the marked points to each measurement (their §4), modeling the top and base points by covariance ellipses and fitting them, by maximum likelihood, to lie on a line through the vanishing point. In their examples, a person’s height was measured as 178.8 cm against a true 180 cm, with the person leaning on one foot (their Fig. 3), and from a security camera image corrected for radial distortion, as 190.6 cm against 190 cm, with a 3-sigma uncertainty of ±4.1 cm from one reference height that fell to ±2.9 cm with three (their Fig. 9). They demonstrate forensic measurement from surveillance images and 3D modeling from single photographs. The journal version (Criminisi, Reid and Zisserman, 2000) could not be read for this article; the details above are from the ICCV paper, as deposited in Oxford’s research archive. Szeliski (2022, §11.1.2) adds that if the focal length is first recovered from orthogonal vanishing points, a known distance on the ground fixes the scale and the reconstruction becomes metric.

Benchmarks and Reported Results

Ground-truth vanishing points are hard to define in real photographs, since which of many sets of parallel lines count is a matter of judgment. Evaluation on real images has therefore concentrated on the horizon (Simon et al., 2018, §2):

• Datasets. The York Urban Dataset (YUD) has 102 images, indoor and outdoor, satisfying the Manhattan assumption (“mostly”, in Simon et al., 2018, §5.2), with camera parameters and ground-truth vanishing point triplets; the Eurasian Cities Dataset (ECD) has 103 urban images in varied architectural styles, taken with different cameras, not all of which satisfy it (Lezama et al., 2014, §3; Zhai et al., 2016, §5.1). Simon et al. (2018) give ECD as 114 images at 1920 × 1080 in their text but 103 in the table of their Fig. 5. Horizon Lines in the Wild (HLW) is newer and harder, with a test set of about 2,000 images (Zhai et al., 2016, §5.1; Simon et al., 2018, count 2,018), many of them holiday-style photographs with wide ranges of camera roll and tilt.
• Protocol. Following the protocol of Barinova et al. (2010), as Lezama et al. describe it, the first 25 images of YUD and of ECD are used to set parameters and the rest to evaluate (Lezama et al., 2014, §3). Zhai et al. (2016, §5) instead use one set of parameters for all datasets.
• Metric. The horizon error is the largest distance, within the image, between the estimated and the true horizon, divided by the image height. Results are plotted as the cumulative distribution of this error over a dataset and summarized by the area under that curve for errors from 0 to 0.25 (Simon et al., 2018, §2; Zhai et al., 2016, §5).

Simon et al. (2018, Fig. 5) collect the area-under-curve values below. Such comparisons reuse published numbers rather than re-running every method; Lezama et al. (2014, §3.1), for example, say they used the values published by Xu et al. (2013).

Method	YUD AUC (%)	ECD AUC (%)	HLW AUC (%)
Košecká and Zhang (2002)	74.34	68.62
Tardif (2009)	84.62	81.33
Lezama et al. (2014)	94.51	89.20	52.59 (from Zhai et al., 2016)
Zhai et al. (2016)	94.78	90.80	58.24
Simon et al. (2018)	95.35	91.10	54.43

On HLW, where the camera roll and tilt vary widely and the horizon is often outside the image, the learned prior of Zhai et al. outperforms the purely geometric method of Simon et al.; Simon et al. trace the difference mainly to how widely the two methods sample candidate horizons. Lezama et al. (2014, §3.1) also report a maximum horizon error on YUD of 0.052, against 0.078 for Xu et al. (2013).

For learned detectors, Zhou et al. (2019) measure the angle between predicted and true vanishing directions and report “angle accuracy”, the area under the curve of the fraction of predictions within each angle, up to a threshold. On the synthetic SU3 wireframe dataset, where vanishing points come from the CAD models, NeurVPS with four conic convolution layers reached a median angular error of 0.09°, against 0.21° for LSD with J-linkage (their Table 1). On ScanNet, with three vanishing points per image derived from reconstructed meshes, NeurVPS’s mean and median errors were 4.5° and 3.0° against 12.6° and 11.8° for LSD with J-linkage (their Table 3); their text describes this as a six-fold and four-fold reduction, while the table gives about 2.8-fold for the mean. They note that ScanNet’s ground truth is itself inaccurate.

Examples

The code renders a 640 × 480 image of a building corner, two facades with window grids standing on paved ground, under a camera with $f = 800$ pixels and the principal point at the image center. Randomly oriented triangles add edges that belong to no vanishing point, and Gaussian noise is added. Line segments are found with scikit-image’s Canny detector and probabilistic Hough transform. Three vanishing points are found by sequential RANSAC, scored by the length of consistent segments, with Lezama et al.’s angular consistency measure and a 2° threshold, and refined by weighted least squares. From them the code computes the focal length with the principal point assumed at the image center and, separately, with the principal point estimated as the orthocenter, since the choice is open and changes the result; then the rotation, by orthogonal Procrustes; and finally the heights of a 1.8 m pole and of the 9 m building from a 3 m reference pole by Criminisi et al.’s eq. 4, with top and base points “clicked” with 1-pixel noise. Three camera poses are tested: an oblique view from above, a view in which one facade is nearly parallel to the image plane, and a level camera, whose vertical vanishing point is at infinity. Each pose is repeated with seeds 0 to 9, which change the clutter, the image noise, the Hough sampling, the RANSAC samples and the clicks, and the median and maximum over the ten runs (over the runs in which a quantity is defined) are printed. The world axes are named by matching the detected vanishing points to the true ones, which only labels them; the errors are angles between back-projected directions computed with the true $\mathbf{K}$. It uses NumPy 2.3.5 and scikit-image 0.25.2 and runs in a few seconds.

import itertools
import numpy as np
from skimage.draw import polygon
from skimage.feature import canny
from skimage.transform import probabilistic_hough_line

W, H, F = 640, 480, 800.0
K = np.array([[F, 0, W / 2], [0, F, H / 2], [0, 0, 1.0]])       # true intrinsics


def camera(C, heading, pitch, roll):
    """World-to-camera rotation (x right, y down, z forward) and centre; world Z is up."""
    h, p, r = np.radians([heading, pitch, roll])
    z = np.array([np.cos(p) * np.cos(h), np.cos(p) * np.sin(h), np.sin(p)])
    x = np.cross(z, [0, 0, 1.0]); x /= np.linalg.norm(x)
    Rr = np.array([[np.cos(r), -np.sin(r), 0], [np.sin(r), np.cos(r), 0], [0, 0, 1]])
    return Rr @ np.array([x, np.cross(z, x), z]), np.array(C, float)


POSES = {"oblique":      dict(C=(-14, -18, 14), heading=50, pitch=-25, roll=3),
         "near-frontal": dict(C=(-3, -24, 14), heading=80, pitch=-25, roll=3),
         "level":        dict(C=(-20, -26, 4.5), heading=50, pitch=0, roll=3)}
POLES = {"reference": ((-3.0, -3.0), 3.0), "person": ((-1.0, -5.0), 1.8)}


def project(R, C, X):
    x = (K @ (R @ (np.atleast_2d(X) - C).T)).T
    return x[:, :2] / x[:, 2:]


def fill(img, R, C, poly3d, value, near=0.5):
    """Clip a planar 3D polygon to depth > near, project it and paint it."""
    P, out = (R @ (np.asarray(poly3d, float) - C).T).T, []
    for a, b in zip(P, np.roll(P, -1, axis=0)):
        if a[2] >= near: out.append(a)
        if (a[2] >= near) != (b[2] >= near):
            out.append(a + (near - a[2]) / (b[2] - a[2]) * (b - a))
    if len(out) >= 3:
        x = (K @ np.array(out).T).T
        img[polygon(x[:, 1] / x[:, 2], x[:, 0] / x[:, 2], img.shape)] = value


def facade(axis, a0, a1, z0, z1):   # rectangle in the plane Y = 0 ("A") or X = 0 ("B")
    q = [(a0, z0), (a1, z0), (a1, z1), (a0, z1)]
    return [(a, 0, z) for a, z in q] if axis == "A" else [(0, a, z) for a, z in q]


def render(R, C, rng):
    """A building corner (facade A along X, B along Y, 9 m tall) on paved ground."""
    img = np.zeros((H, W))
    fill(img, R, C, [(-200, -200, 0), (200, -200, 0), (200, 200, 0), (-200, 200, 0)], 0.55)
    for k in np.arange(-24, 0, 2.0):                               # paving joints
        fill(img, R, C, [(k, -24, 0), (k + .1, -24, 0), (k + .1, -.5, 0), (k, -.5, 0)], 0.42)
        fill(img, R, C, [(-24, k, 0), (-.5, k, 0), (-.5, k + .1, 0), (-24, k + .1, 0)], 0.42)
    if C[2] > 9:
        fill(img, R, C, [(0, 0, 9), (14, 0, 9), (14, 12, 9), (0, 12, 9)], 0.62)   # roof
    fill(img, R, C, facade("A", 0, 14, 0, 9), 0.72)
    fill(img, R, C, facade("B", 0, 12, 0, 9), 0.88)
    for z0 in (1.0, 4.0, 7.0):                                     # windows
        for k in range(5): fill(img, R, C, facade("A", 1 + 2.6 * k, 2.2 + 2.6 * k, z0, z0 + 1.6), 0.25)
        for k in range(4): fill(img, R, C, facade("B", 1 + 2.8 * k, 2.2 + 2.8 * k, z0, z0 + 1.6), 0.30)
    for _ in range(10):                                            # clutter: random triangles
        c, ang, rad = rng.uniform([0, H * .3], [W, H]), rng.uniform(0, 2 * np.pi, 3), rng.uniform(10, 35, 3)
        img[polygon(c[1] + rad * np.sin(ang), c[0] + rad * np.cos(ang), img.shape)] = rng.uniform(.15, .95)
    e = R[0].copy(); e[2] = 0; e /= np.linalg.norm(e)              # poles 0.2 m wide
    for (x, y), h in POLES.values():
        p = np.array([x, y, 0.0])
        fill(img, R, C, [p - .1 * e, p + .1 * e, p + .1 * e + [0, 0, h], p - .1 * e + [0, 0, h]], 0.08)
    return np.clip(img + rng.normal(0, 0.02, img.shape), 0, 1)


def unit_lines(seg):
    L = np.cross(np.c_[seg[:, 0], np.ones(len(seg))], np.c_[seg[:, 1], np.ones(len(seg))])
    return L / np.linalg.norm(L[:, :2], axis=1, keepdims=True)


def consistency(seg, v):
    """Angle (deg) between each segment and the line joining its midpoint to v."""
    d, to_v = seg[:, 1] - seg[:, 0], v[:2] - seg.mean(axis=1) * v[2]
    c = np.abs(np.sum(d * to_v, 1)) / (np.linalg.norm(d, axis=1) * np.linalg.norm(to_v, axis=1) + 1e-300)
    return np.degrees(np.arccos(np.clip(c, 0, 1)))


def refine(seg, w):
    """Homogeneous point closest to all lines (weighted least squares on the unit sphere)."""
    return np.linalg.svd(unit_lines(seg) * np.sqrt(w)[:, None])[2][-1]


def ransac_vps(seg_px, rng, n_vp=3, iters=500, thr=2.0):
    """Sequential RANSAC: two segments propose a VP; the support is the inlier length."""
    T = np.array([[W, 0, W / 2], [0, W, H / 2], [0, 0, 1.0]])     # conditioning
    seg = (seg_px - [W / 2, H / 2]) / W
    length, L = np.linalg.norm(seg[:, 1] - seg[:, 0], axis=1), unit_lines(seg)
    alive, vps = np.ones(len(seg), bool), []
    for _ in range(n_vp):
        idx = np.flatnonzero(alive)
        hyp = [np.cross(*L[rng.choice(idx, 2, replace=False)]) for _ in range(iters)]
        v = max(hyp, key=lambda h: length[idx][consistency(seg[idx], h) < thr].sum())
        for _ in range(2):                                         # refit on the inliers
            inl = idx[consistency(seg[idx], v) < thr]
            v = refine(seg[inl], length[inl])
        alive[idx[consistency(seg[idx], v) < thr]] = False
        vps.append(T @ v / np.linalg.norm(T @ v))
    return np.array(vps)


def angle_err(a, b):          # angle between the 3D directions the VPs imply (true K)
    a, b = np.linalg.solve(K, a), np.linalg.solve(K, b)
    return np.degrees(np.arccos(min(1.0, abs(a @ b) / np.linalg.norm(a) / np.linalg.norm(b))))


def focal(vps, c):
    """f^2 = -(v_i - c).(v_j - c) for each pair of finite orthogonal VPs, averaged."""
    f2 = [-(a[:2] - c * a[2]) @ (b[:2] - c * b[2]) / (a[2] * b[2])
          for a, b in itertools.combinations(vps, 2) if abs(a[2] * b[2]) > 1e-9]
    f2 = [x for x in f2 if x > 0]
    return np.sqrt(np.mean(f2)) if f2 else np.nan


def orthocenter(vps):
    if np.any(np.abs(vps[:, 2]) < 1e-9):
        return np.full(2, np.nan)
    p = vps[:, :2] / vps[:, 2:]
    A = np.array([p[2] - p[1], p[2] - p[0]])
    return np.linalg.solve(A, [(p[2] - p[1]) @ p[0], (p[2] - p[0]) @ p[1]])


def rotation(vps, f, c):
    """Columns K^-1 v (signs: horizontals forward, vertical up), then the nearest rotation."""
    Kh = np.array([[f, 0, c[0]], [0, f, c[1]], [0, 0, 1.0]])
    M = np.array([np.linalg.solve(Kh, v) / np.linalg.norm(np.linalg.solve(Kh, v)) for v in vps]).T
    M[:, :2] *= np.sign(M[2, :2]); M[:, 2] *= -np.sign(M[1, 2])
    U, _, Vt = np.linalg.svd(M)
    return U @ np.diag([1, 1, np.linalg.det(U @ Vt)]) @ Vt


def rot_err(Ra, Rb):
    return np.degrees(np.arccos(np.clip((np.trace(Ra @ Rb.T) - 1) / 2, -1, 1)))


def height(vx, vy, vz, ref, z_ref, target):
    """Single-view metrology (Criminisi et al.): alpha Z = -|b x t| / ((l.b) |v x t|)."""
    l = np.cross(vx, vy); l /= np.linalg.norm(l)                   # vanishing line of the ground
    def aZ(b, t):
        b, t = np.r_[b, 1.0], np.r_[t, 1.0]
        return -np.linalg.norm(np.cross(b, t)) / ((l @ b) * np.linalg.norm(np.cross(vz, t)))
    return aZ(*target) / (aZ(*ref) / z_ref)


for name, pose in POSES.items():
    R, C = camera(**pose)
    true_vps = (K @ R).T                                          # K R e_X, K R e_Y, K R e_Z
    rows = []
    for seed in range(10):
        rng = np.random.default_rng(seed)
        img = render(R, C, rng)
        segs = np.array(probabilistic_hough_line(canny(img, sigma=2.0), threshold=10,
                                                 line_length=25, line_gap=3, rng=rng), float)
        est = ransac_vps(segs, rng)
        perm = min(itertools.permutations(range(3)),               # name the axes (X, Y, Z)
                   key=lambda p: sum(angle_err(est[p[k]], true_vps[k]) for k in range(3)))
        vps = est[list(perm)]
        vp_err = [angle_err(vps[k], true_vps[k]) for k in range(3)]
        c0, cb = np.array([W / 2, H / 2]), orthocenter(vps)
        fa, fb = focal(vps, c0), focal(vps, cb)
        ra = rot_err(rotation(vps, fa, c0), R) if np.isfinite(fa) else np.nan
        rb = rot_err(rotation(vps, fb, cb), R) if np.isfinite(fb) else np.nan
        noisy = lambda X: project(R, C, X) + rng.normal(0, 1.0, (2, 2))   # clicked points
        (xr, yr), hr = POLES["reference"]; (xp, yp), hp = POLES["person"]
        ref = noisy([[xr, yr, 0], [xr, yr, hr]])
        h_person = height(*vps, ref, hr, noisy([[xp, yp, 0], [xp, yp, hp]]))
        h_bldg = height(*vps, ref, hr, noisy([[0, 0, 0], [0, 0, 9.0]]))
        rows.append([len(segs), *vp_err, 100 * abs(fa / F - 1), 100 * abs(fb / F - 1),
                     np.linalg.norm(cb - c0), ra, rb,
                     100 * abs(h_person / hp - 1), 100 * abs(h_bldg / 9 - 1)])
    rows = np.array(rows)
    where = [str(np.round(v[:2] / v[2])) if abs(v[2]) > 1e-9 else "at infinity" for v in true_vps]
    print(f"{name}: true VPs X {where[0]}, Y {where[1]}, Z {where[2]}")
    labels = ["segments", "VP X error (deg)", "VP Y error (deg)", "VP Z error (deg)", "f error, pp at centre (%)",
              "f error, pp at orthocentre (%)", "orthocentre error (px)",
              "rotation error, pp at centre (deg)", "rotation error, orthocentre (deg)",
              "person height error (%)", "building height error (%)"]
    for lab, col in zip(labels, rows.T):
        ok = np.isfinite(col)                                  # nan: quantity undefined
        stats = f"median {np.median(col[ok]):8.2f}   max {np.max(col[ok]):8.2f}" if ok.any() else "undefined in every run"
        print(f"  {lab:36s} {stats}" + (f"   (undefined in {np.sum(~ok)} of 10)" if 0 < np.sum(~ok) else ""))

Output:

oblique: true VPs X [1390.  -77.], Y [-400. -171.], Z [ 230. 1953.]
  segments                             median    86.50   max   106.00
  VP X error (deg)                     median     0.30   max     0.60
  VP Y error (deg)                     median     0.31   max     0.60
  VP Z error (deg)                     median     0.47   max     0.83
  f error, pp at centre (%)            median     0.80   max     1.85
  f error, pp at orthocentre (%)       median     0.78   max     1.17
  orthocentre error (px)               median     7.99   max    14.87
  rotation error, pp at centre (deg)   median     0.36   max     0.49
  rotation error, orthocentre (deg)    median     0.56   max     0.94
  person height error (%)              median     1.92   max     7.17
  building height error (%)            median     1.46   max     2.64
near-frontal: true VPs X [5339.  129.], Y [ 184. -141.], Z [ 230. 1953.]
  segments                             median    91.50   max   100.00
  VP X error (deg)                     median     1.30   max    83.90
  VP Y error (deg)                     median     6.06   max    35.05
  VP Z error (deg)                     median     0.97   max     1.43
  f error, pp at centre (%)            median    15.86   max    44.22   (undefined in 1 of 10)
  f error, pp at orthocentre (%)       median     3.30   max     9.83   (undefined in 5 of 10)
  orthocentre error (px)               median    95.01   max 29741.05
  rotation error, pp at centre (deg)   median     2.08   max    45.55   (undefined in 1 of 10)
  rotation error, orthocentre (deg)    median     2.97   max     5.91   (undefined in 5 of 10)
  person height error (%)              median     9.55   max   418.93
  building height error (%)            median    14.47   max 303771.17
level: true VPs X [1272.  290.], Y [-350.  205.], Z at infinity
  segments                             median    80.00   max    88.00
  VP X error (deg)                     median     1.61   max     3.44
  VP Y error (deg)                     median     0.50   max    39.59
  VP Z error (deg)                     median     0.62   max     1.41
  f error, pp at centre (%)            median     9.65   max    71.57
  f error, pp at orthocentre (%)       median    23.64   max    56.97   (undefined in 5 of 10)
  orthocentre error (px)               median  1120.31   max 11072.88
  rotation error, pp at centre (deg)   median     1.21   max    15.68
  rotation error, orthocentre (deg)    median    25.20   max    37.95   (undefined in 5 of 10)
  person height error (%)              median     2.47   max    56.66
  building height error (%)            median     1.49   max   238.26

In the oblique view all three vanishing points are recovered to within 0.83° in every run, the focal length to within 1.9%, the rotation to within 1°, and the two heights to within about 7%. Choosing the principal point at the center or at the orthocenter makes little difference to $f$ here, but the orthocenter itself is off by 8 pixels in the median run and the rotation computed with it is slightly worse, which is consistent with Szeliski’s advice to assume the principal point or re-estimate it by nonlinear least squares rather than rely on the orthocenter alone.

The near-frontal view shows why the sensitivity property matters. The $X$ vanishing point lies more than 5,000 pixels from the image center, so the horizontal edges of facade A are nearly parallel in the image, and the $Y$ facade is seen at a grazing angle. The median focal length error rises to about 16% with the principal point at the center, and in one run no pair of vanishing points was on opposite sides of the center, so no positive $f^2$ existed. Errors of a few degrees in the horizontal vanishing points also move the vanishing line, and the median height errors rise to 10–14%. In our reading, the extreme maxima come from runs in which a vanishing point was badly misplaced; the largest $X$ error printed is 84°, which means a cluster of clutter or of another family’s segments won the search.

The level camera is a different degeneracy. Its vertical vanishing point is at infinity, so the orthocenter is meaningless: the estimated vertical vanishing point is finite only because of noise, very far away, and the orthocenter computed from it lies a median of about 1,100 pixels from the true principal point. Only the focal length with the principal point at the center is usable, with a median error near 10%; since the code averages $f^2$ over every pair that gives a positive value, pairs with the noisy vertical vanishing point can enter that average along with the horizontal pair. In our reading, part of the trouble is the property noted above: with the camera at 4.5 m, the window edges near that height lie close to the horizon and are consistent with both horizontal vanishing points, so they can contaminate the clusters. Metrology, which needs no calibration, stays within 2.5% in the median run in the oblique and level poses, though not in the near-frontal one, where it depends on the poorly placed horizontal vanishing points; in our reading, its large maxima come from runs in which a horizontal vanishing point was badly placed.

The script was run twice before the version shown. The first run stopped with an error in the near-frontal pose, because in one run no focal length existed and the rotation step was called with an undefined value; the script was changed to mark such quantities as undefined and to report how many runs that happened in. The second run printed only the largest of the three vanishing point errors in each run; the printout was then changed to give each vanishing point’s error separately. The computation is otherwise unchanged from the first run.

Common Misconceptions

• “A vanishing point is where parallel lines meet in the scene.” Parallel lines never meet in the scene. The vanishing point is the image of their common point at infinity, and it is the same for all of them because it depends only on their direction.
• “Vanishing points are always inside the image.” They are often far outside it, and at infinity when the direction is parallel to the image plane. Representing them as unit homogeneous vectors handles all cases.
• “The horizon is where the sky meets the ground.” The geometric horizon is the vanishing line of the ground plane, the set of points at the camera’s height on a level ground. Hills, buildings or a tilted ground plane can put the visible skyline well away from it.
• “Two vanishing points always give the focal length.” Only for finite, orthogonal directions, with the principal point known, and only when the two vanishing points lie on opposite sides of it. Accuracy falls quickly as either moves away from the image center, as the near-frontal example shows.
• “Single-view measurement needs a calibrated camera.” Criminisi, Reid and Zisserman (1999) measure heights and area ratios from a vanishing line, a vanishing point and one reference length, without $\mathbf{K}$. Calibration is needed only to make the whole reconstruction metric.
• “Any set of concurrent image lines marks a vanishing point.” Lines through one image point are not necessarily parallel in space, and clutter produces chance intersections. Detectors guard against this with consistency thresholds, statistical significance tests (Lezama et al., 2014) or priors such as orthogonality and the horizon.

Where It Is Used

• Camera calibration and orientation. Vanishing points give the focal length, principal point and rotation from one image of a structured scene, as above, complementing pattern-based camera calibration, and they give the camera pose’s rotation relative to a building or room. In image sequences they serve as a rotation reference, as in Košecká and Zhang’s “video compass”.
• Single-view reconstruction and measurement. Heights of people in surveillance images, architectural models from single photographs, and reconstructions of paintings (Criminisi, Reid and Zisserman, 1999); Szeliski (2022, §11.1.2) also points to interactive piecewise-planar modeling systems that use vanishing points to fix line directions and plane normals. Once orthogonal vanishing points are found, rectangles on facades and in rooms can be searched for, mainly in architectural scenes (Szeliski, 2022, §7.4.3).
• Rectification of planes. It follows from the projective hierarchy that a homography mapping a plane’s vanishing line back to the line at infinity removes the plane’s projective distortion and leaves an affine one, since fixing the line at infinity recovers affine geometry (see projective geometry); this is a first step in rectifying a facade or a floor for image warping.
• Driving and sports scenes. Lanes converge toward a vanishing point, which lane detection can use to restrict line searches; inverse perspective mapping in BEV perception undoes the same convergence; and sports field registration can estimate vanishing points from field markings (see sports analytics).
• Other uses named in the detection literature. Zhou et al. (2019) list camera calibration, 3D reconstruction, photo forensics, object detection, wireframe parsing and autonomous driving; Simon et al. (2018) add robot navigation and augmented reality.