Use dual basis as a filter Projection

Ψ ijkl (u , v , s , t )= R 2 cos 2 φ ϕ ij (u , v)ϕ kl (s , t ) 〈 ϕ mn ' | ϕ ij 〉 =δ m i δ n j E ijkl = L ijkl L (u , v , s , t )= ∑ ijkl L ijkl Ψ ijkl (u , v , s , t ) Quadratic Basis Quadratic Basis Light Source L ( x , y , θ , ϕ ) ϕ ij ( u , v ) ϕ kl ( s , t ) R ( u , v , s , t ) Filter Support easure/Image Support φ ( u , v ,
Ψ ijkl a (u , v , s , t )=ϕ ij (u , v )ϕ kl (s , t ) E ijkl = ( cos 2 φ R 2 ) ijkl ∫ uvst L (u , v , s , t )ϕ ij ' (u , v )ϕ kl ' (s , t ) E ' ijkl = ∫ uvst L ( u , v , s , t ) ϕ ij ' ( u , v ) ϕ kl ' ( s , t ) L (u , v , s , t )= ∑ ij ϕ ij (u , v) ∑ kl E ' ijkl ϕ kl (s , t ) L (u , v , s , t )= ∑ ij ϕ ij (u , v)T ij (s , t )
Lighting Approximation Lighting Approximation

Constant value over a basis support for :

Geometric term Visibility BRDF

Average energy through a basis

Similar to a mip-map

Texturing + Multi-Pass

I ( x , o )= ∑ ij Γ ij ( x , o , i )V ij ( x )G ij ( x )T ij ( x )
Hardware Implementation (1/2) Hardware Implementation (1/2)

For each measure Shadow map (ie Visibility) Lighting

2N rendering passes Requirement

Floating-point textures Floating-point buffer for the sum

V ij ( x ) Γ ij ( x , o , i )G ij ( x )T ij ( x )
Hardware Implementation (2/2) Hardware Implementation (2/2) Decompression

e b =log 2 ( max (| R f | , | G f | , | B f |) ) +2 [ R ,G , B ]=2 -e b [ R f , G f , B f ] +0.5 E = e+126 255 [ R f , G f , B f ] =2 255⋅E -126 ([ R , G , B ] -0.5 )
Φ Φ ijkl (u , v , s , t ) =Φ ϕ ij (u ,v )⋅ϕ kl ( s ,t ) =Φ ϕ i , bellq (u )⋅ϕ j , bellq (v )⋅ϕ k , const ( s)⋅ϕ l , const (t )

  s , t )

	Approximations Approximations Shift-Invariant Approximation Shift-Invariant Approximation
	Faster rendering Constant geometric term (error ≤ 8%)
	Texture based
	Shift Invariant
	Remove geometric term in the reconstruction function
	Direct lighting New reconstruction function
	Use each measure as a texture
	One lighting evaluation for each measure

synthetic slide projector data set focused at 3m, acquired at 90 cm

Acquired Data Acquired Data Acquired Data Acquired Data "directional information sorted by spatial origin"


Michael Goesele

Xavier Granier

Wolfgang Heidrich

Hans-Peter Seidel

Accurate Light Source Acquisition and Rendering

Keywords: CR Categories: I.3.3 [COMPUTER GRAPHICS]: Picture/Image Generation-Digitizing and Scanning, I.3.3 [COMPUTER GRAPHICS]: Three-Dimensional Graphics and Realism-Raytracing, I.4.1 [IMAGE PROCESSING AND COMPUTER VISION]: Digitization and Image Capture-Radiometry, I.4.1 [IM-AGE PROCESSING AND COMPUTER VISION]: Digitization and Image Capture-Sampling Modeling -Light Sources, Modeling -Object Scanning/Acquisition, Modeling -Physically-based Modelling, Rendering -Global Illumination, Rendering -Image-based Rendering, Rendering -Raytracing

Figure 1: Stages of light source measurement and rendering (from left to right): a) photo of flashlight -b) flashlight in measurement setupc) 2D reconstruction of the measured data in a virtual plane -d) measured light used in a global illumination simulation.

Overview and Related Work

Convincing and realistic light source models are an important prerequisite for photorealistic image synthesis. Given the exact de-sign of a luminaire (including geometry, materials, and emission characteristics of the illuminant), one can apply standard global illumination computations to incorporate complex luminaires into synthetic scenes. This approach, however, is impractical in most real-world scenarios since the precise design data represents an intellectual property of significant value, and is unlikely to be shared for rendering purposes. But even if such design data is made available, it does not account for manufacturing tolerances. However, relatively minor variations of parameters like the glass thickness in halogen light bulbs can have a significant impact on the illumination patterns generated. In addition, rendering based on the geometric model can be computationally very demanding.

As an alternative, one can incorporate measured light data from complex luminaires directly into the rendering system. Such data is provided by many companies in the form of goniometric diagrams [START_REF] Verbeck | A comprehensive light source description for computer graphics[END_REF], which represent measurements of the far field of a luminaire (i.e. directional information of the emitted light from a point light source). Unfortunately, the far field is only a faithful approximation of the emitted light when the light source is sufficiently far away from the object to be illuminated -as a rule of thumb at least five times the maximum luminaire dimension [START_REF] Ashdown | Near-Field Photometry: Measuring and Modeling Complex 3-D Light Sources[END_REF]]. Even improved models using spatial intensity maps or complex emitting geometry [START_REF] Rykowski | Source Modeling for Illumination Design[END_REF]] are still too weak to represent the near field of a complex light source.

On the other hand, a light field [START_REF] Gortler | The Lumigraph[END_REF]; [START_REF] Levoy | Light field rendering[END_REF] completely represents near and far field illumination of a light source, and can thus be used to represent the light source without knowing its geometry and internal structure.

Near Field Photometry

