














Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
The radiosity equation, a method used in computer graphics to calculate the energy transfer and reflection between surfaces in a scene. The equation involves the total energy passing between points, energy emitted directly, reflection coefficients, and visibility terms. The document also covers the relationship between the rendering equation and radiosity equation, and provides formulas for finding the form factor and solving for radiosity values.
Typology: Study notes
1 / 22
This page cannot be seen from the preview
Don't miss anything!















Lecture 23: Radiosity
And he shall be as the light of the morning
1. Radiosity
Radiosity models the transfer of light between surfaces. Recursive ray tracing also models
light bouncing off surfaces, but ray tracing makes several simplifying assumptions that give scenes
a harsh, unnatural look. Ray tracing assumes that all light sources are point sources and that the
ambient light is constant throughout the scene. These assumptions lead to stark images with sharp
shadows. In contrast, radiosity models all surfaces as both emitters and reflectors. This approach
softens the shadows and provides a more realistic model for ambient light.
Informally, radiosity is the rate at which light energy leaves a surface. There are two
contributions to radiosity: emission and reflection. Hence
radiosity = emitted energy + reflected energy.
For the purpose of display, we shall identify intensity with radiosity.
Radiosity computations typically take much longer than recursive ray tracing because the
model of light is much more complex. To simplify these computations, we shall model only diffuse
reflections; we shall not attempt to model specular reflections with radiosity. Since radiosity
replaces ambient and diffuse intensity, radiosity is view independent. Thus we can reuse the same
computation for every viewpoint once we compute the radiosity of all the surfaces. View dependent
calculations are required only to compute hidden surfaces.
2. The Radiosity Equations
We will begin with a very general integral equation called the Rendering Equation based on
energy conservation. We will then repeatedly simplify and discretize this equation till we get a large
system of linear equations -- the Radiosity Equations -- which we can solve numerically for the
radiosity of each surface. Radiosity is identified with intensity, so once we have the radiosity for
each surface we can render the scene.
2.1 The Rendering Equation. Energy conservation for light is equivalent to
Total Illumination = Emitted Energy + Reflected Energy.
We can rewrite these innocent looking words as an integral equation called the Rendering Equation.
Rendering Equation
I ( x ,
x ) = E ( x ,
x ) + ρ( x ,
x ,
x )
S
x ,
x ) d
x (2.1)
where
I ( x ,
x ) is the total energy passing from
x to
x.
E ( x ,
x ) is the energy emitted directly from
x to
x.
ρ( x ,
x ,
x ) is the reflection coefficient -- the percentage of the energy transferred from
x to
x that is passed on to
x.
Essentially all of the computations in Computer Graphics that involve light are summarized in
the Rendering Equation. Notice, in particular, that the Rendering Equation is precisely the setup for
recursive ray tracing!
2.2 The Radiosity Equation -- Continuous Form. The continuous form of the Radiosity
Equation is just the Rendering Equation restricted to diffuse reflections. Once again by
conservation of energy:
Radiosity = Emitted Energy + Reflected Energy.
Now, however, since we are dealing only with diffuse reflections, we can be more specific about the
form of the reflected energy. Restricting to diffuse reflections leads to the following integral
equation for radiosity.
Radiosity Equation -- Continuous Form
B ( x ) = E ( x ) + ρ
d
( x ) B ( y )
S
cos θ cos θ′
π r
2
( x , y )
V ( x , y ) dy (2.2)
where
B ( x ) is the radiosity at the point x , which we identify with the intensity or energy
reflecting off a surface in any direction -- that is, the total power leaving a
surface/unit area/solid angle. This energy is uniform in all directions, since we
are assuming that the scene has only diffuse reflectors.
E ( x ) is the energy emitted directly from a point x in any direction. This energy is
uniform in all directions, since we are assuming that the scene has only diffuse
emitters.
ρ
d
( x ) is the diffuse reflection coefficient -- the percentage of energy reflected in all
directions from the surface at a point x. By definition,
0 ≤ ρ
d
( x ) ≤ 1.
V ( x , y ) is the visibility term:
V ( x , y ) = 0 if x is not visible from y.
V ( x , y ) = 1 if x is visible from y.
θ = angle between the surface normal ( N ) at x and the light ray ( L ) to y.
θ = angle between surface normal (
N ) at y and the light ray ( L ) to x.
r ( x , y ) = distance from x to y.
ρ(
x , x ,
x ) = the percentage of the energy transferred from
x to
x that is passed on to a
single point
x. (Note that here we have reversed the roles of x and
x .)
whereas in the Radiosity Equation
ρ
d
( x ) = the percentage of energy reflected in all directions from the surface at a point x.
Thus we need to get from
ρ(
x , x ,
x ) in the Rendering Equation to
ρ
d
( x ) in the Radiosity
Equation. We shall now show that these two functions differ by a factor of π.
Directions can be represented by points on the unit sphere, and points on the unit sphere can
be parametrized by two angles
θ,φ. The angle
θ represents the angle between the z -axis and the
vector from the center of the sphere to the point
( θ, φ) on the sphere, and the angle
φ represents the
amount of rotation around the equator from the x -axis to the great circle of constant longitude
φ.
Since
ρ
d
( x ) is the percentage of energy reflected in all direction whereas
ρ(
x , x ,
x ) is the
percentage of energy reflected in one fixed direction, it follows that
ρ
d
( x ) is the integral of
ρ(
x , x ,
x ) over the unit hemisphere -- that is, over all the directions that make an acute angle with
the normal vector at x, which we identify with the z -axis. A factor of
cos( θ) must appear in this
integral for the same reason (projection) that this factor appears in Lambert’s Law. Thus
ρ
d
( x ) = ρ(
x , x ,
x )
H
cos( θ) dS ,
where H is a unit hemisphere centered at x. Since we are dealing with diffuse reflection, the
function
ρ(
x , x ,
x ) is the same in all directions. Thus we can pull
ρ(
x , x ,
x ) outside the integral,
so
ρ
d
( x ) = ρ(
x , x ,
x ) cos( θ) dS
H
It remains then to compute
cos( θ) dS
H
To compute this integral, we use the parametrization
( θ, φ) of the unit sphere provided by
spherical coordinates -- that is, by setting
S ( θ,φ) = sin( θ)cos( φ),sin(θ)sin(φ),cos( θ) ( )
With this parametrization, a differential area (parallelogram) on the unit sphere is given by
dS =
∂ S ( θ, φ)
∂ θ
d θ ×
∂ S ( θ,φ)
∂ φ
d φ =
∂ S (θ, φ)
∂ θ
∂ S ( θ,φ)
∂ φ
d θ d φ.
But by Equation (2.4)
∂ S ( θ, φ)
∂ θ
= cos( θ)cos( φ), cos( θ)sin(φ),− sin( θ) ( )
∂ S ( θ, φ)
∂ φ
= − sin( θ)sin(φ),sin(θ)cos( φ), 0 ( )
so by direct computation
∂ S ( θ,φ)
∂ θ
∂ S ( θ,φ)
∂ φ
= sin(θ).
Therefore
cos( θ) dS =
π / 2
∫
2 π
∫
cos( θ)sin( θ) d θ d
H
sin
( θ)
2 π
∫
π /
π /
∫
We conclude then from Equation (2.3) that
ρ
d
( x ) = π ρ(
x , x ,
x ).
Thus when we replace
ρ(
x , x ,
x ) in the Rendering Equation by
ρ
d
( x ) in the Radiosity Equation,
we must divide by π. This accounts for the factor π in the denominator of the second term on the
right hand side of the Radiosity Equation.
2.3 The Radiosity Equations -- Discrete Form. To find the radiosity at any point, we must
solve the Radiosity Equation for
B ( x ) -- that is, we must calculate the integral on the right hand side
of Equation (2.2). But this integral is not easy to compute. Worse yet
B ( u ) appears on both sides
of the Radiosity Equation; thus we must know
B ( u ) to find
B ( u )!
The solution to both problems is to discretize the Radiosity Equation by breaking the surfaces
in the scene into small patches. Since radiosity is approximately constant over a small patch, we
shall be able to replace integrals by sums and continuous functions by discrete values. The integral
equation then reduces to a large system of linear equations which we shall be able to solve by
numerical methods.
To discretize the Radiosity Equation, we begin by breaking the surfaces S into small patches
j
j = 1 ,K, N. Over a small patch
j
, the radiosity
B ( y ) is approximately a constant
j
, and the
integral over S breaks up into a sum of integrals over the patches
j
. Therefore by Equation (2.2)
B ( x ) = E ( x ) + ρ
d
( x ) B ( y )
S
cos θ cos θ′
π r
2
( x , y )
V ( x , y ) dy
≈ E ( x ) + ρ
d
( x ) B
j
j = 1
N
cos θ cos θ′
π r
2
( x , y )
j
V ( x , y ) dA
j
We still need to discretize
B ( x ) on the left hand side and
E ( x ) on the right hand side of
Equation (2.5). We can approximate the radiosity and energy of single patch
i
by an area
weighted average. Thus
i
i
) B ( x ) dA
i
P
i
i
i
) E ( x ) dA
i
i
3. Form Factors
To gain a better understanding of the form factors, we are going to provide both a physical and
a geometric interpretation for these constants. We begin with a physical interpretation.
Theorem 1: The form factor
ij
is the fraction of the energy leaving patch
i
arriving at patch
j
Proof: Let
ij
denote the radiosity transferred from
i
to
j
. From the Radiosity Equation,
ij
ji
i
Since radiosity is energy per unit area, to find the total energy transferred from
i
to
j
, we must
multiply the radiosity transferred from
i
to
j
by the area of
j
. Let
i
be the area of patch
i
and let
j
be the area of patch
j
. Then the total energy transferred from
i
to
j
is
j
ij
j
ji
i
But by Equation (2.6)
j
ji
cos θ cos θ′
π r
2
( x , y )
j
V ( x , y ) dA
j
dA
i P
i
i
ij
Therefore
j
ij
i
i
ij
so
ij
j
ij
i
i
Since the radiosity
i
is the energy radiated in all directions from
i
per unit area, the denominator
i
i
represents the total energy radiated from
i
. But we have already seen that the numerator
j
ij
represents the total energy transferred from
i
to
j
. Therefore, it follows from Equation
(3.1) that the form factor
ij
is the fraction of the energy that leaves
i
and arrives at
j
Corollary 1: For each value of i, the form factors
ij
form a partition of unity. That is,
ij
j
Proof: This result follows from conservation of energy. By Theorem 1 the form factor
ij
is the
fraction of the energy leaving patch
i
arriving at patch
j
. Since energy is conserved, the energy
leaving patch
i
must arrive somewhere -- that is, at one of the patches
j
. Therefore
ij
j
To provide a geometric interpretation for the form factors, we must first discretize still further.
Over a small patch
i
, we can treat the inner integral as roughly constant, so
ij
i
cos θ cos θ′
π r
2
( x , y )
j
V ( x , y ) dA
j
dA
i P
i
cos θ cos θ′
π r
2
( x , y )
V ( x , y )
j
dA
j
We shall now investigate the meaning of the differential
cos θ cos θ′
π r
2
( x , y )
dA
j
The product
cos θ′
r
2
dA
j
is the projection of the differential area
dA
j
onto the unit hemisphere
centered at the patch
dA
i
. Similarly, the product
(cos θ)
cos θ′
r
2
dA
j
is the projection of the
differential area
cos θ′
r
2
dA
j
onto the base of the hemisphere centered at
dA
i
, the plane perpendicular
to the normal of
dA
i
(see Figure 2).
i
dA
j
Figure 2: Projecting a patch
j
first onto the unit hemisphere centered at the patch
dA
i
and then
onto the base of the hemisphere centered at
dA
i
In 3-dimensions the cosine terms in these projections are not so easy to visualize, so to get a
feel for these cosine terms, let us look instead in 2-dimensions. To project a differential length
dL
onto a circle of radius r, where r is the distance between the patches, we multiply
dL by
cos(
θ )
where
θ is the angle between the tangent to the circle of radius r and the tangent to the curve
dL
(see Figure 3a). Of course, the angle between the tangents is the same as the angle between the
normals. If we think of
dL as one patch and if the second patch is located at the center of the circle
Corollary 2: Two surfaces
j
j ′
with the same projection onto the unit hemisphere centered at
a small patch
dA
i
have the same form factor -- that is,
ij
i j ′
Corollary 2 is the main result of all this analysis: two patches with equal projection on the unit
hemisphere centered around a small patch
dA
i
will necessarily have the same form factor
ij
. We
can use this insight to compute the form factors once and for all for some simple surface and then
find the form factors for arbitrary surfaces by projecting onto the known surface.
3.1 Hemi–Cubes. A hemi-cube is the upper half of a cube with sides of length two, centered at a
small patch
i
(see Figure 4). Here we are going to compute explicit formulas for the form factors
for small patches on the surface of the hemi-cube. To find the form factors for arbitrary patches,
we will project these patches onto the hemi-cube and use Corollary 2 from the previous section
which says that patches with equal projections have equal form factors.
Figure 4: A hemi-cube. The lengths of the sides are twice the height of the hemi-cube.
To compute form factors for small patches on the surface of the hemi-cube, recall from
Equation (3.3) that in general the form factor is given by
di , j
cos θ cos θ′
π r
2
( x , y )
V ( x , y )
j
dA
j
For a very small patch with area
j
, the integrand can be approximated by a constant, so
di , dj
cos θ cos θ′
π r
2
j
To find an explicit formula for the form factor
di , dj
, we need to find explicit formulas for
cos θ ,
cos
θ , and r , where
θ is the angle between the normal N to the patch
i
and the vector from the center of
i
to
the center of
j
θ is the angle between the normal
N to the patch
j
and the vector from the center of
i
to the center of
j
i
to the center of
j
For the hemi-cube there are two kinds of patches to consider: patches on the top of the hemi-
cube and patches along the sides of the hemi-cube. To compute the form factors for these patches,
choose a coordinate system with the origin at the center of the hemi-cube and the z -axis aligned
with the normal to the patch
i
at the center of the hemi-cube. For a small patch on the top face of
the hemi-cube centered at the point
P = ( x , y ,1), it follows from simple trigonometry (see Figure 5a)
that
cos θ = cos
θ = 1 / r.
Therefore for small patches with area
j
on the top face of the hemi-cube centered at the point
P = ( x , y ,1),
di , dj
cos θ cos θ′
π r
2
j
j
π r
4
j
π ( x
2
2
2
€
N ′
€
r = x
2
2
€
( x , y ,1)
€
€
€
θ
€
N ′
€
r = 1 + y
2
2
€
( 1 , y , z )
€
€
€
θ
€
1
€
θ ′
z
(a)
cos θ = cos
θ = 1 / r (b)
cos θ = z / r cos
θ = 1 / r
Figure 5: Schematic views of the top and side faces of a hemi-cube. (a) For the top face of the
hemi-cube, the vectors
N are parallel. Therefore by simple trigonometry,
cos θ = cos
θ = 1 / r.
(b) For the side faces of the hemi-cube, the vectors
N are orthogonal. Therefore by simple
trigonometry,
cos θ = z / r and
cos
θ = 1 / r.
Similarly, for a small patch on the side face of the hemi-cube parallel to the yz -plane centered at
the point
P = (± 1 , y , z ) , it again follows by simple trigonometry (see Figure 5b) that
cos θ = z / r cos
θ = 1 / r.
Therefore for small patches on side faces of the hemi-cube parallel to the yz -plane
di , dj
cos θ cos θ′
π r
2
j
z Δ A
j
π r
4
z Δ A
j
π ( y
2
2
2
Finally, for a small patch on the side face of the hemi-cube parallel to the xz -plane centered at the
point
P = ( x , ± 1 , z ) , it follows by an analogous argument that
cos θ = z / r cos
θ = 1 / r
4. The Radiosity Rendering Algorithm
We introduced radiosity in order to develop more realistic looking images using Computer
Graphics. Now let us put together what we now know about the Radiosity Equations, form factors,
hemi-cubes, and shading to develop a rendering algorithm based on radiosity.
Radiosity Rendering Algorithm
Break each surface into small surface patches.
Use the hemi-cube algorithm.
i
i
i
ij
j
j = 1
N
i = 1 ,K, N.
(See below.)
Use any hidden surface algorithm.
The only steps that require further elaboration are step 3 and step 4. In step 3 we must solve a
large system of linear equations. For large linear systems, standard techniques like Gaussian
elimination are slow and unstable. We shall provide instead two alternative robust numerical
methods for solving these equations, but we defer this discussion till the next section. In step 4 we
need to find the radiosity for each vertex so that we can perform Gouraud shading in step 7 in order
to eliminate discontinuities in intensity between adjacent patches.
For regular meshes consisting of rectangular patches, there are three kinds of vertices: interior,
boundary, and corner (see Figure 6). For interior vertices, it is natural simply to set the intensity to
the average of the radiosities at the four adjacent patches:
interior
1
2
3
4
For edges and corners there are two competing strategies: either set
edge
1
2
corner
1
or set
edge
interior
1
2
corner
interior
1
so that
edge
1
2
interior
1
2
3
4
corner
1
interior
1
2
3
4
Both strategies yield reasonable results for Gouraud shading.
€
B
1
€
B
2
€
B
3
€
B
4
€
I
€ interior
I
edge
€
I
corner
€
€
€
Figure 6: The intensities at the vertices depend on the radiosities of the adjacent patches. There
are three kinds of vertices: interior vertices, edge vertices, and corner vertices. Each vertex type has
a slightly different formula for intensity based on the radiosities of the adjacent patches.
5. Solving the Radiosity Equations
The Radiosity Equations are:
i
i
i
ij
j
j = 1
N
i = 1 ,K, N.
Bringing all the radiosities to the left hand side yields
( δ
ij
− ρ
i
ij
j
j = 1
N
i
i = 1 ,K, N ,
or equivalently
( 1 − ρ
1
11
1
− ρ
1
12
2
− L− ρ
1
1 N
N
1
− ρ
2
21
1
2
22
2
−L − ρ
2
2 N
N
2
− ρ
N
N 1
1
− ρ
N
N 2
2
−L + ( 1 − ρ
N
N N
N
N
Jacobi Relaxation
i
p
i
ii
ij
ii
j ≠ i
∑
j
p − 1
i = 1 ,K, N
Gauss-Seidel Relaxation
i
p
i
ii
ij
ii
j = 1
i − 1
∑
j
p
ij
ii
j = i + 1
N
∑
j
p − 1
i = 1 ,K, N
In Jacobi relaxation we use the values of the guess at the previous level to generate the values
of the guess at the next level. In Gauss-Seidel relaxation, we use the values of the guess at the
previous level along with the values of the guess already computed at the current level to compute
the next value of the guess at the current level. Gauss-Seidel relaxation is a bit more complicated
than Jacobi relaxation, but Gauss-Seidel relaxation typically converges faster to the solution of the
linear system. (For additional details, see Lecture 7, Section 3.2.)
Convergence is guaranteed in both relaxation methods for any initial guess when the system is
diagonally dominant. A system of equations such as (5.1) is diagonally dominant if
ii
ij
j ≠ i
It follows easily from Corollary 1 that the Radiosity Equations -- Equations (2.7) -- are diagonally
dominant (see Exercise 1), so we can use these relaxation methods to solve the Radiosity Equations.
Solving for the radiosity using relaxation methods is called gathering because we gather the
radiosity from all the surfaces simultaneously. The disadvantage of gathering is that we must
compute all the form factors to get an intermediate result. Thus we must solve for all the hemi-
cubes radiosities before we can begin to render the scene. Typically solving for all the hemi-cube
radiosities and computing all the form factors takes a long time, so if there is some error in the
scene or in the code we will have to wait a long time to detect the mistake. Therefore we shall
consider another technique called shooting , where we can render intermediate results progressively
without waiting to compute all the form factors and all the radiosities for each hemi-cube.
5.2 Shooting -- Progressive Refinement. Each patch contributes to the radiosity of every other
patch. To compute radiosity progressively, we fix one particular patch and compute the radiosity of
every other patch due to the radiosity of the fixed patch. We can then display the scene and repeat
the process by choosing another patch. In this way we get to see intermediate results quickly
without waiting to finish the entire computation.
From the Radiosity Equations
i
i
i
ij
j
j = 1
N
i = 1 ,K, N ,
so the radiosity
i
due to
j
is
ρ
i
ij
j
. To solve the Radiosity Equations by gathering, we need
to compute all the form factors for every patch; thus we need one hemi-cube for each patch to
initiate the computation.
But in shooting, we are interested initially only in the radiosity
j
due to one fixed radiosity
i
. The radiosity
j
due to
i
is
ρ
j
ji
i
. To compute the form factors
ji
directly for each j ,
we would still need to introduce one hemi-cube computation for each patch. But recall that by the
reciprocity relationships (Equation (3.4))
i
ij
j
ji
therefore
j
due to B
i
= ρ
j
ji
i
= ρ
j
ij
i
i
j
Now to compute
ij
for each j , we need the form factors
ij
only for the single patch
i
in order to
update all the patches
j
! Thus we need to introduce only one hemi-cube at a time to update all the
patch radiosities at once. This reduction can save a lot of time.
In shooting, we think of radiosity as accumulating on each patch until we shoot this radiosity
out to all the other patches. Let
j
denote the unshot radiosity of the patch
j
. Then from the
Radiosity Equations it follows that after we shoot the accumulated radiosity
i
from the patch
i
j
j
old
j
ij
i
i
j
j
j
old
j
ij
i
i
j
i
These equations lead to the following progressive refinement procedure.
Shooting Algorithm
0
1
2
N
0
6. Summary
Radiosity is a classical energy transfer technique adapted to rendering surfaces in Computer
Graphics. The advantages of radiosity over other standard rendering methods such as recursive ray
tracing are that radiosity provides better photorealistic effects such as softer shadows and color
bleeding. Radiosity computations are also view independent, since typically radiosity models only
ambient light and diffuse reflection.
The main disadvantages of radiosity are that radiosity computations are expensive both in time
and in space. Good accuracy demands lots of small patches to model curved surfaces, leading to
lots of form factors and lots of hemi-cubes as well as a very large system of linear equations for the
radiosities. In addition radiosity does not typically model specular reflections, since specular
reflections are view dependent.
Enhancements to standard radiosity include a two pass algorithm for computing specular
reflections, more realistic light sources, participating mediums such as atmospheric fog, adaptive
mesh generation for computing more accurate form factors, and finite element methods or wavelets
for calculating more accurate approximations to radiosity. The interested reader can find these
subjects in the literature; we shall not pursue these topics here.
Below for easy reference we review the Rendering Equation along with the continuous and
discrete forms of the Radiosity Equation and the form factors.
Rendering Equation
I ( x ,
x ) = E ( x ,
x ) + ρ( x ,
x ,
x )
S
x ,
x ) d
x (2.1)
where
I ( x ,
x ) is the total energy passing from
x to
x.
E ( x ,
x ) is the energy emitted directly from
x to
x.
ρ( x ,
x ,
x ) is the reflection coefficient -- the percentage of the energy transferred from
x to
x that is passed on to
x.
Radiosity Equation -- Continuous Form
B ( x ) = E ( x ) + ρ
d
( x ) B ( y )
cos θ cos θ′
π r
2
( x , y )
V ( x , y ) dy (2.2)
where
B ( x ) is the radiosity at the point x , which we identify with the intensity or energy -- that is,
the total power leaving a surface/unit area/solid angle.
E ( x ) is the energy emitted directly from a point x. This energy is uniform in all
directions, since we are assuming that the scene has only diffuse emitters.
ρ
d
( x ) is the diffuse reflection coefficient --
0 ≤ ρ
d
( x ) ≤ 1.
V ( x , y ) is the visibility term:
V ( x , y ) = 0 if x is not visible from y.
V ( x , y ) = 1 if x is visible from y.
θ = angle between the surface normal ( N ) at x and the light ray ( L ) to y.
θ = angle between surface normal (
N ) at y and the light ray ( L ) to x.
r ( x , y ) = distance from x to y.
Radiosity Equations -- Discrete Form
i
i
i
ij
j
j = 1
N
i = 1 ,K, N (2.7)
where
i
is the radiosity of patch
i
, which we identify with intensity.
i
is the energy emitted from patch
i
(uniform in all directions, since by assumption we
are dealing only with diffuse emitters).
ij
are the form factors, which depend only on the geometry of the scene and are
independent of the lighting.
ρ
i
is the diffuse reflection coefficient for patch
i
0 ≤ ρ
i
Form Factors
ij
i
cos θ cos θ′
π r
2
( x , y )
j
V ( x , y ) dA
j
dA
i P
i
Exercises:
i , j ∪ k
i , j
i , k
ij
1
N
( )
T
1
N
( )
T
, and let D be the diagonal matrix defined
by
ij
ii
i = j
= 0 i ≠ j.
Consider the system of linear equations: