IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

1

Recovering Inner Slices of Layered Translucent
Objects by Multi-frequency Illumination
Kenichiro Tanaka, Student Member, IEEE, Yasuhiro Mukaigawa, Member, IEEE,
Hiroyuki Kubo, Yasuyuki Matsushita, Senior Member, IEEE, and Yasushi Yagi, Member, IEEE
Abstract—This paper describes a method for recovering appearance of inner slices of translucent objects. The appearance of alayered translucent object is the summed appearance of all layers, where each layer is blurred by a depth-dependent point spreadfunction (PSF). By exploiting the difference of low-pass characteristics of depth-dependent PSFs, we develop a multi-frequency
illumination method for obtaining the appearance of individual inner slices. Specifically, by observing the target object with varying the
spatial frequency of checker-pattern illumination, our method recovers the appearance of inner slices via computation. We study theeffect of non-uniform transmission due to inhomogeneity of translucent objects and develop a method for recovering clear inner slices
based on the pixel-wise PSF estimates under the assumption of spatial smoothness of inner slice appearances. We quantitatively
evaluate the accuracy of the proposed method by simulations and qualitatively show faithful recovery using real-world scenes.
Index Terms—Descattering, layer separation, image restoration, projector-camera system.

F

1

I NTRODUCTIONS

EEING inside a translucent object is a difficult task
because its appearance is superposition of light rays
emitted from every inner layers that are blurred due tosubsurface scattering. Because observing the internal appearance of objects is of broad interest in medical, art,
and industrial inspections, various imaging techniques havebeen developed for achieving this goal in the past. In particular, since the translucency effect becomes significant for
many materials in near infrared (NIR) wavelengths, infrared
photography is used as one of common techniques for thispurpose. For example, it has been used for observing inner
layers of oil paintings that tell us a lot about the drawing
technique, history, and/or authenticity of old-age painters.
One of the major difficulties in observing inner layersof translucent objects is to separate individual inner appearances with properly dealing with scattering. To overcome this difficulty, we develop a multi-frequency illumination method, which can recover sharp appearance of innerslices at a desired depth with explicitly removing scattering
blur. Compared with conventional techniques that aim ata similar goal, our method is faster and safer than the Xray fluorescence technique [1], and sharper results can beobtained differently than infrared reflectography [2].
Our method exploits the fact that the spread of light
due to scattering has dependency on the depth of inner
layer from which light rays are emitted. By modeling thelight spreads as depth-dependent point spread functions•••

K. Tanaka and Y. Yagi are with the Institute of Scientific and Industrial
Research, Osaka University, Japan,
E-mail: {tanaka, yagi}@am.sanken.osaka-u.ac.jp,
Y. Mukaigawa and H. Kubo are with the Graduate School of Information
Science, Nara Institute of Science and Technology, Japan,
E-mail: {mukaigawa, hkubo}@is.naist.jp,
Y. Matsushita is with the Graduate School of Information Science andTechnology, Osaka University, Japan,
E-mail: yasumat@ist.osaka-u.ac.jp

Manuscript submitted November 30, 2015.

(PSFs) and exploiting their depth-dependent low-pass characteristics, we develop a method for recovering inner layer
appearances from a set of images taken under a variantof high-frequency illumination [3]. Specifically, our method
uses a spatial pattern projection with varying its pattern
pitch – we call it multi-frequency illumination. Our multifrequency illumination method allows us to separate direct(high-frequency) and global (low-frequency) components ina similar manner to [3], yet at various frequency levels thatdefine high- and low-frequencies. The observed direction
components are then used for recovering the appearance of
inner slices, which are related to the direct components via
depth-dependent PSFs.
Furthermore, to deal with inhomogeneous translucent
objects that exhibit spatial PSF variations due to its inhomogeneity, we develop an extended method for determining
the pixel-wise PSFs with an assumption of spatial smoothness of the target depth. In addition, to reduce the visible
shadowing effect that appears in the inner layers, which isdue to occlusion by their upper layers, we develop a shadow
detection method so that the shadowed regions can be later
post-processed for better visual quality.
The key contributions of this paper are threefold. First,
we describe the relationship between depth inside a translucent object and its PSF by a physically motivated scatteringmodel. Second, based on the relationship, we develop amethod for recovering the appearance of inner slices using varying pitch pattern projection. Third, we develop amethod for recovering clear inner slices of inhomogeneous
objects based on the analysis of spatially non-uniform transmittance. We apply our method to the real-world scenes
using a coaxial projector-camera system, and show the effectiveness of the proposed method using oil painting, mural,and paper scenes.
This paper extends its preliminary version [4] with thefollowing major differences: (1) We provide more discussions and details about the PSF of subsurface scatteringIEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

blur and its relationship to depth inside the media andits translucency. (2) The proposed technique is further generalized for handling spatially inhomogeneous translucent
objects. (3) A post-processing method is developed for detecting the shadowing effect due to occlusions.
The rest of the paper is organized as follows. Section 2
reviews related prior works. Section 3 describes the image
formation model for translucent objects and the depthdependent PSF model that represents blur effects due tosubsurface scattering. Based on the model, we describe theproposed method for recovering inner slices of translucent
objects in Sec. 4. Section 5 shows results for both simulation
and real-world data. Finally, we conclude the paper in 6.

2

R ELATED WORK

Our work aims at recovering the appearance of inner slices
from superposed observation. There are threads of works
for layer decomposition from super-imposed images [5], [6],[7]. Since these methods are designed for semi-transparent
layers that do not exhibit multiple scattering, the image formation model is represented by simple alpha blending [8].
Our method, on the other hand, is developed for dealing
with translucent scenes where significant multiple scattering is observed, where alpha blending cannot describe theimage formation.
In computer graphics area, various image formation
models for representing subsurface scattering have been
proposed. Donner et al. [9] propose a multi-pole model for
rendering translucent layered objects such as human skinby extending dipole [10] and Kubelka-Munk theory [11].
d’Eon et al. [12] propose an accurate image formation
model for translucent objects by decoupling single andmultiple scatterings and quantizing the diffusion Green’s
function [13]. Based on [12], Jakob et al. [14] propose aframework of rendering arbitrarily layered materials. Thesemodels enable high-quality rendering of translucent objects
with multiple scattering. Since they are designed for theforward rendering process, it is not straightforward to directly adopt them in our backward inverse problem, i.e.,
recovering appearance of inner layers of translucent objects,
due to the computational complexity. In this work, we use
a simplified version of Kubelka-Munk theory [11] for theimage formation model, which makes the inverse estimation
problem tractable while retaining the expression ability.
Our work is also related to imaging through scattering/occlusion methods in computational photography.
Levoy et al. [15] combine a synthetic aperture technique [16]in remote sensing and a confocal imaging technique [17] inmicroscopy to reconstruct images of scenes behind occlusion
or inside scattering medium. Fuchs et al. [18] also use theconfocal imaging for the purpose of recording a solid target
in scattering media. Because these methods are based on
confocal imaging, for which both the camera and projector
need exact focusing on the same depth plane, it is needed
to have a priori knowledge about the target depth. For
imaging through scattering media, Narasimhan et al. [19]’s
and Gu et al. [20]’s methods sharpen images of a target scene
in muddy liquid by precisely modeling single scattering.
Their methods work well for those scenes that do notexhibit multiple scattering. Differently from these works,

2

our method recovers images of inner slices of a translucent
object, where significant multiple scattering is observed andthe optical thickness of the target is unknown.
Our method can be grouped in a class of active sensing techniques that use high-frequency pattern projection.
The original Nayar et al. [3]’s method and its extended
methods [21], [22], [23], [24], [25], [26] separate direct andglobal components by projecting multiple high-frequency
patterns. Reddy et al. [27] separate light transport into direct,
near-indirect, and far-indirect rays by frequency-domain
modeling and analysis. Our method is also based upon
Nayar’s method [3] and we use a relationship among directcomponents, the size of scattering blur, and the pitch of the
projection pattern to separate depth layers. There are otherpattern projection techniques to decompose light transport.
Gupta et al. [28] acquire scene depths with direct-global
separation by modeling both projector’s defocus and global
light transport. O’Toole et al. [29], [30], [31] illuminate ascene by a pattern while masking the camera by the complementary pattern to spatially probing the light transport of
the target scene. Our method also uses a pattern projection
technique not only for separating scattering effects, butalso for specifically recovering appearance of multiple innerslices.
Time-domain coding for analyzing light transport is another approach to recover images that are not directly measurable. Heide et al. [32] sweep the modulation frequency
and phase of their customized Time-of-Flight (ToF) camera
to recover the light propagation inside scattering medium.
Kadambi et al. [33] build a coded-illumination ToF camera
with a deconvolution technique and use it for recoveringa sharp image by observing through a diffuser. O’Toole etal. [34] combine spatial probing and ToF imaging to separate
direct and indirect light-in-flight images. Tadano et al. [35]propose an imaging system that is capable to select a target
depth using a coded ToF camera. While these methods effectively recover light transport, they require carefully tailored
ToF cameras. Contrary, our method uses a simple projectorcamera system and spatial pattern coding to analyze lighttransport inside translucent objects.
Methods for measuring transparent or translucent media, such as smoke scanning, are also related to our work.
Morris et al. [36] recover the shape of clear transparent
objects that refract light by recording light rays from different viewpoints. Hawkins et al. [37] acquire the density
distribution of participating media, such as smoke, by laser
scanning. Ihrke and Magnor [38] reconstruct the volume
of dynamic semi-transparent media by a visual-hull based
method. While related, our goal is to recover slices insidetranslucent objects instead of reconstructing 3-D shape ofan object’s surface.
There have been independent developments of technologies for imaging internal structures of target objects for special purposes. In art analysis, several techniques have beendeveloped for imaging hidden layers of paintings. Infrared
reflectography [2] and X-ray transmission radiography [1]have been used for visualizing internal layers of paintings,
although the surface texture cannot be separated. X-ray
fluorescence technique [1] uses spectroscopic information
measured over a couple of days and estimates the metallic
atom distribution for determining colored appearance ofIEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

PSFs

3

as illustrated in Fig. 1(b). In this manner, there is a close
relationship between the PSF and depth. By denoting ht asa PSF at depth t, the appearance slice St can be expressed as

St (c) = (Rt ∗ ht )(c),Radiance at each depth
(a)

(b)

Fig. 1: Illustration of the image formation model for translucent objects. (a) Recorded intensity is the summation of alllayer’s appearance. (b) Spatial spread of light varies withdepth.
inner layers. Tera-hertz imaging [39] is another techniquethat can see inner layers of paints. In the medical imaging
and its related areas, optical coherence tomography [40],[41] techniques are widely used, especially for visualizing
retina. These techniques enable to observe inside translucent
objects based on interference of coherent light at the costof expensive wavelength-order optics and mechanics. In
contrast, our method uses a commodity camera and projector for recovering slices inside translucent objects, whichallows low-cost implementation. In microscopy, there aremethods that use pattern projection for visualizing insidesubstances, or reconstructing 3-D shape of a small sample,
such as protein structure [42], [43], [44]. They sharpen the
microscopic image by precisely taking into account both
micro-scale optical blur and scattering blur. Our aim is todevelop a technique that is applicable to a standard scale,where scattering blur becomes more significant.

3

A PPEARANCE OF TRANSLUCENT OBJECTS

When an image of a translucent object is recorded, the observed intensity Lo (c) at camera pixel c ∈ Z2 is a summationof all possible paths ρ ∈ Pt from depth t:
Z τZ
Lo (c) =ι(ρ)dρdt,
(1)
Pt

0

where τ is the maximum thickness of the object, and ι(ρ)
is the intensity of the light path ρ. It can be modeled
as a summation of the appearance of all depth slices as
illustrated in Fig. 1(a), which can be described asZ τ
Lo (c) =St (c)dt,(2)0where St is the appearance slice at depth t.
The appearance slice St is generally blurry due to thescattering effect inside the medium. The spread of radianceat a scene point inside a translucent object varies depending on its depth from the object surface [45]. In general,
the spatial spread of light can be expressed using PSFs.
Let us consider light rays emitted (or returned) from a
specific depth inside a translucent object. When the depth
is shallower, the PSF becomes sharper. It gradually wider
spreads as the depth t becomes deeper inside the medium(3)

where Rt is the sharp appearance slice that we are interested
in estimating, which we call a radiance slice, and ∗ denotes
a convolution operator.
Since the appearance of the translucent object under normal illumination is a superposition of radiance of multiple
slices as Eq. (2), the observation Lo can be re-written asZ τ
Lo (c) =(Rt ∗ ht )(c)dt.
(4)0Specifically, we are interested in recovering a few planar
parallel layers inside the medium; therefore, the image
formation model can be discretized as a multiple layered
model asX
Lo (c) =(Rd ∗ hd )(c),(5)
dwhere Rd and hd are the radiance slice and depthdependent PSF of d-th layer, respectively. Our goal is recovering radiance slices Rd convolved by unknown depthdependent PSFs hd from the composite observation Lo . Before introducing the solution method, we describe a modelof depth-dependent PSFs hd .
Depth-dependent PSFs
Our depth-dependent PSF model is motivated by the workof [46], which is a physically motivated scattering model
designed for the PSF of scattering medium, named theradiative transfer equation (RTE) [47], [48]. In this work, we
represent the depth-dependent PSFs hd by the RTE model.
With the RTE model, the intensity of an arbitrary light ray
in a homogeneous scattering medium can be iteratively
calculated [46] as

I(d, θ) =∞X(gk (T ) + gk+1 (T )) Lk (cos θ),

(6)

k=0

where Lk is the Legendre polynomial of order k , g0 = 0,and
2k + 1k−1(1 − q)T − (k + 1) log T .
gk (T ) = I0 exp −k
The parameter T (= σd) represents optical thickness, whichis the product of scattering coefficient σ and distance d
between the point light source and the scattering point.
The forward scattering parameter q controls how light rays
spread; q = 0 corresponds to isotropic, while positive andnegative values of q indicate forward and back scattering,
respectively. θ is the angle between the depth axis and lightray’s direction as depicted in Fig. 2(a). k is the numberof light bounces, and I0 is the intensity of the point lightsource.
In our setting, a camera is placed outside the media,
where light rays do not scatter. In this case, the depthdependent PSF applied to a point light source inside scattering media corresponds to the intensity of the light rays
emitted from the surface, if the camera is placed sufficiently
far from the object. When we consider that the x-axis lies

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

scattering point

d

θ

x

light sourced

Homogeneous medium
(a)

φ d
I0
(b)

Intensity

Fig. 2: Illustration of RTE model described in [46]. (a) Solidarrow represents a scattering light ray, whose intensity
depends on d and θ. d is the distance from the point lightsource I0 to the point inside the scattering medium, and θ isthe radial direction of the light ray. (b) Depth-dependent PSF
hd can be expressed using RTE model. We consider the xaxis lies along the tangent plane of the surface of translucent
object and assume the direction of light rays emitted fromthe surface becomes parallel to the depth axis.

1.0
0.75
0.5
0.25
0.0

-50

P ROPOSED M ETHOD