In his work on near field photometry [START_REF] Ashdown | Near-Field Photometry: A New Approach[END_REF][START_REF] Ashdown | Near-Field Photometry: Measuring and Modeling Complex 3-D Light Sources[END_REF], Ashdown has presented methods for measuring the light field of luminaires in both its near field and its far field components. A digital camera mounted on a robot arm observes a light source from many positions on a surrounding (hemi-)sphere. Both [START_REF] Rykowski | Source Modeling for Illumination Design[END_REF], and [START_REF] Jenkins | Source Imaging Goniometer Method of Light Source Characterization for Accurate Projection System Design[END_REF] employ a similar setup to acquire the light field of a luminaire, while [START_REF] Siegel | A General Near-Zone Light Source Model and its Application to Computer Automated Reflector Design[END_REF] replaced the camera lens with a pinhole.

On the rendering side, [START_REF] Heidrich | Canned lightsources[END_REF] have described relatively efficient algorithms based on a similar representation they called a canned light source. While their algorithms are faster than a brute-force rendering approach, they are still quadratic in the light field resolution. However, the similarity between the representations used by Ashdown and by Heidrich et al. illustrates that measurement and rendering can work on the same data structures, thereby allowing for an efficient processing pipeline from acquisition to image synthesis.

In near field photometry, a number of cameras are pointed at the luminaire to be measured (or, more practically, a single camera is moved around), and then the irradiance incident to the film plane is recorded using an imaging sensor (e.g. a CCD chip). The camera positions correspond to a sampling of some virtual sampling surface ❙ . In practice, the light source may produce arbitrarily high spatial frequencies on ❙ , in which case we have just introduced aliasing by not applying a low-pass filter before the sampling step.

This issue has been raised by [START_REF] Halle | Holographic stereograms as discrete imaging systems[END_REF] and was noted by [START_REF] Levoy | Light field rendering[END_REF]. They both propose to use the finite aperture of the camera lens as a low-pass filter by choosing the aperture size equal to the size of a sample on the camera plane. The shape of the filter kernel is thus a side effect of the lens system rather than a deliberate choice based on sampling theory.

Contributions

In this paper, we describe a measurement system that significantly improves on near field photometry by projecting the light field emitted by the light source into a finite basis before sampling. This is done using a simple optical system. The shape and support of the basis functions are specifically designed for a particular sampling scheme. Based on the resulting measurements we can exactly reconstruct the least-squares approximation of the true light field in our basis. Alternatively, we can reconstruct with a more efficient, shift-invariant filter, and obtain a close approximation to the leastsquares solution that is suitable for hardware accelerated rendering with the approach proposed by [START_REF] Heidrich | Canned lightsources[END_REF].

In comparison to previous acquisition approaches, our optical setup lets us move in very close to the surface of the light source. In particular, the distance between the light and the aperture of our system is smaller than the distance between the aperture and the imaging plane. This allows us to capture a much wider field of view in a single photograph.

On the rendering side, we introduce a solution to perform importance sampling of our light field representation. This allows us to integrate the acquired light sources into global illumination methods like the Photon Map Algorithm [Jensen 2001b]. In contrast to the work by Heidrich et al., this algorithm has constant time complexity, i.e. it is independent of the light field resolution.

The remainder of this paper is structured as follows: first, we present out conceptual approach to filtered light source measurement (Section 2) before we discuss the theoretical underpinnings of our method in Section 3. We then describe the physical setup we use for measurement in Section 4, and our rendering algorithm in Section 5. Finally, we present results based on both measurements and simulation (Section 6), and conclude with some ideas for future work (Section 7).

Basic Approach

The conceptual setup for our approach to filtered light source measurement is depicted in Figure 2. Light rays are emitted from the light source and hit a filter in a plane ❙ . This plane is opaque except for the finite support area of the filter which can be moved to a number of discrete positions on a regular grid. The filter is a semitransparent film, similar to a slide, containing the 2D image of a basis function ✟ ✵ ✐ ✭✉✁ ✈✮. The light falling through this filter hits a second plane ▼, on which we are able to measure the irradiance with a high spatial resolution.

When comparing this setup with the work on near field photometry [START_REF] Ashdown | Near-Field Photometry: A New Approach[END_REF][START_REF] Siegel | A General Near-Zone Light Source Model and its Application to Computer Automated Reflector Design[END_REF], ❙ roughly corresponds to the plane on which the pinholes or finite apertures of the cameras are located, while ▼ corresponds to the imaging sensor of the camera. The fundamental difference in our approach is that we can select an arbitrary filter ✟ ✵ ✐ ✭✉✁ ✈✮ which projects the light field emitted by the luminaire into an arbitrarily chosen function space. This results in an optical low-pass filtering of the light field before sampling, and thus helps avoiding aliasing artifacts.

In general, we will select a standard reconstruction filter ✟✐ ✭✉✁ ✈✮ for rendering, and its biorthogonal basis ✟ ✵ ✐ ✭✉✁ ✈✮ for measurement. This will make it necessary to perform measurements with filters that contain negative values, and can therefore not be implemented optically as an opacity. In this case, we split the filters into a positive and a negative component:

✟ ✵ ✐ ✂✄☎ ✭✉✁ ✈✮ ✿✆ ✝ ✟ ✵ ✐ ✭✉✁ ✈✮ ✟ ✵ ✐ ✭✉✁ ✈✮ ✕ ✞ ✞ ❡✠✡❡ ✟ ✵ ✐ ✂☛☎ ✭✉✁ ✈✮ ✿✆ ✝ ✞ ✟ ✵ ✐ ✭✉✁ ✈✮ ✕ ✞ ☞✟ ✵ ✐ ✭✉✁ ✈✮ ❡✠✡❡
We then perform separate measurements for both parts, and subtract the results in software. Similarly, we can rescale basis functions with values ❃ ✶, and correct for it in software.

In principle, one also has to perform low-pass filtering on the measurement plane ▼. This, however, is optically very difficult or even impossible to implement if the filter kernels overlap. Fortunately, it is technically easy to achieve a very high measurement resolution on ▼, so that we can consider the measured data as point samples, and implement an arbitrary filter kernel while downsampling to the final resolution.

Note that the filter kernels for ❙ and ▼ can be chosen independently as required by the reconstruction process. In this paper, we are using a piecewise quadratic basis to construct ✟ ✵ ✐ ✭✉✁ ✈✮ (see Section 4.2), and a box filter for ✟ ✵ ❦• ✭s✁ t✮. The latter works best with our importance sampling algorithm (Section 5 and Appendix B).

