1 Introduction

In the paper, we give a method for the reconstruction of the center of a projection. The process is based on distances and incidence relations in the following sense. The distance part means that we know the distance between the point to be projected and the center of the projection. The incidence part means that we know (at least) three collinear points to be projected. Combining this information, we can construct a surface of revolution containing the center of the projection. It is a generalized conic in the sense of [5, 10, 15, 16], because all of its points have a zero weighted distance sum from the elements of the (collinear) triplet of the projected images. Using iterative squares, such a generalized conic can also be represented as an algebraic surface. The rotational symmetry allows us to restrict the investigations to the defining polynomial of the profile curve in the image plane.

In this paper, a necessary and sufficient condition for the boundedness is given in terms of the input parameters, and it is shown that the defining polynomial of the profile curve is always irreducible in \(\mathbb {C}[x,y]\). Much of the calculation is computer-assisted, but rather surprisingly, after introducing a number of cases, the obstruction to reducibility always presents itself in the form of a simple contradiction. We refine the process by using more than three collinear points to substitute the generalized conics with spheres. Finally, the center of the projection can be given as the intersection of three spheres with non-collinear centers lying on the image plane.

Reconstructing the center of a projection is a fundamental concept in photogrammetry. Various methods of space resection are studied intensively based on different assumptions in the formulation of the problem [1, 4, 6,7,8]. The paper [6] provides a thorough review, by dividing these methods into three major categories: approximate solutions (including direct linear transformation, the Church method, and 3D conformal transformation), rigorous solutions based on fundamental conditions (such as using coplanarity conditions to determine the relative position of more cameras), and solutions based on projective geometry (for example using vanishing points for camera orientation and calibration). Our approach belongs to the second category, rigorous solutions based on fundamental conditions, where certain incidence relations are made use of to determine the parameters of one or more cameras. For instance, collinearity methods are based on the collinearity of the object, the center and the image of the object. In contrast, we derive equations for the coordinates of the center based on collinear objects in the space. There is a similar difference between coplanarity-based techniques in [6] and coplanarity conditions in Sect. 3.3.

We show that six collinear points in the space allow us to construct a sphere (instead of the more complicated generalized conics) containing the center of the projection. Such a sphere is centered at the image plane and the process ends by finding three spheres with non-collinear centers. Up to the symmetry about the image plane, their intersection determines the center of the projection. As soon as the center is found, we can reconstruct the coordinates of the original points in the space from the image coordinates.

This information can also be used for the follow-up editing of photographic images including, for example, the imitation of effects in the spatial area (fog, smoke, smog, brightness). Since the physical models use the coordinates of the original positions, we need a process to reconstruct the points in the space from the image coordinates. The results of the paper are partially motivated by a project with the involvement of both authors, where the goal is to simulate the effect of smog on real-life photographs. The case of homogeneous smog is well understood [9, 12,13,14]. In some databases, the pictures are equipped with a distance matrix, and that information is enough to simulate homogeneous fog [9]. However, for the inhomogeneous variant of the problem, these data are insufficient: one needs the original coordinates of each point projected on the photograph to create a realistic foggy image.

2 Reconstruction of the Center of a Projection: A Distance-and-Incidence Process

2.1 Notation and Terminology

Disregarding distortions and aberrations of optical systems [2, 11, 17], photographic images are produced by a central projection

$$\begin{aligned} \pi :M\subset \mathbb {R}^3\rightarrow F\subset \mathbb {R}^2 \end{aligned}$$

of a restricted area in the space into the image plane \(\mathbb {R}^2\), identified with the set \(\{(x,y,z)\in \mathbb {R}^3 \mid z=0\}\) if necessary. In practice, we have that \(\displaystyle {F:=[0,a]\times [0,b]}\) is an image rectangle constituted by a relatively dense grid of pixels. Throughout the paper, this simplified case of the central projection without any aberrations of the optical systems will be used (a.k.a. the pinhole mode).

Let \(\displaystyle {C(x_C, y_C, z_C)}\) be the center of the projection, \(\displaystyle {z_C>0}\) and suppose that \(\displaystyle {M\subset \pi ^{-1} (F)\subset \mathbb {R}^3}\) satisfying the following properties:

  • the center of the projection separates the point P and its projected image for all \(P\in M\) (which we indicate by \(\pi (P) - C - P\)), and

  • the distance \(r_{CP} = |C-P|\) is given for any \(P\in M\).

Fig. 1
figure 1

In our model, photographs are produced by central projections

In what follows, we use the coordinate system illustrated by Fig. 1. Since the following computations are independent of the vertical ordering, we omit the half-turn about the axis t of the camera. Identifying the points in the space with the position vectors with respect to the origin, the reconstruction formula of a point \(P\in M\) is

$$\begin{aligned} P=C+r_{CP}\frac{C-Q}{|C-Q|}, \end{aligned}$$
(1)

where \(Q=\pi (P)\). Its analytic form is

$$\begin{aligned} \begin{aligned}&P(x,y,z)=(x_C, y_C, z_C)\\&\quad +\frac{r_{CP}}{\sqrt{(x_C-x_Q)^2+(y_C-y_Q)^2+z_C^2}}\\&\quad (x_C-x_Q,y_C-y_Q,z_C). \end{aligned} \end{aligned}$$
(2)

In order to apply formula (2), we obviously need the coordinates of the center of the projection in terms of the (pixel) coordinates of the points in \(\displaystyle {\pi (M)\subset F}\) and the distances \(\displaystyle {r_{CP}}\) (\(P\in M\)).

To solve the reconstruction problem, we use the incidence relations among the points in the space. In practical applications, these relations can be detected by the help of images of special objects (facade, roadway, curb-stone, lamp-post, etc.).

2.2 Collinear Triples

Out of many possible incidence relations, we mainly focus on collinearity. Throughout this section, we adhere to the following convention.

Convention 1

Let \(P_1, P_2\) and \(P_3\) be three points projected to a photograph through center C with pairwise different, collinear images \(Q_1, Q_2\) and \(Q_3\), respectively. Assume that the order of this collinear triple is \(Q_1 - Q_3 - Q_2\), and that \(Q_3=(1-\lambda )Q_1+\lambda Q_2\) with some \(0< \lambda <1\). We put \(r_i = r_{CP_i}\) for \(i=1, 2, 3\). Furthermore, let

$$u=\frac{\lambda }{1-\lambda }, \,\,\text {and}\,\, a= \frac{1-\lambda }{r_1}, \,\, b= \frac{\lambda }{r_2}, \,\, c= \frac{1}{r_3}.$$

Note that these \(r_i\) are positive, and \(r_1=r_2=r_3\) is not possible, as three collinear points cannot lie on a sphere. Also note that \(u,a,b,c>0\).

Collinearity is preserved by central projections: thus, if \(P_1, P_2\ and\ P_3\) are collinear as in Fig. 2, then so are \(Q_1, Q_2\) and \(Q_3\), cf. Convention 1. The following result gives a necessary and sufficient condition for the reverse implication.

Fig. 2
figure 2

The images of collinear points under a central projection

Proposition 1

Under Convention 1, the points \(P_1\), \(P_2\) and \(P_3\) are collinear if and only if

$$\begin{aligned} a|C-Q_1|+ b|C-Q_2|- c|C-Q_3|=0. \end{aligned}$$
(3)

Proof

Using formula (1),

$$\begin{aligned}&(P_2-P_1)\times (P_3-P_1)\\&\quad =r_2r_3 \frac{(C-Q_2)\times (C-Q_3)}{|C-Q_2|\cdot |C-Q_3|}\\&\qquad -r_1r_2 \frac{(C-Q_2)\times (C-Q_1)}{|C-Q_2|\cdot |C-Q_1|}\\&\qquad -r_1r_3 \frac{(C-Q_1)\times (C-Q_3)}{|C-Q_1|\cdot |C-Q_3|}. \end{aligned}$$