We are interested in recovering radiance slices Rd from themixed observation Lo . To achieve this goal, we develop amulti-frequency illumination measurement method, which isbuilt upon the high-frequency illumination (HFI) method
proposed by Nayar et al. [3]. To begin with, we briefly
review the original HFI method.

hd (x)I(d, θ)I0

4

camera4

T =1.12

T =1.82

T =6.420

50 -500
50 -50
Radial length [px]050Fig. 3: Simulated PSF variations using the RTE model withvarying optical thickness T

along the tangent plane of the surface of translucent object
as illustrated in Fig. 2(b), the depth-dependent PSF hd (x)can be written as

hd (x) = I (d0 , φ) ,(7)√
where d0 = x2 + d2 and φ = tan−1 (x/d). We assume
that the camera is placed sufficiently far from the object
compared to both d and x so that the direction of emitting
light rays from the surface becomes parallel to the depth
axis. In addition, we ignore the refraction of the surface forsimplicity. Figure 3 shows simulated PSFs by changing theoptical thickness T . When T is small (thin scattering media),a sharply pointed PSF is generated, and as T becomes larger(thicker scattering media), the PSF spatially spreads wider.
These PSFs only represent the spread of light but dot
not consider the absorption by the medium. When theabsorption is homogeneous, it only makes the deeper layers
darker and has no effect on the sharpness of the layer
recovery except for the siganl-to-noise ratio (SNR), hence
it can be practically neglected. However, if the absorption
is inhomogeneous, which is the case in inhomogeneoustranslucent objects, the radiance slice of inner layer exhibits
artifacts if the spatial variation of PSFs is ignored. In thiswork, we explicitly take into account the spatially-varying
PSFs in Sec. 4.3 for dealing with inhomogeneous translucent
objects.

High-frequency illumination method [3]The HFI method separates direct and global components byprojecting small pitch checker patterns. When the phase of
the projection pattern changes slightly, the direct componentD(c) varies accordingly, but the global component G(c)remains stable. Based on this observation, their method
computes direct and global components using the maximum Lmax (c) and minimum Lmin (c) intensities that are
obtained by shifting the projector pattern asD(c) = Lmax (c) − Lmin (c),(8)G(c) = 2Lmin (c).
The direct component D(c) contains high-frequency components, while the global component G(c) contains only lower
frequency components than the frequency of projection
pattern. Therefore, the HFI method can be viewed as a separation technique for high- and low-frequency components.
Pattern pitch of HFI
In our case, when a translucent object is measured under
HFI with pattern pitch p, we can obtain direct componentD(p, c) and global component G(p, c) asP D(p, c) =Pd Dd (p, c)G(p, c) =(9)d Gd (p, c)
Sd (c) = Dd (p, c) + Gd (p, c),where Dd (p, c) and Gd (p, c) are the direct (high-frequency)
and global (low-frequency) components at depth d, respectively, and the sum of direct and global components for each
depth becomes the radiance slice. As mentioned in [3], thepattern pitch p must be sufficiently smaller than the scenetexture for a faithful separation.
We have observed that the direct and global components
vary when the projected pattern pitch p changes. Indeed,
there is a tight relationship between the pitch size andstrength of direct component. Suppose we measure thescene with two distinct pattern pitches pi and pj (pi < pj )independently. When the pattern pitch becomes larger, thewider spread of light rays are included in the direct component, hence the direct components have the following
relationship:

Dd (pi , c) < Dd (pj , c),pi < pj .

(10)This relationship indicates that the separation frequency
varies with the pattern pitch – wider frequency band is regarded as the direct component as the pattern pitch becomeslarger. As low-pass characteristics of depth-dependent PSFs
are also depth-dependent, the difference of direct components Dd (pj , c) − Dd (pi , c) also varies with depth d. Our
method exploits these differences for recovering the appearance of each inner slice by changing the projection pattern
pitch.

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

4.1

Multi-frequency illumination

By measuring the target object using multi-frequency patterns, multiple of corresponding direct components are obtained. Unfortunately, increasing the number of measurements does not make the problem easier as it also increases
the number of variables to solve for. To make the problem
tractable, we assume that the texture of direct components
does not vary drastically when the pattern frequency is highenough and the pitch variation (pj − pi ) is sufficiently small.
These direct components Dd at a certain depth d are supposed to have a similar texture with the original radiance;
therefore, we can expect the following relationship:

Dd (p, c) ≈ α(hd , p)Rd (c),(11)

where α(hd , p) is the relative brightness of Dd (p, c) to Rd (c).
We call α(hd , p) the direct component ratio that represents theratio of direct component’s mean intensity to the radiance
Rd (c)’s mean intensity. Hence, Eq. (9) can be rewritten asX
D(p, c) =α(hd , p)Rd (c).
(12)
dThese assumptions are based on the fact that the diffuse
reflection and subsurface scattering can be regarded as
the same physical phenomena [49], [50]; the light scatterson or beneath the surface and eventually bounces off ofthe material in random directions. Direct components Dd
represent total intensities of lights from all the points inside
the object whose distance from the incident point on thesurface is smaller than the pattern pitch p. Hence the pattern
pitch p controls the scale of the separation of scattered lights
in the direct-global separation scheme, and thus controls
the intensity of direct components. Furthermore, because
p is sufficiently smaller than the scene texture, the textureof direct components and the original texture are largely
similar. Based on these observations, we obtain the original
texture at different brightnesses by changing p.
With these assumptions, a set of direct component images D(p, c) taken under the multi-frequency illuminationof m pitch variations (p = p1 , p2 , . . . , pm ) can be written ina matrix form as

D(c) = AR(c),

(13)

where

D(c)

=

A =



D(p1 , c) D(p2 , c) · · · D(pm , c)α(hd1 , p1 ) · · · α(hdn , p1 )


..
..
..

,
.
.
.
α(hd1 , pm ) · · ·

R(c)

=


m

T

,

α(hdn , pm )Rd1 (c) Rd2 (c) · · ·

Rdn (c)T.

Here, D ∈ R is a vector of direct components measured
under m variations of the pattern pitches at pixel c, R ∈ Rn
is a vector of n layers of radiance slices, and A ∈ Rm×n is
a matrix containing direct component ratios computed fromthe projected pattern pitch and the depth-dependent PSF.
When the number of projected patterns m is no less than
the number of depth layers n (m ≥ n) and rank(A) = n, theradiance slices R can be obtained by a norm approximationof the residual vector, i.e., D(c) − AR(c). For example, with5

a least-squares approximation, the radiance slices R(c) canbe determined using the pseudo-inverse A+ as

R(c) = A+ D(c).

(14)

Computation of direct component ratio
The direct component ratio α(hd , p) can be derived from thedepth-dependent PSF hd and the projected pattern pitch p.
When a checker pattern is projected to a translucent object,
it reaches the depth d with some blur effect and returns
to the surface with the additional blur effect. To obtain the
direct component ratio, we consider the difference between
maximum and minimum intensities in a similar manner tothe original HFI as :α(hd , p) = max ((lp ∗ hd ) ∗ hd )− min ((lp ∗ hd ) ∗ hd ) ,(15)

where lp is normalized projection patterns, whose pitch is p.
The normalized illumination is defined in the range between
0 and 1; 0 being black and 1 being white. The operators max
and min return the maximum and minimum value from all
pixels, respectively.
4.2

Estimation of informative slices

Once we know the target depths (or PSFs) to inspect, we
can set up a matrix A using Eq. (15), and thus can recoverslices corresponding to the depths using Eq. (14). However,such a prior knowledge is difficult to obtain before measurement; therefore, automatically selecting a good set of
depths becomes important for recovering informative slices.
For example, if an arbitrary depth is chosen, it has a chance
to correspond to the middle of distinct texture layers. To
recover informative slices, we use a two-step approach. The
first step is the estimation of a set of informative depths via
optimization. This is equivalent to selecting a small numberof useful PSFs from many other possible PSFs. The second
step is the recovery of slices at the informative depths
determined by the earlier step. The overall procedure of thisstrategy is illustrated in Fig. 4. Now, we explain the detailsof each step.
Step 1: Estimation of informative depths
Elements of matrix A depend on depth-dependent PSFs;
hence estimation of informative depths corresponds to estimating the shape of PSFs. The candidate relative depth T ,which determines the shape of the PSFs, are set in the range
between 1 +  and 10 +  with a step size 0.05. We add
the offset  (= 0.01) because Eq. (6) does not converge if
T ≤ 1. Initially, we set up the matrix A using all candidate
parameters T (thus m < n) in order to estimate informativeslices.
Frequently, there are only a small number of informativeslices inside translucent objects. Hence we can regard suchradiance slices exist sparsely along depth. Our method uses
this sparsity to determine the informative slices by solving
a l1 regularized problem (as known as the lasso [51]) with a
non-negative constraint about R:
2

R̂(c) = argmin kAR(c) − D(c)k2 + λ kR(c)k1
R(c)

subject to R  0.

(16)

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

6

Inhomogeneous upper layer (different PSFs)

Pigments
Different irradiance
Inner layer
(a)

(b)

Fig. 5: Scene of inhomogeneous upper layer. PSFs and irradiance of inner layer vary depending on the upper layer’s
optical property.

(c)Fig. 4: Selection of informative PSFs (corresponding to certain depths). (a) Estimates R̂d obtained via optimization(Eq. (16)). Non-zero pixels indicate the informative regions.
(b) We find the local maxima of the non-zero pixel counts(corresponding to l0 -norm of R̂d ) for all pixel coordinatesacross depth. (c) Selected PSFs corresponding to the local
maximasWe can regard the depth d (= Tσ ), where R̂d (c) has a nonzero value, is informative while others are not. Equation (16)becomes a quadratic programming (QP) problem and thus
can be efficiently solved in a polynomial time. We solve
the optimization in a per-pixel manner. Solving a similar
problem for the entire image at a time instead of computing
in a pixel-wise manner is also a viable option; however,
we have observed that they do not make much difference
because of the following Step 2. Therefore, for efficientparallelization, we choose the per-pixel implementation.
Step 2: Informative slice recovery
This step determines informative depth slices of the whole
image by consolidating all the pixel-wise selections. The
informative depths d̂ are local maximas of the sum of l0
norm of R̂d (c) for all pixels. The evaluation function f (d) is
defined asX
kR̂d (c)k0 ,(17)
f (d) =
cand we find all local maxima of f (d) along d in the rage ofinterest as shown in Fig. 4(b). Once the depths of interest
are selected, we can set up a small matrix Â. Finally, theappearances of informative slices are recovered using thematrix Â using the least-squares approximation as

R(c) = Â+ D(c)(18)

in the same manner to Eq. (14). We call the method describedso far a baseline method.
4.3 Recovering inner slices inside inhomogeneous upper layer
Now we describe an extension of the baseline methodfor dealing with inhomogeneous translucent objects. Whenthe upper layer is inhomogeneous, due to either different
materials or thicknesses, or even both, recovered inner slices
by the baseline method are affected by the upper layer’s
non-uniform scattering and transmission. Because of thespatially non-uniform scattering and absorption property,
both PSF and irradiance of the inner layer vary even if theinner layer is homogeneous as illustrated in Fig. 5. Under
this situation, the baseline method suffers from artifacts
because it assumes a single PSF per slice for the entire image
region.
To deal with the inhomogeneous upper layer, we employ
a pixel-wise PSF selection approach. By applying Step 2of the baseline method in a sliding window manner, a setof PSFs can be obtained for each pixel. We apply a 2D
Gaussian-weighted filter for the sliding window for suppressing high-frequency observation noise. The evaluation
function fc (d) for pixel location c is defined asX(19)
fc (d) =Nσf (kcn − ck2 ) R̂d (cn ) ,cn ∈Wc0where Nσf is a zero-mean Gaussian distribution with standard deviation σf centered at c, and cn is a pixel locationin the window Wc . By selecting the informative depths foreach pixel via finding the local maxima of fc , we can recoverpixel-wise informative slices in a similar manner to Eq. (18)as

R(c) = Â+ (c)D(c),(20)

where Â+ (c) is the pseudo-inverse of a small matrix containing direct component ratios computed for pixel c using
the selected PSFs.
The recovered inner slices may still contain the brightness discontinuity due to the non-uniform transmission,
because PSFs considered so far neglect absorption by themedium. Inhomogeneous transmittance causes inhomogeneous irradiance on the inner layer, hence the recovered
inner layer may have brightness variations. Indeed, thisbrightness discontinuity stands out at the boarder of distinct
PSFs. In this work, we adopt a gradient-based image repairing approach to reduce the brightness discontinuity. Specifically, we suppress the gradient of edges at the brightness
discontinuity in the inner with retaining other edges so thata cleaner appearance of the inner slice can be obtained. Suchedges E can be identified by observing the change of PSF
estimates because the discontinuity of transmission occurs
at the boarder of different materials or thicknesses in theupper layer. Hence, the brightness corrected slice Rd0 can berecovered by setting the gradient of such edges to zero as(0
if c ∈ E
0∇Rd (c) =(21)∇Rd (c) otherwise,&URVVFRUUHODWLRQ

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

1.000
0.995
0.990
0.985
0.980
T =1.08
T =3.240.975
0.970
T =2.16
T =4.860.965
0.960
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

3DWWHUQVL]HpSL[HOV

Fig. 6: Evaluation of the approximation in Eq. (11). We
calculate the cross correlation between direct component
slice Dd (p) and the ground truth slice Rd . It shows highcorrelation even for the worst case.

where ∇ is the gradient operator and E is the set of edges
corresponding to the material (or thickness) borders. The
integration is computed by a standard sparse linear systemsolver as done in the work of Poisson image editing [52].

5

E XPERIMENTS

We evaluate our method numerically by simulation and also
show qualitative results using real-world scenes.
5.1

Evaluation by simulation

We first assess the appropriateness of the approximationin Eq. (11), and then evaluate the accuracy of the slice
recovery.
Validity of approximation in Eq. (11)In this simulation, we change the pattern pitch p from 3
to 20 pixels, the depth d from 4 to 18 [mm] (in opticaldepths T , which correspond to 4σ to 18σ ). The scene is onelayered for the purpose of assessing the approximation in
Eq. (11). The coefficient σ is set to 0.27 [mm−1 ], and forward
scattering parameter q is set to 0.9. These parameters arechosen according to [53]. We generate the appearance slices
ϕSd,pusing checker pattern lpϕ as
ϕSd,p= ((lpϕ ∗ hd ) ◦ Rd ) ∗ hd .

(22)p and ϕ are the pitch and phase of the pattern, and ◦ is the
Hadamard (element-wise) product operator. We compute
direct component Dd (p) from these synthetic images withchanging the phase ϕ.
We assess correlations between the recovered direct components Dd (p) and the ground truth slices Rd using zeromean cross correlation (ZNCC). The ZNCC value falls in
the range from −1 to 1, where a negative value indicatesnegative correlation, and a greater value indicates higher
similarity. The evaluation results are summarized in Fig. 6.
The plots correspond to the average of ZNCC scores obtained from 15 different textures Rd . ZNCC values decreases
as the depth becomes deeper (greater T ) and pattern pitch
becomes larger. However, they are all consistently highly
correlated (minimum ZNCC value is 0.982 and the meanvalue is 0.992), and it shows the accuracy of the approximation.

7

LayersTop layer
Bottom layer

ZNCC values
Max. Mean Min.
0.99
0.93
0.89
0.92
0.84
0.67

TABLE 1: Two layers recovery result for 20 sets of scenes.
We compare recovered slices with the ground truth slices
by ZNCC. Higher ZNCC scores indicate more accurate
recovery.

Evaluation of slice recovery
Here we evaluate the overall accuracy of our method via
simulation. We generate images Lϕ
p under projection patterns lpϕ as
X
Lϕ((lpϕ ∗ hd ) ◦ Rd ) ∗ hd .
(23)
p =
d

We use 20 different two-layered scenes that have distinct
textures R5 and R15 at two depth ranges 5 ≤ d < 15 and15 ≤ d. The textures are randomly paired from the ones
that are used in the previous experiment. We change thepitch p of checker pattern and shift the pattern with ϕ. Fromthe generated images, we compute D(p, c) for each pattern
pitch p, and apply our method to recover slices at d = 5 and15.
We again use ZNCC values between the recovered andground truth slices. The experimental results are summarized in Table 1. In all data sets, the ZNCC values of upper
layer is higher than the lower layer as expected. The result of
this simulation indicates that our method can recover slices
of various textures with high accuracy.
Figure 7 shows a synthetic example of three-layer recovery, where slices at at d = 1, 5, and 15 are recovered. The
ZNCC scores for the recovery results are 0.98, 0.83, and0.51, respectively. While the result is generally consistent
with the ground truth, negative intensities and ringing
artifacts are observed due to the discrete pitch variations
and convolution.
5.2

Real-world experiment

We develop a coaxial projector-camera setup for realizing
the measurement setup as shown in Fig. 8. The coaxial setup
has a favorable property; correspondence between projector
and camera pixels becomes invariant with respect to depths.
Unlike non-coaxial settings, with which a illumination ray
inside the translucent object forms a line in the image coordinates [24], the coaxial setting allows us to easily separate
the direct rays. We use a LightCommander projector, which
is a DMD projector development kit by Texas Instruments,
and use near infrared (NIR) light for measurements. The
lenses of both camera and projector are set equivalent (AiMicro-Nikkor 105mm f/2.8S) for making the alignment easy.
In the experiment, we use 18 variations of checker patterns(3px to 20px with 1px interval), and shift the pattern for
one-third of square size for each pattern.
Experimental results
First, we use an oil painting as a target scene as shown in
Fig. 9(a), which has draft and surface pigment layers as

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY(a) Simulated scene. three star-shaped pillars with different
heights are placed in scattering media. An example slice isillustrated by the green plane.
Ground truth slices

8(a) The scene(b) Draft layer(c) Oil painted

Recovered slices

d=1
d=5
d = 15(b) Comparison of the recovered and ground truth slices(d) NIR photo(f) Li and Brown [7](h) Ours(reflectography)(sharp layer)(surface slice)(e) Enhanced(g) Li and Brown [7]Fig. 7: (a) Target scene (b) Result of three-layer recovery at
d = 1, 5, and 15. The ground truth radiances slices (upper)
and recovered slices (lower) are shown. ZNCC scores are
0.98, 0.83, and 0.51, respectively.

cameraa pair of same lenses

beam splitter

projector

target object
Fig. 8: Measurement setup. The coaxial system allows us
to maintain the correspondences between projector andcamera pixels.

depicted in Figs. 9(b) and 9(c). By taking a standard photograph under the near infrared light (as done in infrared
reflectography [2] in the art field), we can only vaguely
observe the draft layer as shown in Fig. 9(d). Since it is thesuperposition of draft and surface pigment layers, it naturally results in a blurry image. Even with a manual contrast
adjustment, it is difficult to clearly observe the shape of draft
tree as depicted in Fig. 9(e). A simple layer decomposition
still suffers from the blur as shown in Figs. 9(f) and 9(g) . By
applying our method to this scene, two PSFs are estimatedas depicted in Fig. 10(a) and two slices are recovered asin Figs. 9(h) and 9(i). The upper surface layer corresponds
to the surface pigment layer. Because the yellow pigment
is almost transparent in the infrared wavelength, the corre-(blurred layer)(i) Ours(inner slice)

Fig. 9: Experimental result of oil painting using the baseline
method. (a) Target scene. We draw a colored round tree
on top of the draft of spiny tree. (b) Inner layer (draft) of
the painting. (c) Painted scene. Red rectangle region is measured. (d) Normal photo using infrared light. (e) Enhanced
image. Intensity range and contrast of (d) are manually
adjusted. (f, g) Layer separation results of Li and Brown [7].
Because their method separates sharp and blurred layers, it
suffers from global components. (h, i) Results of our method.
Layer of surface texture and hidden drawing, respectively.
Range of the intensities is adjusted for visualization.

sponding painting regions become dark in the surface slice.
The lower layer shows the inner layer, where the texture of
the tree is clearly observed. The variation of the remaining
high-frequency components of each layer, which can becomputed from its PSF using Eq. (15), along with the pattern
pitch is shown in Fig. 10(b). It shows that the high-frequency
components of the upper layer remain more significantly
than that of the inner layer, and it indicates the separabilityof the layers.
The second target scene is stacked translucent sheets thatconsist of layered translucent sheets and printed texture
films as shown in Fig. 11(a). With a conventional NIR
photograph, we can observe mixed textures as shown in
Fig. 11(b), where textures are blurry and appearances acrossdepths are superposed. Our method correctly selects two

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

9l0 -norm

8e4

two selected PSFs

6e4
4e4
2e4
0
1.0

3.0

5.0

7.0

9.0

optical depth ( T )

11.0

13.0

15.0

Remaining high-freq.

(a) Selected PSFs corresponding to local maxima

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

(a) The scene(b) NIR photo(c) Upper slice(d) Lower slice

Upper layer
Inner layer

4

6

8

10
12
14
Pattern pitch [px]

16

18

20

(b) Remaining high-frequency components for each layer

Fig. 10: Selected PSFs and remaining high-frequency components for each layer. (a) Selected PSFs. There are two
peaks in the plot of Eq. (17), hence two corresponding
PSFs are selected to recover. (b) Remaining high-frequency
components for each layer. Upper layer remains more highfrequency components than that of the inner layer.

informative depths and recovers their radiance slices asshown in Figs. 11(c) and 11(d). In the upper slice, only
‘ABCD’ texture is visible, and ‘1234’ texture appears in
the lower slice. Due to the shadowing effect caused by
the opaque material on the upper layer (the texture of‘ABCD’), the lower slice contains the ‘ABCD’ texture as anunobserved shadowed region.
Additional results are shown in Fig. 12. The top row isa piece of painting, where the painter’s signature is hidden
under pigment. Our method clearly recovers surface texture
and inner signature slices. The middle row is a mural
painting covered by white mold. In this example, we used
RGB light sources instead of NIR, and record each color
channel separately. Our method can recover the slices inthis example as well, and the inner appearance is clearlyvisible in the result of Fig. 12(d). The bottom row is a doublesided ancient document with RGB light sources. While anordinary digital scanner suffers from the back-side texture,
our method can faithfully separate surface and back-side
textures.
Figure 13 shows the result of the scene with an inhomogeneous upper layer. We use a 21 × 21 pixels Gaussian window for selecting PSFs. The target object is an oil painting,where the upper layer consists of a wide variety of pigmentsin their thicknesses and materials. The result of the baseline
method shown in Fig. 13(b) exhibits brightness discontinuity at the middle of image region due to inhomogeneity of
the upper pigments. By selecting PSFs per pixel using theGaussian-window, the inner layer is recovered in a pixelwise manner as shown in Fig. 13(f), although it still suffersfrom brightness discontinuity. Finally, with the gradient-Fig. 11: Experimental result of layered scene using thebaseline method. (a) Scene is a composite object of texture
and translucent sheets (tracing paper). (b) Normal photo of
the scene. Textures inside the object can be seen, althoughblurry. (c) Recovered slice of upper layer. (d) Recovered
slice of lower layer. Textures in sub-millimeter gap can beseparated. Further analysis of shadowed region is shown in
Fig. 14domain filtering to suppress the brightness discontinuity,
the inner layer that has uniform brightness can be recovered
as shown in Fig. 13(h). It shows improvement in the visualquality, especially in the left part of the image and overall
contrast.
Shadow detection from recovered inner layers
The recovered appearance of inner layers suffers from shadows caused by their upper layers. If the upper layer’stransmittance is low, it casts shadow to the inner layer, e.g.,
as ‘ABCD’ scene shown in Fig. 11(d). Here we describe
a simple technique for dealing with shadows by postprocessing. To simplify the discussion, we assume that there
are only two layers in the scene, i.e., upper and inner layers.
The shadowed regions have two common properties:(1) The shape of shadows observed in the inner slice issimilar to the texture of upper slice. (2) The shadowed
region becomes darker due to low irradiance. Based on theseobservations, we define a shadow likelihood measure Ps for
identifying shadowed regions.
The similarity S of the texture shapes in the recovered
slices R1 and R2 can be obtained using the absolute valueof their cross correlation as

S (c) = |Cw,c (R1 , R2 )| ,(24)

where Cw,c is the cross-correlation within a small window
w centered at pixel location c, and values in R1 and R2 arenormalized in the range of [0, 1]. Using this similarity, theshadow likelihood Ps is defined using the darkness of the
shadowed region in R2 as

Ps (c) = S (c)(1 − R2 (c)).

(25)

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY(a) Target scene(b) Normal photo(c) Surface slice

10(d) Inner slice

Fig. 12: Additional results for a painting, a mural, and an ancient document using the baseline method. (a) Target scenes.
Top: author’s signature covered by red pigment. Middle: A mural painting covered by white mold. Bottom: An ancient
document of double side print. Rectangle regions are measured. (b) Normal photograph. Both upper and lower slices, andglobal components are composed. (c) Recovered slices of surface. (d) Recovered inner slices.

10σ

100σ

(a) Normal IR photo(b) Inner slice (baseline)(c) Upper layer’s depth(d) Inner layer’s depth(e) Upper slice (pixel-wise)(f) Inner slice (pixel-wise)(g) Gradient mask(h) Inner slice (post-processed)

Fig. 13: Result of a scene with an inhomogeneous upper layer using the extended method. (a) Captured image usingnear infrared lights. (b) Recovered inner slice by the baseline method. (c, d) Selected depths for upper and inner layers,respectively. (e) Recovered upper slice by pixel-wise recovery. (f) Recovered inner slice by pixel-wise recovery. (g) Gradientmask computed from (d) that is used for gradient-domain filtering. (h) Recovered inner slice from (f) and (g) by image
enhancement in the gradient domain.

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY(a) Inner slice
including shadow(b) Likelihood of shadow

11(c) Shadow region mask(d) Inpainted inner layer

It yields a likelihood score for each pixel being in shadow,
and by a simple thresholding, a shadowed region can be
determined.
Once the shadow regions are identified, we can use ashadow removal technique. For example, the lost information within the shadow region can be filled in by arbitrary
image inpainting methods. In this paper, we use a patchmatch based image inpainting method [54].
The result of shadow detection and removal for the‘ABCD’ and ancient document scenes are shown in Fig. 14.
For this experiment, we use 21×21 window size for computing the similarity measure. The shadow likelihood shownin Fig. 14(b) qualitatively corresponds to the shadowed
regions in the inner layer. The shadow mask generated bythresholding the shadow likelihood is shown in Fig. 14(c).
Figure 14(d) shows the image inpainting result using thedetected shadow mask and the method of [54]. It shows
improvement of visual quality of the inner layer.
Finally, we show a recovery and shadow removal result
for a three-layer scene. Figure 15(a) shows the target scene,
where three printed papers are overlaid. By applying ourbaseline method, three slices with shadows are recovered.
With shadow detection and inpainting, we obtain the recovery of three layers as shown in Figs. 15(b)-15(d) For
the shadow detection of the third slice, we compute thecorrelation with the original second slice which contains
the shadows casted by the top slice. Figure 15(e) shows
how high-frequency components remain along with varying
pattern pitches. The lowest layer only retains a limited
amount of high-frequency components hence its recovery
becomes slightly blurry. However, it shows the different
characteristics among layers, which indicates the separability of the layers.

6

C ONCLUSION

This paper described a method for recovering inner slicesof translucent objects based on multi-frequency pattern pro-

Remaining high-freq.

Fig. 14: Shadow detection and inpainting results. The scene are the same as Fig. 11 and Fig. 12. (a) The inner layer
including shadowed region. (b) Shadow likelihood from recovered slices. White pixels indicate higher likelihood of being
in shadow. (c) Shadow region is indicated by red by thresholding the likelihood. (d) Shadow region is recovered by image
inpainting [54]. Shadow artifact caused by the upper layer is suppressed.

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

(a) The scene(b) Uppermost slice(c) Middle slice(d) Lowest slice

UppermostMiddleLowest

4

6

8

10
12
14
Pattern pitch [px]

16

18

20

(e) Remaining high-frequency components

Fig. 15: The result of seeing through pages. (a) The target
object. Three printed thin papers are superposed. (b) Recovered uppermost slice. Only the first page can be seen.
(c) Recovered second page. Shadows from the upper layer
is inpainted. (d) Recovered third page. Shadows from thefirst and second layers are inpainted. (e) Remaining highfreqency components for each pixel.

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY

jection. The proposed method is built upon the observationthat the PSF inside translucent objects varies according to
the depth of slices. Based on that, we have shown that innerradiance slices can be recovered by estimating PSFs using
varying pitches of projection patterns. We also developed
a method for automatically selecting informative slices via
a sparse representation, i.e., determining sparse coefficientsthat corresponds to radiance slices. We further extended themethod for dealing with inhomogeneous translucent objects
based on a combination of pixel-wise appearance recovery
and gradient-based image repairing. The effectiveness ofthe proposed method are shown by several experiments on
simulation and real-world translucent objects.

ACKNOWLEDGMENTS
We thank all anonymous reviewers who gave us various
insightful and constructive comments. We thank Norihiko
Kawai for providing his inpainting implememtation. The
part of this work is supported by JSPS KAKENHI Grant
Number JP26-6433, Grant-in-Aid for JSPS Fellows.

R EFERENCES[1]

J. Dik, K. Janssens, G. Van der Snickt, L. van der Loeff, K. Rickers,
and M. Cotte, “Visualization of a lost painting by Vincent van
Gogh using synchrotron radiation based X-ray fluorescence elemental mapping.” Analytical chemistry, vol. 80, no. 16, pp. 6436–42,
Aug. 2008.
[2] D. Gavrilov, E. Maeva, O. Grube, X. Maldague, and R. Maev,“Infrared Methods in Noninvasive Inspection of Artwork,” inInternational Conference on NDT of Art, no. May, 2008, pp. 1–5.
[3] S. K. Nayar, G. Krishnan, M. D. Grossberg, and R. Raskar, “FastSeparation of Direct and Global Components of a Scene using
High Frequency Illumination,” in Proc. SIGGRAPH, 2006, pp. 935–
944.
[4] K. Tanaka, Y. Mukaigawa, H. Kubo, Y. Matsushita, and Y. Yagi,“Recovering Inner Slices of Translucent Objects by Multifrequency Illumination,” in Proc. IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2015, pp. 5464–5472.
[5] R. Szeliski, S. Avidan, and P. Anandan, “Layer Extraction from
Multiple Images Containing Reflections and Transparency,” inProc. IEEE Conference on Computer Vision and Pattern Recognition(CVPR), 2000, pp. 246–253.
[6] B. Sarel and M. Irani, “Separating Transparent Layers through
Layer Information Exchange,” in Proc. European Conference on
Computer Vision (ECCV), 2004.
[7] Y. Li and M. S. Brown, “Single Image Layer Separation using
Relative Smoothness,” in Proc. IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2014.
[8] T. Poter and T. Duff, “Compositing Digital Images,” Computer
Graphics, vol. 18, no. 3, pp. 253–259, 1984.
[9] C. Donner and H. W. Jensen, “Light Diffusion in Multi-Layered
Translucent Materials,” ACM Transactions on Graphics (ToG),
vol. 24, no. 3, p. 1032, jul 2005.
[10] H. W. Jensen, S. R. Marschner, M. Levoy, and P. Hanrahan, “A
Practical Model for Subsurface Light Transport,” in Proc. SIGGRAPH. ACM, 2001, pp. 511–518.
[11] P. Kubelka and F. Munk, “Ein Beitrag Zur Optik Der Farbanstriche,” Zeitschrift fr Technische Physik, vol. 12, pp. 593–601,
1931.
[12] E. d’Eon and G. Irving, “A Quantized-Diffusion Model for Rendering Translucent Materials,” in Proc. SIGGRAPH. ACM Press,
2011, p. 1.
[13] “Green Function,” in Encyclopedia of Mathematics, M. Hazewinkel,
Ed. Springer, 2001.
[14] W. Jakob, E. d’Eon, O. Jakob, and S. Marschner, “A Comprehensive
Framework for Rendering Layered Materials,” ACM Transactions
on Graphics (ToG), vol. 33, no. 4, pp. 1–14, jul 2014.
[15] M. Levoy, B. Chen, V. Vaish, M. Horowitz, I. McDowall, and M. Bolas, “Synthetic Aperture Confocal Imaging,” ACM Transactions on
Graphics (ToG), vol. 23, no. 3, pp. 825–834, 2004.

12[16] M. Levoy and P. Hanrahan, “Light Field Rendering,” in Proc.
SIGGRAPH, 1996.
[17] T. R. Corle and G. S. Kino, Eds., Confocal Scanning Optical Microscopy and Related Imaging Systems. Academic Press, 1996.
[18] C. Fuchs, M. Heinz, M. Levoy, H.-P. Seidel, and H. P. A.
Lensch, “Combining Confocal Imaging and Descattering,” Computer Graphics Forum, Special Issue for the Eurographics Symposium
on Rendering 2008, vol. 27, no. 4, pp. 1245–1253, 6 2008.
[19] S. G. Narasimhan, S. K. Nayar, B. Sun, and S. J. Koppal, “Structured Light in Scattering Media,” in Proc. International Conference
on Computer Vision (ICCV), 2005, pp. 420–427.
[20] J. Gu, S. K. Nayar, E. Grinspun, P. Belhumeur, and R. Ramamoorthi, “Compressive Structured Light for Recovering Inhomogeneous Participating Media,” IEEE Transactions on Pattern Analysis
and Machine Intelligence (TPAMI), vol. 35, no. 3, pp. 555–567, 2013.
[21] B. Lamond, P. Peers, and P. Debevec, “Fast Image-based Separationof Diffuse and Specular Reflections,” in Proc. ACM SIGGRAPH
sketches, 2007.
[22] Y. Mukaigawa, R. Raskar, and Y. Yagi, “Analysis of Light Transport
in Scattering Media,” in Proc. IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2010, pp. 153–160.
[23] ——, “Analysis of Scattering Light Transport in Translucent Media,” IPSJ Transactions on Computer Vision and Applications, vol. 3,
pp. 122–133, 2011.
[24] K. Tanaka, Y. Mukaigawa, Y. Matsushita, and Y. Yagi, “Descattering of Transmissive Observation using Parallel High-Frequency
Illumination,” in Proc. International Conference on Computational
Photography, 2013.
[25] S. Achar, S. T. Nuske, and S. G. Narasimhan, “Compensatingfor Motion during Direct-Global Separation,” in Proc. International
Conference on Computer Vision (ICCV), 2013.
[26] S. Achar and S. G. Narasimhan, “Multi Focus Structured Lightfor Recovering Scene Shape and Global Illumination,” in Proc.
European Conference on Computer Vision (ECCV), 2014.
[27] D. Reddy, R. Ramamoorthi, and B. Curless, “Frequency-space
Decomposition and Acquisition of Light Transport under Spatially
Varying Illumination,” in Proc. European Conference on Computer
Vision (ECCV), 2012, pp. 596–610.
[28] M. Gupta, Y. Tian, S. G. Narasimhan, and L. Zhang, “A Combined
Theory of Defocused Illumination and Global Light Transport,”
International Journal of Computer Vision (IJCV), vol. 98, no. 2, pp.
146–167, Oct. 2011.
[29] M. O’Toole, R. Raskar, and K. N. Kutulakos, “Primal-Dual Coding
to Probe Light Transport,” ACM Transactions on Graphics (ToG),
vol. 31, no. 4, pp. 1–11, jul 2012.
[30] M. O’Toole, J. Mather, and K. N. Kutulakos, “3D Shape andIndirect Appearance by Structured Light Transport,” in Proc. IEEE
Conference on Computer Vision and Pattern Recognition (CVPR).
IEEE, 2014, pp. 3246–3253.
[31] M. O’Toole, S. Achar, S. G. Narasimhan, and K. N. Kutulakos, “Homogeneous Codes for Energy-Efficient Illumination and Imaging,”
ACM Transactions on Graphics (ToG), vol. 34, no. 4, pp. 35:1–35:13,
jul 2015.
[32] F. Heide, L. Xiao, A. Kolb, M. B. Hullin, and W. Heidrich, “Imaging
in Scattering Media using Correlation Image Sensors and Sparse
Convolutional Coding.” Optics express, vol. 22, no. 21, pp. 26 338–
50, oct 2014.
[33] A. Kadambi, R. Whyte, A. Bhandari, L. Streeter, C. Barsi, A. Dorrington, and R. Raskar, “Coded Time of Flight Cameras: Sparse
Deconvolution to Address Multipath Interference and Recover
Time Profiles,” ACM Transactions on Graphics (ToG), vol. 32, no. 6,
pp. 1–10, nov 2013.
[34] M. O’Toole, F. Heide, L. Xiao, M. B. Hullin, W. Heidrich, andK. N. Kutulakos, “Temporal Frequency Probing for 5D Transient
Analysis of Global Light Transport,” ACM Transactions on Graphics(ToG), vol. 33, no. 4, pp. 1–11, jul 2014.
[35] R. Tadano, A. Kumar Pediredla, and A. Veeraraghavan, “Depth Selective Camera: A Direct, On-Chip, Programmable Technique for
Depth Selectivity in Photography,” in Proc. International Conference
on Computer Vision (ICCV), 2015, pp. 3595–3603.
[36] N. J. W. Morris and K. N. Kutulakos, “Reconstructing the Surfaceof Inhomogeneous Transparent Scenes by Scatter Trace Photography,” in Proc. International Conference on Computer Vision (ICCV),
2007.
[37] T. Hawkins, P. Einarsson, and P. Debevec, “Acquisition of Timevarying Participating Media,” ACM Transactions on Graphics (ToG),
vol. 24, no. 3, pp. 812–815, Jul. 2005.

IEEE TRANSACTION ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XX, NO. X, MM YYYY[38] I. Ihrke and M. Magnor, “Adaptive grid optical tomography,”
Graphical Models, vol. 68, no. 5, pp. 484 – 495, 2006, special Issue
on the Vision, Video and Graphics Conference 2005.
[39] A. J. L. Adam, P. C. M. Planken, S. Meloni, and J. Dik, “TeraHertz
imaging of hidden paint layers on canvas,” Optics Express, vol. 17,
no. 5, p. 3407, Feb. 2009.
[40] D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang,
M. Hee, T. Flotte, K. Gregory, C. Puliafito, and A. Et, “Opticalcoherence tomography,” Science, vol. 254, no. 5035, pp. 1178–1181,
Nov. 1991.
[41] A. F. Fercher, “Optical coherence tomography,” Journal of Biomedical Optics, vol. 1, no. 2, pp. 157–173, 1996.
[42] M. G. L. Gustafsson, L. Shao, P. M. Carlton, C. J. R. Wang, I. N.
Golubovskaya, W. Z. Cande, D. A. Agard, and J. W. Sedat, “Threedimensional resolution doubling in wide-field fluorescence microscopy by structured illumination,” Biophysical journal, vol. 94,
no. 12, pp. 4957–70, Jun. 2008.
[43] P. Kner, B. B. Chhun, E. R. Griffis, L. Winoto, and M. G. L.
Gustafsson, “Super-resolution video microscopy of live cells by
structured illumination,” Nature methods, vol. 6, no. 5, pp. 339–42,
May 2009.
[44] O. D. Therrien, B. Aubé, S. Pagès, P. D. Koninck, and D. Côté,“Wide-field multiphoton imaging of cellular dynamics in thick tissue by temporal focusing and patterned illumination,” Biomedical
optics express, vol. 2, no. 3, pp. 696–704, Jan. 2011.
[45] K. Shimizu, K. Tochio, and Y. Kato, “Improvement of transcutaneous fluorescent images with a depth-dependent point-spreadfunction,” Applied Optics, vol. 44, no. 11, pp. 2154–2161, 2005.
[46] S. Narasimhan and S. Nayar, “Shedding Light on the Weather,”in Proc. IEEE Conference on Computer Vision and Pattern Recognition(CVPR), 2003, pp. 665–672.
[47] S. Chandrasekhar, Radiative Transfer. Dover Publications, Inc.,
1960.
[48] A. Ishimaru, Wave Propagation and Scattering in Random Media.
IEEE Press, 1997.
[49] H. W. Jensen, S. R. Marschner, M. Levoy, and P. Hanrahan, “A
Practical Model for Subsurface Light Transport,” in Proc. SIGGRAPH. ACM Press, aug 2001, pp. 511–518.
[50] P. Hanrahan and W. Krueger, “Reflection from Layered Surfaces
due to Subsurface Scattering,” in Proc. SIGGRAPH. ACM Press,
sep 1993, pp. 165–174.
[51] R. Tibshirani, “Regression Shrinkage and Selection Via the Lasso,”
Journal of the Royal Statistical Society, Series B, vol. 58, pp. 267–288,
1994.
[52] P. Pérez, M. Gangnet, and A. Blake, “Poisson Image Editing,” ACM
Transactions on Graphics (ToG), vol. 22, no. 3, pp. 313–318, Jul. 2003.
[53] S. G. Narasimhan, M. Gupta, C. Donner, R. Ramamoorthi, S. K.
Nayar, and H. W. Jensen, “Acquiring Scattering Properties ofParticipating Media by Dilution,” ACM Transactions on Graphics(ToG), vol. 25, pp. 1003–1012, 2006.
[54] N. Kawai and N. Yokoya, “Image Inpainting Considering Symmetric Patterns,” in Proc. International Conference on Pattern Recognition, 2012, pp. 2744–2747.

Kenichiro Tanaka received his B.S. and M.S.
degrees in CS from Osaka University in 2012
and 2014, respectively. He is in pursuit of Ph.D.
degree under supervision of Prof. Yagi, Prof.
Mukaigawa, and Prof. Matsushita at Osaka University. He is a research fellow (DC1) of Japan
Society for the Promotion of Science (JSPS).
His research interests include computer vision,
computational photography, and computational
illumination. He is a member of CVF and a student member of IEEE.

13

Yasuhiro Mukaigawa received his M.E. and
Ph.D. degrees from University of Tsukuba in
1994 and 1997, respectively. He became a research associate at Okayama University in 1997,
an assistant professor at University of Tsukuba
in 2003, an associate professor at Osaka University in 2004, and a professor at Nara Instituteof Science and Technology (NAIST) in 2014. His
current research interests include photometricanalysis and computational photography. He isa member of IEEE.

Hiroyuki Kubo is an assistant professor at Nara
Institute of Science and Technology (NAIST),
Japan since 2014. His research interests include computer graphics and computer animation. He received the MS and Ph.D degrees fromWaseda University, in 2008 and 2012, respectively. He is a member of ACM.

Yasuyuki Matsushita received his B.S., M.S.
and Ph.D. degrees in EECS from the Universityof Tokyo in 1998, 2000, and 2003, respectively.
From April 2003 to March 2015, he was with
Visual Computing group at Microsoft Research
Asia. In April 2015, he joined Osaka University
as a professor. His research area includes computer vision, machine learning and optimization.
He is on the editorial board of IEEE Transactions
on Pattern Analysis and Machine Intelligence(TPAMI), International Journal of Computer Vision (IJCV), IPSJ Journal of Computer Vision and Applications (CVA),The Visual Computer Journal, and Encyclopedia of Computer Vision.
He served/is serving as a Program Co-Chair of PSIVT 2010, 3DIMPVT
2011, ACCV 2012, ICCV 2017, and a General Co-Chair for ACCV 2014.
He is a senior member of IEEE.

Yasushi Yagi is the Executive Vice Presidentof Osaka university. He received his Ph.D. degrees from Osaka University in 1991. In 1985, hejoined the Product Development Laboratory, Mitsubishi Electric Corporation, where he worked
on robotics and inspections. He became a Research Associate in 1990, a Lecturer in 1993,
an Associate Professor in 1996, and a Professor in 2003 at Osaka University. He was also
Director of the Institute of Scientific and Industrial Research, Osaka university from 2012 to2015. International conferences for which he has served as Chair include: FG1998 (Financial Chair), OMINVIS2003 (Organizing chair), ROBIO2006 (Program co-chair), ACCV2007 (Program chair), PSVIT2009(Financial chair), ICRA2009 (Technical Visit Chair), ACCV2009 (Generalchair), ACPR2011 (Program co-chair) and ACPR2013 (General chair).
He has also served as the Editor of IEEE ICRA Conference Editorial
Board (2007–2011). He is the Editorial member of IJCV and the Editorin-Chief of IPSJ Transactions on Computer Vision & Applications. He
was awarded ACM VRST2003 Honorable Mention Award, IEEE ROBIO2006 Finalist of T.J. Tan Best Paper in Robotics, IEEE ICRA2008
Finalist for Best Vision Paper, MIRU2008 Nagao Award, and PSIVT2010
Best Paper Award. His research interests are computer vision, medical
engineering and robotics. He is a fellow of IPSJ and a member of IEICE,
RSJ, and IEEE.