Theoretical Framework

Before we discuss the theory behind the proposed method in detail, we first introduce the mathematical notation used throughout this document. This notation is summarized in Table 1 andFigure For the measurement, we assume that the light field emitted by the light source is well represented by a projection into a basis ❢✠✐ ❦•✭✉✁ ✈✁ s✁ t✮❣ ✐ ❦•✷) ) :

▲✭✉✁ ✈✁ s✁ t✮ ✙ ⑦ ▲✭✉✁ ✈✁ s✁ t✮ ✿☎ ❳ ✐✆ ✆❦✆• ✠✐ ❦•✭✉✁ ✈✁ s✁ t✮ ✄ ▲✐ ❦•✝ (1)
We assume that ✠✐ ❦• has local support, and ✞✁ •✁ ✡✁ and ☛ roughly correspond to translations in ✉✁ ✈✁ s✁ and t, respectively. Note, however, that the translated basis functions will not in all cases have the same shape, i.e. ✠ ✐ ☞ ☞ ❦ ☞ • ☞ ✭✉✁ ✈✁ s✁ t✮ may not be an exact copy of ✠✐ ❦•✭✉✁ ✈✁ s✁ t✮ in general.

We also define two additional sets of basis functions, one for measuring and one for reconstruction. For reconstruction we use a 1D basis ❢✟✐❣ ✐✷) ) with the property ✟✐ ✭①✮ ☎ ✟✭① ✰ ✞✮. The 4D reconstruction basis is then given as the tensor product basis

✟ ✐ ❦•✭✉✁ ✈✁ s✁ t✮ ✿☎✟✐ ✭✉✁ ✈✮ ✄ ✟ ❦•✭s✁ t✮ (2) ☎✟✐ ✭✉✮ ✄ ✟ ✭✈✮ ✄ ✟ ❦ ✭s✮ ✄ ✟ •✭t✮✝ (3)
For measurement, we use the biorthogonal (or dual) ❢✟ ✵ ✐ ✭①✮❣ ✐✷) )

of the reconstruction basis with

✌ ✶ ✍✶ ✟ ✵ ✐ ✭①✮ ✄ ✟ ✐ ☞ ✭①✮ ❞① ☎ ✚ ✎ ❀ ✏✑ ✞ ☎ ✞ ✵ ✒ ❀ ❡✓✔❡ (4)
and again we use a tensor-product construction for the 4D basis.

Measured Irradiance

Our approach is based on measuring the irradiance ❊♠✂✭s✁ t✮ in the ✭s✁ t✮-plane that is caused by the incident radiance

▲♠✂✭✉✁ ✈✁ s✁ t✮ ☎ ✟ ✵ ♠✂ ✭✉✁ ✈✮ ✄ ▲✭✉✁ ✈✁ s✁ t✮.
In Section 4.1 we will discuss two physical setups for performing this kind of measurement. The result of such a measurement is

❊♠✂✭s✁ t✮ ☎ ✌ ✶ ✍✶ ✌ ✶ ✍✶ ❝✕✔ ✖ ✗ | ✖ ✄ ✟ ✵ ♠✂ ✭✉✁ ✈✮ ✄ ▲✭✉✁ ✈✁ s✁ t✮ ❞✉ ❞✈ ✙ ✌ ✶ ✍✶ ✌ ✶ ✍✶ ❝✕✔ ✖ ✗ | ✖ ✄ ✟ ✵ ♠✂ ✭✉✁ ✈✮ ✄ ⑦ ▲✭✉✁ ✈✁ s✁ t✮ ❞✉ ❞✈ ☎ ❳ ✐✆ ✆❦✆• ✌ ✶ ✍✶ ✌ ✶ ✍✶ ❝✕✔ ✖ ✗ | ✖ ✄ ✟ ✵ ♠✂ ✭✉✁ ✈✮ (5) ✄ ✠✐ ❦•✭✉✁ ✈✁ s✁ t✮ ✄ ▲✐ ❦• ❞✉ ❞✈✝
The geometric term ❝✕✔ ✖ ✗❂| ✖ is composed of the distance | of the point on the ✭✉✁ ✈✮-plane from the point on the ✭s✁ t✮ plane, as well as the cosine of the angle ✗ between the plane normals and the vector connecting the two points. Note that this term also accounts for any differences in the parameterizations on the two planes (i.e. different grid spacings).

Exact Reconstruction

We now describe an exact reconstruction algorithm given the measurements ❊♠✂ . To this end, we first define what the relationship between the basis functions ✠✐ ❦• and the reconstruction and measurement bases should be. We define

✠✐ ❦•✭✉✁ ✈✁ s✁ t✮ ✿☎ | ✖ ❝✕✔ ✖ ✗ ✄ ✟✐ ✭✉✁ ✈✮ ✄ ✟ ❦•✭s✁ t✮✝ (6)
Inserting this definition into Equation 5, and using the biorthogonality relationship (Equation 4) yields

❊♠✂✭s✁ t✮ ☎ ❳ ❦ ❳ • ✟❦•✭s✁ t✮ ✄ ▲♠✂❦•✝ (7)
To determine which reconstruction filter to use, we now rewrite Equation 1 using Equations 6 and 7:

⑦ ▲✭✉✁ ✈✁ s✁ t✮ ☎ ❳ ♠✆✂✆❦✆• ✠♠✂❦•✭✉✁ ✈✁ s✁ t✮ ✄ ▲♠✂❦• ☎ ❳ ♠✆✂✆❦✆• | ✖ ❝✕✔ ✖ ✗ ✄ ✟♠✂✭✉✁ ✈✮ ✄ ✟ ❦•✭s✁ t✮ ✄ ▲♠✂❦• (8) ☎ ❳ ♠✆✂ | ✖ ❝✕✔ ✖ ✗ ✄ ✟♠✂✭✉✁ ✈✮ ✄ ❊♠✂✭s✁ t✮
This indicates that we can exactly reconstruct ⑦ ▲, the projection of ▲ into the basis ❢✠✐ ❦•❣ by using the reconstruction filter