Since \(Q_1\), \(Q_2\) and \(Q_3\) are collinear, the vector products are parallel vectors. Its analytic expression is based on the affine combination \(C-Q_3=(1-\lambda )(C-Q_1)+\lambda (C-Q_2)\). Hence, \((P_2-P_1)\times (P_3-P_1)\) is a scalar multiple of \((C-Q_2)\times (C-Q_1)\), where the scalar is

$$\begin{aligned}&\frac{(1-\lambda ) r_2r_3}{|C-Q_2|\cdot |C-Q_3|}-\frac{r_1r_2 }{|C-Q_2|\cdot |C-Q_1|}\\&\quad +\frac{\lambda r_1r_3}{|C-Q_1|\cdot |C-Q_3|}\\&\quad =\frac{\left( a|C-Q_1|+ b|C-Q_2|- c|C-Q_3|\right) r_1r_2r_3}{|C-Q_1|\cdot |C-Q_2|\cdot |C-Q_3|}. \end{aligned}$$

\(\square \)

Equation (3) provides a condition for the unknown center C in terms of the fixed quantities abc and the coordinates of \(Q_1, Q_2\) and \(Q_3\). It means the vanishing of a certain weighted distance sum of the point C from \(Q_1\), \(Q_2\) and \(Q_3\), respectively. It is a generalized conic [5]; see also [10, 15, 16]. Since the focal points are collinear, Eq. (3) is invariant under the rotation about the common line of \(Q_1\), \(Q_2\) and \(Q_3\). Therefore, it is a surface of revolution. This equation is of central importance for the rest of the paper.

2.3 A Sufficient and Necessary Condition for the Boundedness

We provide a simple equivalent condition for the boundedness of the generalized conic surface (3) in terms of the parameters in Convention 1.

Proposition 2

Under Convention 1, the generalized conic surface (3) is unbounded if and only if \(c=a+b\).

Proof

Suppose that there is an unbounded sequence \(C_n\) of points on the surface (3). Since \(C_n \ne Q_3\) for all but finitely many indices, we can write Eq. (3) into the form

$$\begin{aligned} c=a \frac{|C_n-Q_1|}{|C_n-Q_3|}+b \frac{|C_n-Q_2|}{|C_n-Q_3|} \end{aligned}$$
(4)

By using the triangle inequality, we obtain the upper estimates

$$\begin{aligned}&\frac{|C_n-Q_1|}{|C_n-Q_3|}\le \frac{|C_n-Q_3|+|Q_3-Q_1|}{|C_n-Q_3|}=1+\frac{|Q_3-Q_1|}{|C_n-Q_3|}, \\&\frac{|C_n-Q_1|}{|C_n-Q_3|}\ge \frac{|C_n-Q_1|}{|C_n-Q_1|+|Q_1-Q_3|}=\frac{1}{1+\frac{|Q_1-Q_3|}{|C_n-Q_1|}}. \end{aligned}$$

Therefore, we have \(\lim _{n\rightarrow \infty }\frac{|C_n-Q_1|}{|C_n-Q_3|}=1\), and similarly \(\lim _{n\rightarrow \infty }\frac{|C_n-Q_2|}{|C_n-Q_3|}=1\). Hence, the condition \(c=a+b\) is necessary.

For the reverse implication, let us restrict the investigations to the plane of the profile curve. If \(c=a+b\), then Eq. (3) can be written into the form

$$\begin{aligned} a \bigg (|C-Q_1|-|C-Q_3|\bigg ) = b \bigg (|C-Q_3|-|C-Q_2|\bigg ). \end{aligned}$$
(5)

In other words, there is a \(K\in \mathbb {R}\) such that

$$\begin{aligned}&|C-Q_1|-|C-Q_3|=\frac{K}{a}, \,\,\text {and} \end{aligned}$$
(6)
$$\begin{aligned}&|C-Q_3|-|C-Q_2|=\frac{K}{b}. \end{aligned}$$
(7)

The solutions of this system of equations are therefore the intersection points of one branch from each of two hyperbolae parameterized by K. If \(K\rightarrow 0\), then the branches tend to the perpendicular bisectors of the segments \(Q_1Q_3\) and \(Q_3Q_2\), respectively. These bisectors are distinct parallel lines since \(Q_3\) lies between \(Q_1\) and \(Q_2\). Hence, it is enough to show that some branches of these hyperbolae intersect for any \(K\ne 0\), as the set of common points of the branches cannot be contained in any bounded region of the plane. Therefore, the generalized conic cannot be bounded as was to be proved. Assuming that \(K>0\),

  • the real axes of the branches (6) and (7) coincide (the common line of the collinear points \(Q_1\), \(Q_2\) and \(Q_3\)); recall that the ordering is \(Q_1 - Q_3 - Q_2\).

  • The points of the branch (6) are closer to \(Q_3\) than to \(Q_1\).

  • The points of the branch (7) are closer to \(Q_2\) than to \(Q_3\).

In classical notation, the parameters of the hyperbolae depending on K are

$$\begin{aligned} a_1=\frac{K}{2a}, \ c_1=\frac{|Q_1-Q_3|}{2}, \ b_1=\sqrt{c_1^2-a_1^2} \end{aligned}$$

and

$$\begin{aligned} a_2=\frac{K}{2b}, \ c_2=\frac{|Q_2-Q_3|}{2}, \ b_2=\sqrt{c_2^2-a_2^2}. \end{aligned}$$

Here, \(a_1\) and \(a_2\) are the semimajor axes, \(c_1\) and \(c_2\) are the linear eccentricities, and \(b_1\) and \(b_2\) are the semiminor axes of the hyperbolae. The slopes of the asymptotic lines are \(m_1=\pm \frac{b_1}{a_1} \ \ \text {and}\ \ m_2=\pm \frac{b_2}{a_2},\) respectively. We show that these two slopes cannot be equal. According to the symmetry about the (common) real axis, it is enough to concentrate on the case of positive slopes. If \(m_1=m_2\), then \(\frac{b_1}{a_1}=\frac{b_2}{a_2}\) yields \(a^2\frac{|Q_1-Q_3|^2}{K^2}=b^2\frac{|Q_2-Q_3|^2}{K^2}\). By using the notation introduced in Convention 1, we have \(|Q_1-Q_3|^2 = \lambda ^2 |Q_2-Q_1|^2\) and \(|Q_2-Q_3|^2 = (1-\lambda )^2 |Q_2-Q_1|^2\); thus, we obtain

$$\begin{aligned} \frac{(1-\lambda )^2}{r_1^2}\frac{\lambda ^2 |Q_2-Q_1|^2}{K^2} = \frac{\lambda ^2}{r_2^2}\frac{(1-\lambda )^2 |Q_2-Q_1|^2}{K^2}, \end{aligned}$$

and consequently, \(r_1=r_2\). Then by \(c=a+b\) we obtain the contradiction \(r_1=r_2=r_3\), cf. Convention 1.

Hence, the slopes of the asymptotic lines must be different.

Fig. 3
figure 3

Illustration to the proof of Proposition 2

If \(m_1 < m_2\), then we have the intersection of the branches (6) and (7); see the branch with default linestyle in Fig. 3. If \(m_2 < m_1\), then we are looking for the intersection of the branches corresponding to \(K<0\); see the branch with dash linestyle in Fig. 3. \(\square \)

Example 1

In Fig. 4, an unbounded solution set of Eq. (3) with \(z=0\) is illustrated. The input data are

$$\begin{aligned} Q_1= & {} (1,2,0), \ r_1=14, \ Q_2=(1,3,0),\ r_2=21, \ \lambda =\frac{3}{10},\\ Q_3= & {} (1-\lambda )Q_1+\lambda Q_2=(1, \frac{3}{10}, 0), r_3= \frac{1}{\frac{1-\lambda }{r_1}+\frac{\lambda }{r_2}}=\frac{140}{9}. \end{aligned}$$
Fig. 4
figure 4

An unbounded profile curve in Example 1

2.4 Irreducibility of the Defining Polynomial in the Image Plane

Squaring iteratively, we can extend (3) to an algebraic surface consisting of the zeros of a polynomial of degree four:

$$\begin{aligned}&a |C-Q_1|+ b |C-Q_2|- c |C-Q_3| = 0\\&a |C-Q_1|+ b |C-Q_2| = c |C-Q_3|\\&c^2 |C-Q_3|^2- a^2 |C-Q_1|^2- b^2 |C-Q_2|^2 \\&\quad =2ab |C-Q_1|\cdot |C-Q_2| \end{aligned}$$

leading to the algebraic condition

$$\begin{aligned} \begin{aligned}&\left( c^2 |C-Q_3|^2- a^2 |C-Q_1|^2- b^2 |C-Q_2|^2\right) ^2\\&\quad -4 a^2b^2 |C-Q_1|^2\cdot |C-Q_2|^2=0. \end{aligned} \end{aligned}$$
(8)

Example 2

In case of the input data

$$\begin{aligned}&Q_1(1,2,0), \ Q_2(3,5,0), \ Q_3=(1-\lambda _3)Q_1+\lambda _3 Q_2,\\&\lambda _3=989\,{\frac{\sqrt{2}}{60\,\sqrt{23}+989\,\sqrt{2}}}, \ r_1=100,\\&r_2=215, \ r_3={\frac{10}{207}}\,\sqrt{9279189+543950\,\sqrt{23}\sqrt{2}} \end{aligned}$$

and

$$\begin{aligned}&Q_1(1,2,0), \ Q_2(3,5,0), \ Q_4=(1-\lambda _4)Q_1+\lambda _4 Q_2,\\&\lambda _4=989\,{\frac{\sqrt{2}}{240\,\sqrt{23}+989\,\sqrt{2}}}, \ r_1=100,\\&r_2=215, \ r_4={\frac{5}{207}}\,\sqrt{16420689+2175800\,\sqrt{23}\sqrt{2}}, \end{aligned}$$

the profile curves of the revolution surfaces are illustrated in Fig. 5.

Fig. 5
figure 5

Bounded profile curves

The input data in Example 2 have been obtained by maple-assisted computations (see Online Resource 1): setting the center of the projection to C(5, 4, 7), we

  • compute \(P_1\) and \(P_2\),

  • taking

    $$\begin{aligned} P_3=\frac{1}{3}P_1+\frac{2}{3}P_2 \end{aligned}$$

    compute its projected pair \(Q_3\),

  • compute \(\lambda _3\) and \(r_3\),

  • define the surface (3) for the plot command.

We investigate the profile curve of the surface given by (8) in the image plane, i.e., \(z=0\). The main result in this subsection is motivated by Bézout’s classical theorem [3]: by showing that the defining polynomial of the profile curve is always irreducible, we obtain that two such curves can intersect in at most 16 points.

Theorem 1

Under Convention 1, the polynomial obtained from the left-hand side of Eq. (8) by the substitution \(z=0\) is irreducible over \(\mathbb {C}\).

By applying an invertible affine transformation, if necessary, we may assume that the coordinates of the given points are \(Q_1(-u, 0)\), \(Q_2(1,0)\), \(Q_3(0,0)\); cf. Convention 1. Indeed, such a transformation preserves the irreducibility of polynomials.

The profile curve of the algebraic surface (8) is then defined by \(f(x,y)=0,\) where f(xy) is the quartic polynomial

$$\begin{aligned}&d_4(x^2+y^2)^2+d_3x(x^2+y^2)+d_{2,x}x^2\\&\quad +d_{2,y}y^2+d_1x+d_0 \end{aligned}$$

such that the coefficients are

$$\begin{aligned} \begin{aligned}&[(x^2+y^2)^2] : d_4\\&\quad =a^4 - 2a^2b^2 + b^4 - 2a^2c^2 - 2b^2c^2 + c^4, \\&[x(x^2+y^2)] : d_3\\&\quad =4(ua^4 + (1-u)a^2b^2 - ua^2c^2 - b^4 + b^2c^2), \\&[x^2] : d_{2,x}\\&\quad =2(3u^2a^4 + (-u^2+4u-1)a^2b^2\\&\qquad - u^2a^2c^2 + 3b^4 - b^2c^2),\\&[y^2] : d_{2,y}\\&\quad =2(u^2a^4 - (u^2+1)a^2b^2 - u^2a^2c^2 + b^4 - b^2c^2),\\&[x] : d_1=4(u^3a^4 + (u^2-u)a^2b^2 - b^4),\\&[1] : d_0=u^4a^4 - 2u^2a^2b^2 + b^4. \end{aligned} \end{aligned}$$

As we shall see, it is more convenient to switch the expression \(d_{2,x}\) in the above list to

$$\begin{aligned} d_{2,x}-d_{2,y}= 4(u^2a^4 + 2ua^2b^2 + b^4)= 4(ua^2+b^2)^2 \end{aligned}$$

in some cases. This leads to the observation

$$\begin{aligned} d_{2,x}-d_{2,y}\ne 0, \end{aligned}$$

which is frequently used in the following arguments. Under the comparison process, the system of equations generated by the original list of the coefficients is equivalent to the system of equations generated by the modified list containing \(d_{2,x}-d_{2,y}\) instead of \(d_{2,x}\). Finally, we note that

$$\begin{aligned} d_1=4(u^2a^2-b^2)(ua^2+b^2), \ d_0=(u^2a^2-b^2)^2 \end{aligned}$$

and

$$\begin{aligned} 4d_1^2-16(d_{2,x}-d_{2,y})d_0 = 0. \end{aligned}$$

In particular, \(d_1=0\) if and only if \(d_0=0\), i.e., \(b=ua\).

Lemma 1

The polynomial

$$\begin{aligned}&f(x,y)\\&\quad =d_4(x^2+y^2)^2+d_3x(x^2+y^2)\\&\quad +\,d_{2,x}x^2+d_{2,y}y^2+d_1x+d_0 \end{aligned}$$

has degree four or three in the variables x and y.

Proof

The leading coefficient \(d_4\) can be factorized:

$$\begin{aligned}&a^4+b^4+c^4-2a^2b^2-2b^2c^2-2c^2a^2\\&\quad = -(a+b+c)(a+b-c)(a-b+c)(-a+b+c). \end{aligned}$$

As ab and c are positive, \(c=-a-b\) is impossible. Hence, the degree of the polynomial f(xy) is less than four if and only if \(c=a+b\), \(c=a-b\) or \(c= b-a\). Since every power of c is even in the coefficients of f(xy), we obtain the same expressions by substituting \(c=a-b\) and \(c=b-a\). In this case, the cubic coefficient \(d_3\) simplifies to \(8ab(a-b)(ua+b)\), which cannot be zero as \(c\ne 0\). On the other hand, if f(xy) is not quartic and \(c=a+b\), then the coefficients of f(xy) can be written as

$$\begin{aligned} \begin{aligned}&[x(x^2+y^2)] : d_3 =8ab(-ua+b)(a+b),\\&[x^2] : d_{2,x}\\&\quad =4(ua-b)(ua^3+(1-u)a^2b+(1-u)ab^2-b^3),\\&[y^2] : d_{2,y}\\&\quad =-4ab(u^2a+b)(a+b),\\&[x] : d_1\\&\quad =4(u^2a^2-b^2)(ua^2+b^2),\\&[1] : d_0\\&\quad =(u^2a^2-b^2)^2. \end{aligned} \end{aligned}$$

Hence, the cubic coefficient \(d_3\) vanishes if and only if \(b=ua\). However, by the definition of u this yields \(r_1=r_2\), and then by \(c=a+b\) we obtain \(r_1=r_2=r_3\), a contradiction (cf. Convention 1). \(\square \)

We break the analysis down to six cases based on the next two lemmas.

Lemma 2