| ✖ ❂ ❝✕✔ ✖ ✗ ✄ ✟♠✂✭✉✁ ✈✮.
Note that this reconstruction filter does contain a shift-variant component in form of the geometric term. Interestingly, this geometric term will cancel out for the Photon Map particle emission algorithm we discuss in Section 5. However, for other rendering algorithms such as the ones proposed by Heidrich et al. [START_REF] Heidrich | Canned lightsources[END_REF]], the reconstruction step should ideally only involve a convolution with a shift-invariant filter kernel. In Appendix A we discuss an approximate reconstruction scheme that satisfies this property.

Physical Setup and Measurements

In this section we explain how the theoretical framework can be transformed into a practical acquisition method. We propose two different hardware setups and discuss their advantages and disadvantages. Finally, we describe the filter kernels used in our measurements.

Possible Implementations

The first setup we propose -Setup A -is based on a standard photographic camera where the lens is removed and replaced by the filter ✟ ✵ ✐ (see left side of Figure 3). As with conventional light fields, the camera has to be positioned at the vertices of a regular grid within the ❙ plane and images can be captured. Although this setup is compact and easy to handle, it does have some drawbacks. The fixed and often quite small size of the imaging sensor (e.g. a CCD chip) imposes a strong limitation on the size of the filter kernel and the incident angles of illumination covered by the measurement. In addition, the distance between the filter plane and the imaging plane determines the focal length of the camera (similar to a pinhole camera). A large distance leads to a tele-lens effect and further limits the useful incident angles. The filter ✟ ✵ ✐ ✭✉✁ ✈✮ projects an image on a diffuse reflector which is captured by a camera. These problems can be avoided with Setup B. Here, the filter

✟ ✵
✐ is used to project an image of the light source onto a diffusely reflecting screen. If a high quality diffuse (Lambertian) reflector material is available, the filter size and viewing angles are (almost) only limited by the size of the reflector. On the other hand, the larger setup makes it potentially harder to use, and errors due to indirect light reflecting from parts of the equipment such as the reflection screen or the environment tend to increase. Furthermore, the calibration effort increases since the intrinsic and extrinsic parameters of the camera system relative to the projected image need to be recovered.

Filter Design

The theoretical framework and the implementations detailed above allow for the use of arbitrary filters ✟ ✐ and ✟ ✵ ✐ with finite support. It is therefore possible to design filters adapted to the specific properties of the light source at hand, for example taking the structure of the light source or the intended use of the data into account. For our measurements we use a piecewise quadratic basis function

✟ ✐✄bellq ✟ ✐✄bellq ❂ ✽ ❃ ❁ ❃ ✿ ✶ ☎ ✷ ① ✆ •①• ✔ ✝ ✆ ✷ ✭•①• ☎ ✶✮ ✆ ✝ ✆ ✞ •①• ✔ ✶ ✠ else ✡ (9)
Although it would also be possible to use the classical quadrilinear basis, we choose this one due to its ❈ ✝ continuity that helps preventing Mach bands. Figure 4 shows a plot of this function and its dual.

Note that some parts of the dual basis are negative. Using this function as a filter will result in negative coefficients. This will occur in places with strong intensity gradients in the original light field. There are two ways of dealing with this situation. The first one is to clamp the value of the coefficients after measurement, guaranteeing that each query in the light field will be positive. However, this approach will remove fine details in the light field and contradicts our filtering approach. We therefore choose to keep the negative coefficients in the data structures, and to only clamp negative results after reconstruction. With this approach, the positive parts of the reconstructed light field remain unaltered.

Rendering

In this section, we discuss the integration of our measured light sources into a particle based global illumination system. We chose the Photon Map algorithm [START_REF] Jensen | Global Illumination Using Photon Maps[END_REF]Jensen 2001a;Jensen 2001b] for rendering, as it is one of the most common and simple to use methods available. We present a solution to emit particles from our measured data, and we also describe some adjustments to the classical Photon Map algorithm. Note that our importance sampling approach does not only apply to Photon Maps, but to any algorithm that requires shooting from the light source such as particle tracing [START_REF] Walter | Global illumination using local linear density estimation[END_REF]] or bidirectional path-tracing [START_REF] Lafortune | Bi-directional Path Tracing[END_REF].

Importance Sampling

For the scene description, we attach the light field to a virtual geometric object, representing an impostor for the area from which light is emitted. This can be, for example, the glass surface of a car head light, or simply a planar polygon representing one of the light field planes. This geometry helps in positioning the light source in the 3D scene and can also be used to aid the sampling in our rendering algorithm.

To emit a particle from our light source representation, we could then use the standard method of selecting a random point on the impostor geometry and then a random direction. While this approach might be sufficient for a relatively uniform light source, it can be inefficient when high spatial and directional frequencies are present.

✒ ✓ Φ i ( u) Φ i+1 ( u) Φ j ( v) Φ j+1 ( v) Φ l ( t) Φ k ( s)
Figure 5: A light field cell for PDF computation formed by a single pixel on the measurement plane ▼ and by ✶❂✹ of a filter on the filter plane ❙ . This corresponds to the region on ❙ where ✟✐✭✉✮, ✟✐✰ ✭✉✮, ✟• ✭✈✮, and ✟•✰ ✭✈✮ overlap.

We therefore develop an importance sampling strategy for our light field representation.

The particles that we want to emit have to estimate the total energy ❊ coming from the light field:

❊ ✁ ) ✂✄☎✄❀✄✆ ⑦ ▲✭✉✝ ✈✝ s✝ t✮ ❝✞s ✷ ✒ | ✷ (10)
In order to achieve a constant energy for each photon, the probability density function (PDF) of the photon emission has to be:

♣✭✉✝ ✈✝ s✝ t✮ ✁ ✶ ❊ ⑦ ▲✭✉✝ ✈✝ s✝ t✮ ❝✞s ✷ ✒ | ✷ (11)
To generate random photons distributed according to this PDF, we use the standard approach of generating uniformly distributed random points in 4D, and passing them through the inverted cumulative density function (CDF, [START_REF] Pitman | Probability[END_REF]). If we have preinverted and tabulated the CDF, then generating a new random photon according to this importance is possible in constant time, i.e. independent of the light field resolution.

One complication in our case is, however, that we have the PDF represented in terms of some complicated and overlapping basis functions instead of the usual point sampled (i.e. tabulated) representation. We therefore do not directly generate and invert the CDF from the basis representation, but rather choose a two-steps approach. We first compute a 4D table of constant probabilities for selecting a 4D cell from the light field (see Figure 5). A single 4D cell corresponds to all the directions coming from one pixel on the measurment plane to ✶❂✹ of a filter on the filter plane (corresponding to the minimum overlap between two neighboring filters in our basis). The probability on this cell can be computed by some simple summations over the light field coefficients. The details are described in Appendix B. Then, as we now have a constant probability for each cell, we can invert this table more easily than the continuous representation of the PDF. Note that all these steps are performed only once in a preprocessing stage. The resulting table is part of the light source description that we save on disk.

Direct Illumination

As pointed out in the previous section, the light field is attached to a virtual base geometry. The irradiance at a point ① can easily be computed by sampling this geometry. But as in the case of indirect illumination, this approach can be quite inefficient depending on the directional variation of the light source.

Our approach is therefore to separate the Global Map of the traditional Photon Map algorithm into two: the Direct Map (i.e. the direct impacts from the lights sources), and a new Global Map (i.e. impacts after one diffuse reflection). With this separation, we have two possibilities for computing the direct illumination: we can either choose to reconstruct it straight from the Direct Map, or we can use the Direct Map for generating an importance to guide a standard Monte Carlo lighting integration.

The former method promises to be quite efficient for highly directional light sources. The quality of this reconstruction can also be increased by storing a larger photon number in the Direct Map, while keeping a low number for the Global Map. This approach involves either limiting the integration to regions where some direct photons are present in the neighborhood, or building an importance sampling table from the nearest particles in the Direct Map (similar to the approach used for the Global Map by [START_REF] Jensen | Importance Driven Path Tracing Using the Photon Map[END_REF]).

The idea of identifying the direct impacts was already introduced in [START_REF] Keller | Efficient importance sampling techniques for the photon map[END_REF]. However, instead of using two different maps, they use the classical Global Map and extend the photon data structure by storing a light source identifier. While this approach is well suited for determining which light source the particle is coming from, it is less efficient for finding all direct impacts, which is what we require here.

In general, the separation of the Global Map into two maps gives us a lot more flexibility in the reconstruction process. As pointed out above, the user can choose between different methods for both the direct and the indirect illumination. Note that this algorithm is independent of our light source representation and can also be used in other Photon Map implementations. This is also true for the following section.

Other Photon Map Improvements

Our implementation uses a slightly different algorithm to find the ◆✠ closest photons in one of the two maps. We find our approach to be better in determining which photons actually contribute to the irradiance estimate.

In the following, ✡ ❞✐ will denote the impact direction of a photon ♣✐, and ①☛ will denote its position. We then search for the particles closest to a position ① with a normal ✡ ♥. During the search, we only consider the particles that have a negative dot product ✭ ✡ ❞✐ ☞ ✡ ♥✮, meaning that the particle arrives from above the surface. This guarantees that only the incident photons on one side of a wall, for example, will contribute to the reconstructed illumination.

Then, with an approach similar to the one described by [START_REF] Lastra | A density estimation technique for radiosity[END_REF], the density is estimated on the tangent plane. The particles ♣✐ are projected into the tangent plane at ①, and the Euclidean distance of that projection from ① is used as the distance metric. The distance is thus defined as:

✍✭①✝ ♣✐✮ ✷ ✁ ✌ ✷ ✎ ✏✌✭① ✎ ①✐ ☞ ✡ ❞✐✮ ✑ ❦① ✎ ①✐❦ ✷ ✌ ✁ ✭① ✎ ①✐ ☞ ✡ ♥✮ ✭ ✡ ❞✐ ☞ ✡ ♥✮
This allows us to increase the accuracy of the reconstruction, especially for non-planar surfaces and corners.

Results

We first present results from a simulation of the complete process, from acquisition (with Setup A) to rendering as a validation of our approach. We then discuss our measurement setups and show renderings of acquired data sets.

Simulations

Before building the real setup, we validate our approach by simulating the acquisition of a strongly directional light field such as a slide projector. It has an aperture of 50 mm and a focal length of 59 mm. We use an inverted SIGGRAPH 2003 logo as a slide, to create strong directional variations (see Figure 6). Simulations are done using Setup A with a filter size of 5.068 mm✂5.068 mm and the quadratic basis. The simulated camera has an imaging sensor with a size of 27.15 mm✂18.1 mm and a resolution of ✶ ✁ ✂ ✶✁✁ pixels. The distance between the filter and the imaging sensor is ✸✄ mm, and the distance between the slide projector and the cam- era is ✷✁ mm. We simulate the acquisition of a ✹✶ ✂ ✹✶ array of images. We store the data in a signed version of the RGBE format [START_REF] Ward | Graphics Gems II[END_REF]] to have a more efficient usage of the memory. The resulting simulation uses about 96 MB of memory.

Original configuration

Moving plane

In-focus plane Based on this data, we then computed a sequence of projections on a plane moving from the lens of the slide projector to the projector's focus plane, to show the directional changes (see Figure 6). The in-focus image shows that we manage to capture the whole shape of the light source, even if the small details are not present which is expected for this setup geometry and the optical filtering. We also used this data to compute a global illumination solution that validates the particle emission approach (see Figure 7). This scene has about 15,000 primitives, and it took about 2 minutes to compute the inverted probability table (size ✸✷ ✂ ✸✷ ✂ ✶✷☎ ✂ ✶✷☎), 4 seconds for emitting 40,000 photons, and 52 minutes to render one image. 1

From light field

Measurements

In order to measure real light sources we built systems for each of the two setups proposed in Section 4.1. The filters were printed with 1016 dpi resolution on high density film using a calibrated LightJet 2080 Digital Film Recorder. This produces slides with a dynamic range of more than 1 : 50,000, which are then mounted into suitable frames. The light sources to be measured -we used a bike light as well as several different flashlight configurations -are mounted on manually operated translation stages so that they can be shifted in front of the camera with high precision and repeatability (see Figure 8). The images were captured with a Kodak DCS 560 1 All timings are for a 1GHz Pentium III Workstation. digital camera in our measurement lab which is covered with dark material to avoid interference with the measurements [START_REF] Goesele | Building a Photo Studio for Measurement Purposes[END_REF].