Assume that \(d_4\ne 0\) and f(xy) is reducible in \(\mathbb {C}[x,y]\). By applying the invertible linear substitution \(x\mapsto x\), \(y\mapsto -y\) if necessary, the polynomial f can be written as a product of not necessarily irreducible complex polynomials in at least one of the following four ways:

  1. 1.

    \(f(x,y)= d_4(x^2+2ixy-y^2+ Ax + iAy +C)\cdot \) \((x^2-2ixy-y^2+ Ax - iAy +C), \text {or}\)

  2. 2.

    \(f(x,y)= d_4(x^2+y^2+ Ax + C)(x^2+y^2+ Dx + F)\), or

  3. 3.

    \(f(x,y)= d_4(x^2+y^2+ Ax + By + C)\cdot \) \((x^2+y^2+ Ax - By + C), \text {or}\)

  4. 4.

    \(f(x,y)= d_4(x^3 - ix^2y + xy^2 - iy^3 + Ax^2 +\) \(i(C-A)xy + Cy^2 + Dx + i(AC-C^2-D)y+\) \(C(C^2-AC+D))(x + iy + C),\)

where A, B, C, D, \(F\in \mathbb {C}\).

Proof

By our assumption, the polynomial

$$\begin{aligned}&\frac{1}{d_4}f(x,y)\\&\quad = (x^2+y^2)^2+\frac{d_3}{d_4}x(x^2+y^2)+\frac{d_{2,x}}{d_4}x^2\\&\qquad +\frac{d_{2,y}}{d_4}y^2+\frac{d_1}{d_4}x+\frac{d_0}{d_4} \end{aligned}$$

is reducible. The sum of all degree four terms in the variables x and y is \((x^2+y^2)^2\), which is obtained as the product of the highest degree monomials in the factors. In the first case suppose that both factors are quadratic. Up to multiplication by nonzero constants, there are two ways to write \((x^2+y^2)^2\) as a product of quadratic polynomials, namely \((x+iy)^2\cdot (x-iy)^2\) and \((x^2+y^2)\cdot (x^2+y^2)\). If the degree two parts of the factors are \((x+iy)^2\) and \((x-iy)^2\), then

$$\begin{aligned} \frac{1}{d_4}f(x,y)= & {} (x^2+2ixy-y^2+ Ax+ By +C)(x^2 \\&-2ixy-y^2+ Dx + Ey +F). \end{aligned}$$

By \([y^3]=0\), we have \(E=-B\), and \([x^2y]=0\) (together with \(E=-B\)) yields \(D=A\). Comparing the coefficients of \(x^3\) and \(xy^2\), it follows that \(B=iA\). Finally, as \([xy] = 2i(F-C) =0\) we have \(C=F\) and item (1). If the degree two parts of the factors are both \((x^2+y^2)\), then

$$\begin{aligned} \frac{1}{d_4}f(x,y)= & {} (x^2+y^2+ Ax + By +C)(x^2\\&+y^2+ Dx + Ey +F). \end{aligned}$$

Since \([x^2y]=0\), we have \(E=-B\). Therefore, \([xy]=0\) and \([y]=0\) yield \(B(D-A)=0\) and \(B(F-C)=0\). Finally, item (2) and item (3) follow immediately depending on \(B=0\) or \(B\ne 0\).

In the second case suppose that there is a cubic and a linear factor. Since we are allowed to apply the linear substitution \(x\mapsto x\), \(y\mapsto -y\), there are no essentially different ways to write the degree four part as a product of a cubic and a linear expressions, except \((x^2+y^2)^2 = ((x-iy)^2(x+iy))\cdot (x+iy)\). Therefore,

$$\begin{aligned}&\frac{1}{d_4}f(x,y) =(x^3 - ix^2y + xy^2 - iy^3 \\&\quad + Ax^2 + Bxy + Cy^2 + Dx + Ey + F)\cdot \\&\quad (x + iy + G). \end{aligned}$$

Since \([y^3]=0\), we have \(C=G\). Hence, \([x^2y]=0\) yields \(B=i(C-A)\). The equation \([xy]=0\) implies that \(E=i(AC-C^2-D)\) and, by \([y]=0\), we have \(F=iEC = C(C^2-AC+D)\) as was to be proved in item (4). \(\square \)

Lemma 3

Assume that \(d_4 = 0\), \(d_3\ne 0\), and f(xy) is reducible in \(\mathbb {C}[x,y]\). By applying the invertible linear substitution \(x\mapsto x\), \(y\mapsto -y\) if necessary, the polynomial f can be written as a product of not necessarily irreducible complex polynomials in at least one of the following two ways:

  1. (5)

    \(f(x,y) = d_3(x^2+y^2+ Ax +C)(x + D)\), or

  2. (6)

    \(f(x,y)= d_3(x^2+ixy+ Ax + By +B(B-iA))\) \((x -iy + (A + iB)),\)

where A, B, C, \(D\in \mathbb {C}\).

Proof

By our assumption, the polynomial

$$\begin{aligned}&\frac{1}{d_3}f(x,y)\\&\quad = x(x^2+y^2)+\frac{d_{2,x}}{d_3}x^2+\frac{d_{2,y}}{d_3}y^2+\frac{d_1}{d_3}x+\frac{d_0}{d_3} \end{aligned}$$

is reducible. The sum of all degree three terms in x and y is \(x(x^2+y^2)\), which is obtained as the product of the highest degree monomials in the factors. Since one of the factors is quadratic, the other is linear, and we are allowed to apply the linear substitution \(x\mapsto x\), \(y\mapsto -y\), there are essentially two ways to write the degree three part as such a product, namely \((x^2+y^2)\cdot x\) and \((x(x+iy))\cdot (x-iy)\). In the first case

$$\begin{aligned} \frac{1}{d_3}f(x,y) = (x^2+y^2+ Ax + By +C)(x + D) \end{aligned}$$

and equation \([xy] = 0\) yields \(B=0\). In the second case

$$\begin{aligned} \frac{1}{d_3}f(x,y) = (x^2+ixy+ Ax + By +C)(x -iy + D) \end{aligned}$$

and equation \([xy] = 0\) yields \(D=A + iB\). Finally, we obtain that \(C= B(B-iA)\) because of \([y]=0\). \(\square \)

We are ready to prove that the defining polynomial is always irreducible. The complete calculation run by sagemath is provided in Online Resource 2.

Proof of Theorem 1

According to Lemma 1, we need to show that if the degree of f(xy) is four or three, then f(xy) is irreducible in \(\mathbb {C}[x,y]\). The proof follows the cases (1)–(6) based on Lemmas 2 and 3.

  1. (1)

    By comparison of coefficients,

    $$\begin{aligned} d_3/d_4= & {} 2A, \ d_{2,x}/d_4 = A^2 + 2C, \ d_{2,y}/d_4 = A^2 - 2C,\\ d_1/d_4= & {} 2AC, \ d_0/d_4 = C^2. \end{aligned}$$

Therefore, \(2(d_{2,x} + d_{2,y})d_0 - d_1^2=0\). A sagemath-assisted computation shows that

$$\begin{aligned}&2(d_{2,x} + d_{2,y})d_0 - d_1^2\\&\quad = -8(u^2a^2-b^2)^2((u^2a^2+b^2)c^2+(u+1)^2a^2b^2). \end{aligned}$$

Its vanishing is equivalent to \(b=ua\), i.e., \(d_0=0\), \(C=0\) and \(d_{2,x}=d_{2,y}\). This is a contradiction because of

$$\begin{aligned} d_{2,x}-d_{2,y} = 4(u^2a^4 + 2ua^2b^2 + b^4)= 4(ua^2+b^2)^2 >0. \end{aligned}$$
  1. (2)

    By comparison of coefficients,

    $$\begin{aligned} d_3/d_4= & {} A+D, \ d_{2,x}/d_4 = AD + C + F,\\ d_{2,y}/d_4= & {} C + F, \ d_1/d_4 = CD + AF,\\ d_0/d_4= & {} CF. \end{aligned}$$

In particular,

$$\begin{aligned} A+D=d_3/d_4 \ \ \text {and}\ \ AD = \frac{d_{2,x}-d_{2,y}}{d_4}. \end{aligned}$$

Therefore, A and D are the solutions of the quadratic equation

$$\begin{aligned} z^2-\frac{d_3}{d_4}z+\frac{d_{2,x}-d_{2,y}}{d_4} = 0, \end{aligned}$$

i.e.,

$$\begin{aligned} A,D = \frac{d_3\pm \sqrt{d_3^2-4(d_{2,x}-d_{2,y})d_4}}{2d_4}. \end{aligned}$$

Similarly, C and F are the solutions of the quadratic equation

$$\begin{aligned} z^2-\frac{d_{2,y}}{d_4}z+\frac{d_0}{d_4} = 0, \end{aligned}$$

i.e.,

$$\begin{aligned} C,F = \frac{d_{2,y}\pm \sqrt{d_{2,y}^2-4d_0d_4}}{2d_4}. \end{aligned}$$

We obtain that

$$\begin{aligned}&2d_4^2(CD+AF)-d_3d_{2,y}\\&\quad = \pm \sqrt{d_3^2-4(d_{2,x}-d_{2,y})d_4}\sqrt{d_{2,y}^2-4d_0d_4}. \end{aligned}$$

Putting \(CD+AF = d_1/d_4\) and taking the square of both sides, we have that the expression

$$\begin{aligned}&d_1^2d_4 - d_1d_3d_{2,y} + (d_{2,x}-d_{2,y})d_{2,y}^2\\&\quad + d_0d_3^2 - 4d_0(d_{2,x}-d_{2,y})d_4 \end{aligned}$$

vanishes. Using a sagemath-assisted computation, we get that this expression is

$$\begin{aligned} 64(u^4+2u^3+u^2)a^4b^4c^4>0, \end{aligned}$$

which is a contradiction.

  1. (3)

    By comparison of coefficients,

    $$\begin{aligned} d_3/d_4= & {} 2A, \ d_{2,x}/d_4 = \ A^2 + 2C\\ d_{2,y}/d_4= & {} -B^2+2C, \ d_1/d_4 = 2AC, \ d_0/d_4 = C^2. \end{aligned}$$

Clearly, \(d_1^2d_4 - d_3^2d_0 = 0\). Using a sagemath-assisted computation,

$$\begin{aligned} d_1^2d_4 - d_3^2d_0= & {} 64a^2b^2c^2(u^2a^2 - b^2)^2(uc^2\\&-((u^2+u)a^2 + (u+1)b^2)). \end{aligned}$$

Therefore, \(b=ua\), or \(c^2=(u+1)\left( a^2 + \frac{1}{u}b^2\right) \). We show that the condition \(c^2=(u+1)\left( a^2 + \frac{1}{u}b^2\right) \) also implies the equation \(b=ua\). First of all observe that

$$\begin{aligned} (4d_{2,x}d_4 - d_3^2)^2 - 64d_0d_4^3 = 0. \end{aligned}$$

The left-hand side is a polynomial expression containing the even powers of c. By applying the above condition \(c^2=(u+1)\left( a^2 + \frac{1}{u}b^2\right) \), a sagemath-assisted computation shows that

$$\begin{aligned}&(4d_{2,x}d_4 - d_3^2)^2 - 64d_0d_4^3\\&\quad = 256a^2b^2(u+1)\left( a^2 + \frac{1}{u}b^2\right) (u^2a^2-b^2)^4\cdot \\&\qquad \left( \left( 1 + \frac{1}{u}\right) a^2 + \left( \frac{1}{u^3} + \frac{1}{u^4}\right) b^2\right) \end{aligned}$$

and its vanishing implies that \(b=ua\). To complete the argument, observe that \(d_0=(u^2a^2-b^2)^2=0\), i.e., \(C=0\) and \(d_3^2-4d_{2,x}d_4=0\). If \(b=ua\), then

$$\begin{aligned} d_3^2-4d_{2,x}d_4=16u^2a^2c^2(c^2-(u+1)^2a^2)^2 \end{aligned}$$

and, consequently, \(c=(u+1)a=a+b\). In particular, \(d_4=0\), which is a contradiction.

  1. (4)

    By comparison of coefficients,

    $$\begin{aligned} d_3/d_4= & {} A+C, \ d_{2,x}/d_4 = AC + D,\\ d_{2,y}/d_4= & {} -AC+2C^2+D, \ d_1/d_4\\= & {} -AC^2+C^3+2CD,\\ d_0/d_4= & {} -AC^3 + C^4 + C^2D. \end{aligned}$$

Using the identity \((2AC - 2C^2)C^2 + 2(AC^2 - C^3 - 2CD)C + 4(-AC^3 + C^4 + C^2D)=0\), we have that C is a root of the quadratic polynomial \((d_{2,x}-d_{2,y})z^2 - 2d_1z+4d_0\); recall that \(d_{2,x}-d_{2,y}\ne 0\).

According to the vanishing discriminant \(4d_1^2-16(d_{2,x}-d_{2,y})d_0 = 0\),

$$\begin{aligned} C = \frac{d_1}{d_{2,x}-d_{2,y}}. \end{aligned}$$

Using that \(A= \frac{d_3}{d_4}-C\) and \(2C(A-C)=\frac{d_{2,x}-d_{2,y}}{d_4}\), we have

$$\begin{aligned} 4C^2 - \frac{2d_3}{d_4}C + \frac{d_{2,x}-d_{2,y}}{d_4} =0, \end{aligned}$$

that is, \(4d_4C^2 - 2d_3C + (d_{2,x}-d_{2,y}) =0\). Putting \(C = \frac{d_1}{d_{2,x}-d_{2,y}}\) yields

$$\begin{aligned} 4d_4d_1^2 - 2d_3d_1(d_{2,x}-d_{2,y}) + (d_{2,x}-d_{2,y})^3 =0. \end{aligned}$$

Computing the left-hand side expression with sagemath, we obtain a quadratic polynomial in \(c^2\) such that the coefficients are real polynomials in uab. Since \(c^2\in \mathbb {R}\) is a real root, the discriminant must be a non-negative real number. In a more detailed form, the discriminant is

$$\begin{aligned} -65536(u^3+2u^2+u)a^4b^4(b^2-u^2a^2)^2(ua^2+b^2)^4 \end{aligned}$$

as a sagemath-assisted computation shows. Therefore, \(b=ua\), \(d_0=d_1=0\) and \(C^2(-AC+C^2+D)=C(-AC+C^2+2D)=0\). If \(C\ne 0\) then \(-AC+C^2+D = -AC+C^2+2D=0\), i.e., \(D=0\) and \(A=C\). Another consequence is that

$$\begin{aligned} \frac{d_{2,x}}{d_4}=\frac{d_{2,y}}{d_4}=A^2, \end{aligned}$$

and \(d_{2,x}-d_{2,y}=0\), which is a contradiction. The same argument is working in case of \(C=0\).

  1. (5)

    By comparison of coefficients,

    $$\begin{aligned}&d_{2,x}/d_3 = A + D, \ d_{2,y}/d_3 = D,\\&d_1/d_3 = AD + C, \ d_0/d_3 = CD. \end{aligned}$$

The identity \(CD = ((AD+C)-D((A+D)-D))D\) yields

$$\begin{aligned} \frac{d_0}{d_3} = \left( \frac{d_1}{d_3}-\frac{d_{2,y}}{d_3}\frac{d_{2,x}-d_{2,y}}{d_3}\right) \frac{d_{2,y}}{d_3}. \end{aligned}$$

Hence, \((d_1d_3-d_{2,y}(d_{2,x}-d_{2,y}))d_{2,y} - d_0d_3^2 =0\). However,

$$\begin{aligned}&(d_1d_3-d_{2,y}(d_{2,x}-d_{2,y}))d_{2,y} - d_0d_3^2\\&\quad = -64(u^4 + 2u^3 + u^2)a^4b^4c^4<0 \end{aligned}$$

as a sagemath-assisted computation shows.

  1. (6)

    By comparison of coefficients,

    $$\begin{aligned} d_{2,x}/d_3= & {} 2A + iB, \ d_{2,y}/d_3 = -iB, \ d_1/d_3 = A^2 + B^2\\ d_0/d_3= & {} -iA^2B + 2AB^2 + iB^3 \end{aligned}$$

We can express the variables A and B from the first and the second equations:

$$\begin{aligned} A=\frac{d_{2,x}+d_{2,x}}{2d_3}, \ B=\frac{id_{2,y}}{d_3}. \end{aligned}$$

By substituting them into the third and the fourth equations, we obtain that

$$\begin{aligned} g:= & {} d_{2,x}^2 + 2d_{2,x}d_{2,y} - 3d_{2,y}^2 - 4d_1d_3=0,\\ h:= & {} 4d_0d_3^2 - (d_{2,x} + d_{2,y})^2d_{2,y} + 4(d_{2,x}\\&+ d_{2,y})d_{2,y}^2 - 4d_{2,y}^3=0. \end{aligned}$$

Using sagemath,

$$\begin{aligned}&h-4ab(u^2a^2+(u^2+1)ab+b^2)g \\&\quad = 256u^2(u + 1)^2a^4b^4(a+b)^4 \,\, \text {if}\,\, c=a+b, \, \text {and}\\&h+4ab(u^2a^2-(u^2+1)ab+b^2)g = \\&\quad 256u^2(u + 1)^2a^4b^4(a-b)^4 \,\, \text {if}\,\, c=a-b \,\, \text {or} \,\, c=b-a \end{aligned}$$

The contradiction is obvious in the case that \(c=a+b\). Finally, if \(c=a-b\) or \(c=b-a\), then \(a-b=0\) and, consequently, \(c=0\), which is a contradiction. \(\square \)

The case \(c=a-b\) is illustrated by the following input data:

$$\begin{aligned}&Q_1(-1/9,0,0), \ Q_2(1,0,0), \ Q_3(0,0,0),\\&\quad r_1=1/60, \ r_2=1/40, \ r_3=1/50, \end{aligned}$$

i.e., \(\lambda =1/10\) and \(a=54\), \(b=4\), \(c=50\), \(u=1/9\). Figure 6 shows both the zeros of the polynomial f(xy) (left) and the profile curve of the generalized conic in the image plane (right).

Fig. 6
figure 6

The case \(c=a-b\)

3 Reconstruction of the Center from More than Three Collinear Points

3.1 Collinear Quadruples

Consider a collinear quadruple of the points \(P_1\), \(P_2\), \(P_3\) and \(P_4\) in the space together with the collinear quadruple of the points \(Q_1\), \(Q_2\), \(Q_3\) and \(Q_4\) on the image plane. Introducing the notation

$$\begin{aligned} r_i:=r_{CP_i} \quad (i=1, 2, 3, 4) \end{aligned}$$

for the distance parameters, let us investigate the system of equations

$$\begin{aligned} \begin{aligned} \frac{1-\lambda _3}{r_1} |C-Q_1|+\frac{\lambda _3}{r_2} |C-Q_2|-\frac{1}{r_3} |C-Q_3|=0,\\ \frac{1-\lambda _4}{r_1} |C-Q_1|+\frac{\lambda _4}{r_2} |C-Q_2|-\frac{1}{r_4} |C-Q_4|=0, \end{aligned} \end{aligned}$$
(9)

where

$$\begin{aligned} Q_3=(1-\lambda _3)Q_1+\lambda _3 Q_2, \ Q_4=(1-\lambda _4)Q_1+\lambda _4 Q_2. \end{aligned}$$

By increasing the number of the input data points, we are going to substitute the generalized conics with simpler revolution surfaces containing the center of the projection. We have

$$\begin{aligned} \begin{aligned} \frac{1-\lambda _3}{r_1} |C-Q_1|+\frac{\lambda _3}{r_2} |C-Q_2|=\frac{1}{r_3} |C-Q_3|,\\ \frac{1-\lambda _4}{r_1} |C-Q_1|+\frac{\lambda _4}{r_2} |C-Q_2|=\frac{1}{r_4} |C-Q_4|. \end{aligned} \end{aligned}$$
(10)

Taking the square of both sides, we find

$$\begin{aligned}&2 |C-Q_1|\cdot |C-Q_2|=\frac{r_1 r_2}{\lambda _3(1-\lambda _3)}\cdot \\&\quad \left( \frac{1}{r_3^2}|C-Q_3|^2-\frac{(1-\lambda _3)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _3^2}{r_2^2} |C-Q_2|^2\right) ,\\&2 |C-Q_1|\cdot |C-Q_2|=\frac{r_1 r_2}{\lambda _4(1-\lambda _4)}\cdot \\&\quad \left( \frac{1}{r_4^2} |C-Q_4|^2-\frac{(1-\lambda _4)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _4^2}{r_2^2} |C-Q_2|^2\right) \end{aligned}$$

and, consequently,

$$\begin{aligned} \begin{aligned}&\frac{r_1 r_2}{\lambda _3(1-\lambda _3)}\cdot \\&\left( \frac{1}{r_3^2}|C-Q_3|^2-\frac{(1-\lambda _3)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _3^2}{r_2^2} |C-Q_2|^2\right) \\&=\frac{r_1 r_2}{\lambda _4(1-\lambda _4)}\cdot \\&\left( \frac{1}{r_4^2} |C-Q_4|^2-\frac{(1-\lambda _4)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _4^2}{r_2^2} |C-Q_2|^2\right) . \end{aligned} \end{aligned}$$
(11)

Equation (11) defines an algebraic revolution surface in the space containing the center of the projection and it is of degree at most two.

Theorem 2

The profile curve of the revolution surface given by equation

$$\begin{aligned} \begin{aligned}&\frac{1}{\lambda _3(1-\lambda _3)}\cdot \\&\left( \frac{1}{r_3^2}|C-Q_3|^2-\frac{(1-\lambda _3)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _3^2}{r_2^2} |C-Q_2|^2\right) \\&=\frac{1}{\lambda _4(1-\lambda _4)}\cdot \\&\left( \frac{1}{r_4^2} |C-Q_4|^2-\frac{(1-\lambda _4)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _4^2}{r_2^2} |C-Q_2|^2\right) \end{aligned} \end{aligned}$$
(12)

for \(Q_3\ne Q_4\) is either a circle such that the center lies on the common line of \(Q_1\), \(Q_2\), \(Q_3\) and \(Q_4\) or the profile curve reduces to a line perpendicular to the common line of the points \(Q_1\), \(Q_2\), \(Q_3\) and \(Q_4\).

Proof

First of all, we show that the profile curve is non-trivial. Using a similarity transformation if necessary, we can suppose that \(Q_1(0,0)\), \(Q_2(1,0)\), \(Q_3(\lambda _3,0)\) and \(Q(\lambda _4,0)\). Taking C(xy), Eq. (12) yields

$$\begin{aligned} \begin{aligned} \frac{1}{\lambda _3(1-\lambda _3)}&\biggr (\frac{1}{r_3^2}((x-\lambda _3)^2+y^2)-\\&\frac{(1-\lambda _3)^2}{r_1^2}(x^2+y^2)-\frac{\lambda _3^2}{r_2^2} ((x-1)^2+y^2)\biggr )\\ =\frac{1}{\lambda _4(1-\lambda _4)}&\biggr (\frac{1}{r_4^2}((x-\lambda _4)^2+y^2)-\\&\frac{(1-\lambda _4)^2}{r_1^2}(x^2+y^2)-\frac{\lambda _4^2}{r_2^2} ((x-1)^2+y^2)\biggr ). \end{aligned} \end{aligned}$$

We collect the coefficients of \((x^2+y^2)\), x and the constant term on both sides of the equation. These coefficients coincide if and only if

$$\begin{aligned}&\frac{1}{\lambda _3(1-\lambda _3)}\frac{1}{r_3^2} -\frac{1-\lambda _3}{\lambda _3}\frac{1}{r_1^2} -\frac{\lambda _3}{1-\lambda _3}\frac{1}{r_2^2}\nonumber \\&\quad =\frac{1}{\lambda _4(1-\lambda _4)}\frac{1}{r_4^2}- \frac{1-\lambda _4}{\lambda _4}\frac{1}{r_1^2} -\frac{\lambda _4}{1-\lambda _4}\frac{1}{r_2^2}, \end{aligned}$$
(13)
$$\begin{aligned}&\frac{1}{1-\lambda _3}\frac{1}{r_3^2} -\frac{\lambda _3}{1-\lambda _3}\frac{1}{r_2^2}\nonumber \\&\quad =\frac{1}{1-\lambda _4}\frac{1}{r_4^2}-\frac{\lambda _4}{1-\lambda _4}\frac{1}{r_2^2}, \text {and} \end{aligned}$$
(14)
$$\begin{aligned}&\frac{\lambda _3}{1-\lambda _3}\frac{1}{r_3^2} -\frac{\lambda _3}{1-\lambda _3}\frac{1}{r_2^2}\nonumber \\&\quad =\frac{\lambda _4}{1-\lambda _4}\frac{1}{r_4^2}-\frac{\lambda _4}{1-\lambda _4}\frac{1}{r_2^2}. \end{aligned}$$
(15)

Subtracting Eqs. (14) and (15) yields

$$\begin{aligned} \frac{1-\lambda _3}{1-\lambda _3}\frac{1}{r_3^2} = \frac{1-\lambda _4}{1-\lambda _4}\frac{1}{r_4^2}, \end{aligned}$$

and it follows that \(r_3=r_4\). Therefore, Eq. (14) implies that \(r_2=r_3=r_4\) or \(\lambda _3=\lambda _4\), i.e., \(Q_3=Q_4\). Since both cases are impossible due to the collinearity of \(P_2\), \(P_3\) and \(P_4\), we have that the profile curve is an algebraic curve of degree one or two. More precisely, if Eq. (13) holds then it is a line of the form

$$\begin{aligned} Bx+D=0. \end{aligned}$$

Otherwise, it is a circle of the form

$$\begin{aligned} A(x^2+y^2)+Bx+D=0 \end{aligned}$$

as was to be proved. \(\square \)

Remark 1 In the sense of Theorem 2, the set of common points of the profile curves of the conics

$$\begin{aligned} \frac{1-\lambda _3}{r_1} |C-Q_1|+\frac{\lambda _3}{r_2} |C-Q_2|-\frac{1}{r_3} |C-Q_3|=0 \end{aligned}$$

and

$$\begin{aligned} \frac{1-\lambda _4}{r_1} |C-Q_1|+\frac{\lambda _4}{r_2} |C-Q_2|-\frac{1}{r_4} |C-Q_4|=0 \end{aligned}$$

must be on the same line or a circle. Substituting the rational parameterization x(t), y(t) and \(z(t)=0\) of the intersection points into the algebraic expression (8) of one of the conics, we have only finitely many solutions according to Theorem 1.

3.2 Finding Generic Collinear Quadruples

In what follows, we prove a theorem to provide the circle as the locus of the solutions of system (9) in the image plane. This means that the center of the projection is on a sphere centered at the common line of the points \(Q_1\), \(Q_2\), \(Q_3\) and \(Q_4\).

Theorem 3

If the points \(P_1\), \(P_2\), \(P_3\), \(P_4\), \(P_5\), \(P_6\) are collinear then we can find a quadruple \((Q_1, Q_2, Q_i, Q_j)\) among the image points such that the profile curve of the revolution surface

$$\begin{aligned} \begin{aligned}&\frac{1}{\lambda _i(1-\lambda _i)}\cdot \\&\left( \frac{1}{r_i^2}|C-Q_i|^2-\frac{(1-\lambda _i)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _i^2}{r_2^2} |C-Q_2|^2\right) \\&=\frac{1}{\lambda _j(1-\lambda _j)}\cdot \\&\left( \frac{1}{r_j^2} |C-Q_j|^2-\frac{(1-\lambda _j)^2}{r_1^2} |C-Q_1|^2-\frac{\lambda _j^2}{r_2^2} |C-Q_2|^2\right) \end{aligned} \end{aligned}$$
(16)

is a circle, where \(i, j=3, 4, 5, 6\) and \(i\ne j\).

Proof

To prove the theorem, consider the left-hand side of (13) as a function of \(\lambda _3\) in the first step. Since \(r_3\) also depends on \(\lambda _3\), we have to clarify the correspondence \(r_3(\lambda _3)\) before investigating the extended expression. Without loss of generality, suppose that

$$\begin{aligned} Q_1 - Q_3 - Q_2, \ Q_1(0,0), \ Q_2(1,0) \ \text {and}\ Q_3(\lambda _3,0). \end{aligned}$$

Under the notation

$$\begin{aligned} \alpha =\angle Q_1 Q_3 C, \quad c_i=|C-Q_i| \quad (i=1, 2, 3), \end{aligned}$$

let us apply the cosine rule to the triangles \(\triangle Q_1Q_3C\) and \(\triangle Q_2Q_3C\), respectively. It follows that

$$\begin{aligned} c_3^2+\lambda _3^2-2\lambda _3 c_3\cos \alpha =c_1^2 \end{aligned}$$

and

$$\begin{aligned} c_3^2+(1-\lambda _3)^2+2(1-\lambda _3) c_3\cos \alpha =c_2^2. \end{aligned}$$

Therefore,

$$\begin{aligned} (1-\lambda _3)\left( c_3^2+\lambda _3^2-2\lambda _3 c_3\cos \alpha \right) =(1-\lambda _3)c_1^2 \end{aligned}$$

and

$$\begin{aligned} \lambda _3 \left( c_3^2+(1-\lambda _3)^2+2(1-\lambda _3) c_3\cos \alpha \right) =\lambda _3c_2^2. \end{aligned}$$

Taking the sum of the equations, one finds

$$\begin{aligned} c_3^2=(1-\lambda _3)c_1^2+\lambda _3 c_2^2-\lambda _3(1-\lambda _3) \end{aligned}$$

and, by Eq. (3),

$$\begin{aligned} \frac{1}{r_3^2}= & {} \frac{1}{(1-\lambda _3)c_1^2+\lambda _3 c_2^2-\lambda _3(1-\lambda _3)}\cdot \\&\left( (1-\lambda _3)^2\frac{c_1^2}{r_1^2} +\lambda _3^2\frac{c_2^2}{r_2^2}+2(1-\lambda _3)\lambda _3\frac{c_1 c_2}{r_1 r_2}\right) . \end{aligned}$$

Introducing the function

$$\begin{aligned} \begin{aligned}&V(t):= \frac{1}{(1-t)c_1^2+t c_2^2-t(1-t)} \cdot \\&\left( \frac{1-t}{t}\frac{c_1^2}{r_1^2} +\frac{t}{1-t}\frac{c_2^2}{r_2^2}+2\frac{c_1c_2}{r_1 r_2} \right) -\frac{1-t}{t}\frac{1}{r_1^2} -\frac{t}{1-t}\frac{1}{r_2^2}, \end{aligned} \end{aligned}$$
(17)

the profile curve is not a circle if and only if Eq. (13) holds, i.e.,

$$\begin{aligned} V(\lambda _3)=V(\lambda _4). \end{aligned}$$

Using Rolle’s theorem, there must be a stationary point satisfying \(V'(t)=0\) between the parameters. In the second step, we are going to prove that there are at most two stationary points of the function (17). The maple-assisted computation

figure c

gives the roots

$$\begin{aligned}&{\frac{{{c_1}}^{2}{r_1}+{c_1}{c_2}{r_1}- {c_1}{c_2}{r_2}-{{c_2}}^{2}{r_2}+{ r_2}}{{r_1}+{r_2}}},\\&\quad {\frac{{{c_1}}^{2}{r_1}-{c_1}{c_2}{r_1}- {c_1}{c_2}{r_2}+{{c_2}}^{2}{r_2}-{r_2}}{{r_1}-{r_2}}}. \end{aligned}$$

Therefore, one of the equations

$$\begin{aligned} V(\lambda _3)=V(\lambda _4), \ V(\lambda _4)=V(\lambda _5), \ V(\lambda _5)=V(\lambda _6) \end{aligned}$$

cannot be true. \(\square \)

3.2.1 The Steps of the Reconstructing Process

Input data: The image points \(Q_1, \ldots , Q_6\) of the collinear points \(P_1, \ldots , P_6\in M\), and \(r_i=r_{CP_i} \,\, \text {for}\,\, i=1 \ldots , 6\). Suppose that the indices follow the ordering \(Q_1 - Q_3 - Q_4 - Q_5 - Q_6 - Q_2\).

  1. Step 1

    Consider the quadruples \((Q_1, Q_ 2, Q_ i, Q_j)\) (\(i\ne j\) and \(i, j=3, 4, 5, 6\)) in all possible ways to find at least one such that the profile curve of the revolution surface (given by (16)) is a circle (see Theorem 3). Without loss of generality, we may assume that the quadruple \((Q_1, Q_2, Q_3, Q_4)\) is one of them. The output is the center O of the profile circle.

Up to a similarity transformation, we can suppose that \(Q_1(0,0)\), \(Q_2(1,0)\), \(Q_3(\lambda _3,0)\) and \(Q(\lambda _4,0)\) and the equation of the circle is of type

$$\begin{aligned} A(x^2+y^2)+Bx+D=0, \end{aligned}$$

where

$$\begin{aligned} A= & {} \frac{1}{\lambda _3(1-\lambda _3)}\frac{1}{r_3^2} -\frac{1-\lambda _3}{\lambda _3}\frac{1}{r_1^2} -\frac{\lambda _3}{1-\lambda _3}\frac{1}{r_2^2}-\\&-\left( \frac{1}{\lambda _4(1-\lambda _4)}\frac{1}{r_4^2}- \frac{1-\lambda _4}{\lambda _4}\frac{1}{r_1^2} -\frac{\lambda _4}{1-\lambda _4}\frac{1}{r_2^2}\right) ,\\&-\frac{1}{2}B=\frac{1}{1-\lambda _3}\frac{1}{r_3^2} -\frac{\lambda _3}{1-\lambda _3}\frac{1}{r_2^2}-\\&\left( \frac{1}{1-\lambda _4}\frac{1}{r_4^2}-\frac{\lambda _4}{1-\lambda _4}\frac{1}{r_2^2}\right) , \text {and}\\ D= & {} \frac{\lambda _3}{1-\lambda _3}\frac{1}{r_3^2} -\frac{\lambda _3}{1-\lambda _3}\frac{1}{r_2^2}-\left( \frac{1}{1-\lambda _4}\frac{1}{r_4^2}-\frac{\lambda _4}{1-\lambda _4}\frac{1}{r_2^2}\right) , \end{aligned}$$

because of (13), (14) and (15). The coordinates of the center of the circle in the image plane are \(O\left( -\frac{B}{2A},0\right) .\)

  1. Step 2

    Iterate Step 1 to two other 6-tuples to obtain two more spheres with centers \(O'\) and \(O''\).

  2. Step 3

    Finally, compute the intersection of these three spheres with coplanar centers.

If the points O, \(O'\) and \(O''\) are not collinear then, up to a reflection about the image plane, the center of the projection is uniquely determined. In practice, the points O, \(O'\) and \(O''\) are almost never collinear. If the input data provide us with collinear centers then the set of possible solutions is a circle.

3.3 A Note About Coplanar Quadruples in the Space

The coplanarity of the points \(P_1\), \(P_2\), \(P_3\) and \(P_4\) can be detected by the help of images of special objects (facade, roadway, traffic signs etc.). The sufficient and necessary condition can be given as the vanishing of the triple product as follows:

$$\begin{aligned} 0=\langle (P_2-P_1)\times (P_3-P_1), P_4-P_1\rangle . \end{aligned}$$
(18)

Using formula (1),

$$\begin{aligned}&(P_2-P_1)\times (P_3-P_1)\\&\quad = r_2r_3 \frac{(C-Q_2)\times (C-Q_3)}{|C-Q_2|\cdot |C-Q_3|}\\&\qquad -r_1r_2\frac{(C-Q_2)\times (C-Q_1)}{|C-Q_2|\cdot |C-Q_1|}-\\&\quad r_1r_3 \frac{(C-Q_1)\times (C-Q_3)}{|C-Q_1|\cdot |C-Q_3|},\\&P_4-P_1=r_4\frac{C-Q_4}{|C-Q_4|}-r_1\frac{C-Q_1}{|C-Q_1|}, \end{aligned}$$

where \(Q_i=\pi (P_i)\) and \(r_i=r_{CP_i}\), for \(i=1, \ldots , 4\). Therefore, we can formulate Eq. (18) in terms of the projected points:

$$\begin{aligned} 0= & {} r_2 r_3 r_4\frac{\langle (C-Q_2)\times (C-Q_3), C-Q_4\rangle }{|C-Q_2|\cdot |C-Q_3|\cdot |C-Q_4|} -\\&\quad r_1 r_2 r_4\frac{\langle (C-Q_2)\times (C-Q_1), C-Q_4\rangle }{|C-Q_2|\cdot |C-Q_1|\cdot |C-Q_4|} -\\&\quad r_1 r_3 r_4\frac{\langle (C-Q_1)\times (C-Q_3), C-Q_4\rangle }{|C-Q_1|\cdot |C-Q_3|\cdot |C-Q_4|} -\\&\quad r_1 r_2 r_3\frac{\langle (C-Q_2)\times (C-Q_3), C-Q_1\rangle }{|C-Q_2|\cdot |C-Q_3|\cdot |C-Q_1|}. \end{aligned}$$

We also have that

$$\begin{aligned}&\langle (C-Q_2)\times (C-Q_3), C-Q_4\rangle :\\&\quad \text {the distance of { C} from the image plane}=-\frac{\mu _{234}}{3}, \end{aligned}$$

where

$$\begin{aligned} \mu _{234}=\left\{ \begin{array}{rl} > 0&{}\text {if the triangle} \triangle Q_2 Q_3 Q_4 \text {is positively}\\ &{} \text {oriented with respect to} C \\ < 0&{}\text {otherwise } \end{array} \right. \end{aligned}$$

is the signed area of the triangle \(\triangle Q_2 Q_3 Q_ 4\). Therefore, Eq. (18) is equivalent to

$$\begin{aligned} 0= & {} -\mu _{234} \frac{r_2 r_3 r_4 }{|C-Q_2|\cdot |C-Q_3|\cdot |C-Q_4|} +\\&\quad \mu _{214}\frac{r_1 r_2 r_4}{|C-Q_2|\cdot |C-Q_1|\cdot |C-Q_4|}+\\&\quad \mu _{134}\frac{r_1 r_3 r_4}{|C-Q_1|\cdot |C-Q_3|\cdot |C-Q_4|}+\\&\quad \mu _{231}\frac{r_1 r_2 r_3 }{|C-Q_2|\cdot |C-Q_3|\cdot |C-Q_1|} \end{aligned}$$

and, consequently, the center of the projection must be on the generalized conic surface given by

$$\begin{aligned} \begin{aligned} 0=W_{234} |C-Q_1| +W_{134} |C-Q_2|+\\ W_{124} |C-Q_3|+W_{123} |C-Q_4|, \end{aligned} \end{aligned}$$
(19)

where the weights are given by

$$\begin{aligned} W_{234}= & {} -\mu _{234} r_2 r_3 r_4, \ W_{134}= \mu _{134} r_1 r_3 r_4,\\ W_{124}= & {} - \mu _{124} r_1 r_2 r_4, \ W_{123}= \mu _{123} r_1 r_2 r_3. \end{aligned}$$

4 Conclusions

Our explicit methods presented in Sect. 3 can complement existing computational techniques; see the introduction and [6] for a survey. The distances from the center and the collinearity of the points in the space turned out to be very effective input data for the reconstruction of the center of a projection. It is demonstrated that a small number of collinear triplets determine the center via relatively simple equations. We applied classical geometrical and algebraic techniques, such as (generalized) conics and Bézout’s theorem, together with computer-assisted calculations in maple and sagemath. The theoretical results guarantee rigorous solutions based on fundamental conditions. Some perspectives are presented in Sect. 3.3 for future research to establish new combined techniques.