For Setup A, a custom-built filter holder was attached to the camera instead of a lens, so that the filters can easily be exchanged without affecting the geometry of the setup (see Figure 9, left). Highdynamic range acquisition techniques [START_REF] Debevec | Recovering High Dynamic Range Radiance Maps from Photographs[END_REF][START_REF] Robertson | Dynamic Range Improvement Through Multiple Exposures[END_REF] were used to acquire the images.

For Setup B, we used a 1.2 m ✂ m flat board that was spray- painted with ultra matte paint as a projection screen. A black board with an integrated filter holder was used to mount the filter and shield the light source from the projection screen. To avoid occlusion problems we mounted the camera above the whole setup and equipped it with a 24 mm shift lens allowing us to capture the whole projection screen. The camera parameters were recovered using standard camera calibration techniques and the images were rectified and registered to the filter setup. Figure 9 (right) shows an overview of Setup B. A 5✂5 median filter was used to suppress dark current noise resulting from longer exposure times in this setup; remaining noise appears as background illumination in the light field.

Measured Data Sets

We acquired several data sets with these setups. Table 2 summarizes the acquisition geometry, Table 3 gives an overview over the number of views captured and the filters used. All acquired images were downsampled to achieve a resolution of 300 ✂ 300 pixels on the measurement plane ▼. Acquisition times are on the order of Table 3: Parameters for the acquired data sets. Note that the width of the reconstruction filter is ✶❂✷ of the width of the dual filter, the filter spacing on the filter plane ❙ is ✶❂✹ of the width of the dual filter. Only the Mag-Lite data set was captured with Setup A. Due to geometry constraints imposed by the dimensions of the camera (see Table 2), we were only able to capture a small part of the light field.

Figure 10 shows on the left the light source pointed at a projection screen and on the right onto the crypt. The captured part corresponds roughly to the central dark region in the projection of the real light source (see top of Figure 11). The low resolution is due to the comparably large filter size. Using a smaller filter would increase both the resolution and the number of views required. In this configuration, Setup A serves as a proof-of-concept but can probably not be used to acquire a complete model of a real light source. In Figure 11 we show the light pattern that a Mini Mag-Lite flash light causes on a plane at 75 cm distance for a specific focus setting on that flashlight. The top image is a photograph of the flashlight pointed at a diffuse surface. The bottom row shows two synthetic images that are reconstructed from measurements. The left image was generated from a dataset where a dual filter width of 30 mm has been used, while a dual filter of 20 mm width was used for the right image. As expected, the larger filter size results in a slightly more blurry image on the left. In comparing the reconstructed images with the photograph, we can see that there are some tone-mapping issues that result in slightly different colors. However, the main features of the light cone, such as the black ring in the center, are represented very well in both reconstructions.

In Figure 12 very well in the renderings. Furthermore, due to the effective size of the luminaire (5 cm✂2.5 cm) and its strongly varying emission characteristics, the near field is quite extended and inhomogeneous. In the rendered image of the Bike Light model projected at different distances, variations of the near field are clearly visible. Figure 14 shows a global illumination solution based on this dataset. Note both the distinctive patterns on the wall and the stripe patterns on floor and ceiling.

Finally, we have performed the same test for a flashlight with 3 LEDs instead of a light bulb (see left of Figure 13 -the fourth, smaller LED is a low battery light, and is not used in normal operation). Although the actual flashlight is identical to the one from Figure 11 except for the light bulb and the reflector, the cone of light produced by this setup is quite different from the one before. In particular, the illumination is much more uniform, although roughly triangular due to the LED layout. This illustrates the importance of using actual measurements from real world light sources in image synthesis, since even small changes in the reflectors can cause major differences in appearance.

Conclusions

In this paper we presented a novel approach for accurate light acquisition. A model is measured using a simple optical system to project the real light source into a finite element basis. This projection allows a correct optical pre-filtering of the acquired light field. We introduced a theoretical framework for this approach and two possible setups for real acquisitions. The tests we performed to validate this approach both on a simulated process and on the two physical setups, have shown that our setup can faithfully represent rather complex near and far field information with a relatively small number of measurements. The fact that we perform optical lowpass filtering avoids aliasing artifacts such as missing directional peaks. In our evaluations we found that the second setup, although larger and harder to calibrate, is more practical for real-world light sources, since it covers a larger field-of-view.

All our acquisitions where done manually. As only linear movements are involved, the entire measurement can be speed up easily by using a motor driven translation stage. This will lead to an entirely automatic process, that would allow us to further increase the resolution of the measurements.

We also presented an efficient algorithm to integrate our new light representation into particle-based global illumination systems. Using a pre-inverted cumulative density function, this method allows a constant time importance sampling to emit those particles. Combined with the acquisition step, we therefore a complete pipeline to bring real light into virtual worlds.

A An Approximate Shift-Invariant Reconstruction

For an approximate shift-invariant reconstruction, we start with the measured irradiance from Section 3.1, but we define the relationship between the ✠✐ ❦• and the measurement and reconstruction bases differently than we did for the exact reconstruction (Section 3.2). Here, we define ✠✐ ❦•✭✉✁ ✈✁ s✁ t✮ ✿✂ ✟✐ ✭✉✁ ✈✮ ✄ ✟❦•✭s✁ t✮.

From this, we get

❊♠☎✭s✁ t✮ ✂ ❳ ✐✆ ✆❦✆• ) ✶ ✝✶ ) ✶ ✝✶ ❝✞✡ ✷ ✒ | ✷ ✄ ✟ ✵ ♠☎ ✭✉✁ ✈✮ (12) ✄ ✟✐ ✭✉✁ ✈✮ ✄ ✟❦•✭s✁ t✮ ✄ ▲✐ ❦• ❞✉ ❞✈
This is still a shift variant filter, since the geometric term depends on the integration variables ✉ and ✈, and can therefore not be moved outside the integral. However, if the distance ❞ between the ✭✉✁ ✈✮-plane and the ✭s✁ t✮-plane is large compared to the support of ✟✐ ✭✉✁ ✈✮, and if ✒ is small, then this term is well approximated by one constant per point on the ✭s✁ t✮-plane:

❝✞✡ ✷ ✒ ☛ ✉✁ ✈✁ s✁ t✮ |✭✉✁ ✈✁ s✁ t✮ ✷ ✙ ❝✞✡ ✷ ✒❦•✭s✁ t✮ |❦•✭s✁ t✮ ✷ ☞ (13)
This yields an approximation of the measured irradiance:

❊♠☎✭s✁ t✮ ✙ ❳ ❦✆• ❝✞✡ ✷ ✒❦•✭s✁ t✮ | ❦•✭s✁ t✮ ✷ ✄ ✟❦•✭s✁ t✮ ✄ ▲♠☎❦•☞ (14)
Since both the approximate geometric term and ✟❦•✭s✁ t✮ are known, it is in principle possible to compute

❊ ✵ ♠☎ ✭s✁ t✮ ✙ ❳ ❦✆• ✟❦• ✭s✁ t✮ ✄ ▲♠☎❦• (15)
by de-convolution. In practice, this is only feasible for basis functions ✟❦•✭s✁ t✮ with a small support. We do not expect this to be a major problem, however, since the practical measurement setups discussed in Section 4.1 have a very high resolution on the ✭s✁ t✮-plane, so that a bilinear or even box filter can be used. Again, we apply the definition of ✠✐ ❦• to determine the appropriate reconstruction filter:

⑦ ▲✭✉✁ ✈✁ s✁ t✮ ✂ ❳ ♠✆☎✆❦✆• ✠♠☎❦•✭✉✁ ✈✁ s✁ t✮ ✄ ▲ ♠☎❦• ✂ ❳ ♠✆☎✆❦✆• ✟♠☎✭✉✁ ✈✮ ✄ ✟❦•✭s✁ t✮ ✄ ▲♠☎❦• (16) ✙ ❳ ♠ ❳ ☎ ✟♠☎✭✉✁ ✈✮ ✄ ❊ ✵ ♠☎ ✭s✁ t✮☞
Thus, an approximate reconstruction of ⑦ ▲, the projection of ▲ into the basis ❢✠✐ ❦•❣ is obtained using the shift-invariant reconstruction filter ✟♠☎✭✉✁ ✈✮. For example, if ✟ is a hat function, the reconstruction process is reduced to quadri-linear interpolation, which is frequently used for light field rendering algorithms.

The quality of this approximation depends on the error introduced by assuming the geometric term constant over the support of basis function ✟✐ ✭✉✁ ✈✮ in Equation 13. An analysis of this error [START_REF] Heidrich | Image-based measurement of light sources with correct filtering[END_REF] shows that it is below 10% if the filter is no larger than 10% of the distance between ▼ and ❙ . That error is reduced to about 1% for filter sizes no larger than 2% of that distance.

B Importance Sampling Details

This appendix details the computation of the 4D probability table introduced in Section 5.1.

To simplify the following equations and the creation of the 4D table, we make the following assumptions: (a) -the configuration is centered; (b) -the overlap between two consecutive basis functions on the sampling plane ❙ correspond to half of its size for width and height; (c) -the width and height shift (s ✇ and s❤, in pixel number) between two measurements are integers; (d) -the basis function ✟✐ is symmetric; (e) -the basis function on the measurement plane ▼ is the box filter; (f) -the support of the original reconstruction basis function ✟ is ❬✂✶✁ ✶✄ and its integral is ✶ on this support. With this assumption, the integrals of the reconstruction basis function are (if their dimension are ☎✆ ✝ ✞✆ ):

) ✟✐ ✭✉✮ ❂ ☎✆ ✷ ✁ ) ✟ ✭✈✮ ❂ ✞✆ ✷
With this we can compute the total energy emitted by the light field (see Equation 10):

❊ ❂ ☎✆ ✝ ✞✆ ✝ ☎ ✠ ✝ ✞ ✠ ✹ ❳ ✐✡ ✡❀✡☛ ▲ ✐ ❀☛
Here, ☎ ✠ and ✞ ✠ represent the dimensions of a pixel on the measurement plane. The 1D-intervals defining a cell can be described by the following equations (■ ✆ ❲✡☞ and ■ ✆ ❍✡✌ on the sampling plane, and ■❲ ✡☞✡✍ and ■❍ ✡✌✡✎ on the measurement plane -these depend on the selected intervals on the sampling plane): The first probability concerns the choice of one cell, while the others represent the conditional probabilities of selecting one interval only.

Figure 2 :

 2 Figure 2: Cross section through the conceptual acquisition setup.

Figure 3 :

 3 Figure 3: Possible setups for the acquisition. Setup A (left): The camera lens is replaced by a filter and the light field is directly projected onto the imaging sensor inside the camera. Setup B (right): The filter ✟ ✵✐ ✭✉✁ ✈✮ projects an image on a diffuse reflector which is captured by a camera.

Figure 4 :

 4 Figure 4: A quadratic reconstruction basis and one dual. Left: function ✟ ✐ (green) and a dual ✟ ✵ ✐ ✂ (red) used in this paper. Center and right: 2D filter images used for the acquisition (split into positive and negative components and normalized).

Figure 6 :

 6 Figure 6: Projections from the simulated light field of a slide projector. The first image (left-up) is a projection on a plane in focus, the second one (left-down) is a projection on a plane closer to the slide projector.

Figure 7 :

 7 Figure 7: Global illumination solution for the simulated acquisition data.

Figure 8 :

 8 Figure 8: Left: The Mini Mag-Lite and the bike light mounted to the translation stage. Right: The negative part of the filter illuminated by the bike light. Both images show parts of Setup B.

Figure 9 :

 9 Figure 9: Left: Setup A. The filter holder mounted to the camera in front of the lens of a light source. Right: Setup B. The projection screen is observed by the camera located above the setup. The black board contains the filter holder. The light sources can be moved on the table parallel to the filter and the projection screen.

Figure 10 :

 10 Figure 10: Mag-Lite acquisition with Setup A: Left, a projection on a distant (3 m) plane, Right: a global illumination solution (Importance table construction -20 s, Particle emission -40 s, Rendering 20 min). The acquired data set was converted to grayscale.

Figure 11 :

 11 Figure 11: Top: photograph of the Mag-Lite illuminating a plane at 75 cm distance. Bottom left: reconstruction of the same setup from data measured with a filter of 30 mm width. Bottom right: reconstruction from data measured with a filter of 20 mm width.With Setup B, we can generally capture a much larger portion of the light field of the light source due to the ✻✼ ✍ viewing angle of the configuration. Additionally, a much larger filter kernel can be used for the filter plane while still improving the resolution of the data set.In Figure11we show the light pattern that a Mini Mag-Lite flash light causes on a plane at 75 cm distance for a specific focus setting on that flashlight. The top image is a photograph of the flashlight pointed at a diffuse surface. The bottom row shows two synthetic images that are reconstructed from measurements. The left image was generated from a dataset where a dual filter width of 30 mm has been used, while a dual filter of 20 mm width was used for the right image. As expected, the larger filter size results in a slightly more blurry image on the left. In comparing the reconstructed images with the photograph, we can see that there are some tone-mapping issues that result in slightly different colors. However, the main features of the light cone, such as the black ring in the center, are represented very well in both reconstructions.In Figure12we show similar images for the light field emitted by a bike light (the light source itself is depicted in the top left of the figure). The top right of the figure again shows a photograph of the light source illuminating a plane at 75 cm distance. The image on the lower right is a reconstruction for the same distance, while the image on the lower left shows the reconstruction at a distance of only 18.7 cm.The narrow vertical stripes are an interesting feature of the bike light that should be preserved in a model. They are mainly emitted form the upper and lower edges of the luminaire and are reproduced

  Figure 12: Acquisition on the Bike Light model. The upper right image shows a plane illuminated by the real light source, at the distance of 75 cm. The lower right image shows the same plane illuminated by the acquired data, while the lower left image shows the result for a plane at a quarter of that distance.

Figure 13 :

 13 Figure 13: Photograph and reconstruction of the LED Flashlight.

Figure 14 :

 14 Figure 14: Exploring the cloister with the Bike Light. Importance table construction -37 s, Particle emission -12 s, Rendering 30 min.

  ✰ ✶✮✑✷ ✰ ❬♥✁ ♥ ✰ ✶✄✮ ✝ ☎✆ ✽♥ ✒ ❬✓✿✿✏✆ ✄ ■ ✆ ❍✡✌ ❂ ✭✂✭✔✆ ✰ ✶✮✑✷ ✰ ❬♦✁ ♦ ✰ ✶✄✮ ✝ ☎✆ ✽♦ ✒ ✠ ✰ s❤ ✂ ✶✄Note that ✏✆ ✝ ✔✆ represents the resolution on the sampling plane (width✝height), and ✏✠ ✝ ✔ ✠ on the measurement plane. A set of light directions ■☞✌✍ ✎ on the sampling plane with the dimension ☎✆ ✑✷✝✞✆ ✑✷ from a existing pixel in the dataset ■ ✠ ❲✡☞✡✍ ✝■ ✠ ❍✡✌✡✎ . On these intervals, we define the events ❯☞ ❂ ✭✉ ✒ ■ ✆ ❲✡☞ ✮, ❱✌ ❂ ✭✈ ✒ ■ ✆ ❍✡✌ ✮, ✕☞✍ ❂ ✭s ✒ ■ ✠ ❲✡☞✡✍ ✮ and ❚✌ ✎ ❂ ✭t ✒ ■ ✠ ❍✡✌✡✎ ✮. We can then compute the probabilities the conditional probabilities of these events: with ❦ ✐☞✍ ❂ ♣✂ ✭✤✂ ♥✰ ✶✮s ✇ and • ✌✎ ❂ q✂ ✭✥ ✂ ♦✰ ✶✮s❤.

Table 1 :

 1 2. Notation (overview) -see also Figure2.

	Symbol	Meaning
	✠✐ ❦•✭✉✁ ✈✁ s✁ t✮	basis function for approximating the light field
	✟✐	1D basis used for reconstruction
	✟✐ ❦•✭✉✁ ✈✁ s✁ t✮ ✟ ✵ ✐ ✟ ✵ ✐ ❦• ✭✉✁ ✈✁ s✁ t✮	4D tensor product basis for reconstruction biorthogonal 1D basis used for measurement 4D tensor product basis for measurement
	▼	surface on which the irradiance is measured (measurement surface, or ✭s✁ t✮-plane)
	❙	surface on which the optical filters are placed (sampling surface, or ✭✉✁ ✈✮-plane)
	▲✭✉✁ ✈✁ s✁ t✮	radiance passing through ✭✉✁ ✈✮ on ❙ and ✭s✁ t✮ on ▼
	⑦ ▲✭✉✁ ✈✁ s✁ t✮ ▲♠✂✭✉✁ ✈✁ s✁ t✮	projection of ▲✭✉✁ ✈✁ s✁ t✮ into basis ❢✠✐ ❦•✭✉✁ ✈✁ s✁ t✮❣ radiance ✟ ✵ ♠✂ ✭✉✁ ✈✮ ✄ ▲✭✉✁ ✈✁ s✁ t✮ projected through one filter ✟ ✵ ♠✂ ✭✉✁ ✈✮
	❊♠✂✭s✁ t✮	irradiance caused by ▲♠✂✭✉✁ ✈✁ s✁ t✮ on the measurement surface ▼

Acknowledgements

We would like to thank Ian Ashdown and the anonymous Siggraph reviewers for valuable comments on earlier versions of the paper. Thanks to REVES/INRIA Sophia-Antipolis for allowing us to use the GIS Global Illumination platform for development in this project, and to UBC's Institute of Applied Mathematics for the access to their PC-Cluster for the measurement simulation. Special thanks to the MPI technical staff for building most of the measurement apparatus and to Oliver Ashoff who helped to perform the measurements. This work was funded by the PIMS Post-doctoral Fellowship program, the BC Advanced Systems Institute, and by the DFG Schwerpunktprogramm V3D2.