Curve Analogies

Thirteenth Eurographics Workshop on Rendering (2002)
P. Debevec and S. Gibson (Editors)
Curve Analogies
Aaron Hertzmann1
Nuria Oliver2
Brian Curless1
Steven M. Seitz1
1University of Washington
2Microsoft Research
Abstract
This paper describes a method for learning statistical models of 2D curves, and shows how these models can be
used to design line art rendering styles by example. A user can create a new style by providing an example of
the style, e.g. by sketching a curve in a drawing program. Our method can then synthesize random new curves in
this style, and modify existing curves to have the same style as the example. This method can incorporate position
constraints on the resulting curves.
1. Introduction
correspondence in advance, and then iteratively reestimate
curve shape, neighborhood sampling patterns and transla-
Designing new styles is one of the most difﬁcult tasks in
tional/rotational alignments.
non-photorealistic rendering (NPR). Typically, one is given a
choice of a few rendering algorithms and the ability to adjust
This paper has three main contributions: ﬁrst, We show
different parameters to these algorithms. While this gives the
how learning curve styles may be formulated as a texture
user many options, it may be difﬁcult to ﬁnd a desired style
synthesis problem on a parametric domain. Second, we show
in this space. Moreover, it is unlikely that highly sophisti-
how curve analogies (i.e. transformations from input to out-
cated or personal styles will be available.
put curves) may be learned from data, by learning spe-
ciﬁc statistics within curves and relationships between them.
In this paper, we present methods for learning line styles
Third, we show how these algorithms can be used to capture
from examples. With this approach, an artist or end-user may
hand-drawn styles of curves. Although we have focused on
simply draw strokes in the desired style; our system will cap-
2D curve synthesis for non-photorealistic rendering, we ex-
ture aspects of their style and can generate new drawings in
pect that this approach can be extended to other 2D signals,
that style. For example, to design an outline style for a “ner-
and even generalized to 3D surfaces.
vous” character, one may draw a jittery stroke, and to design
a style for a ﬁeld of grass, one may draw an example of its
In this paper, we use the term “curve style” to refer to the
jagged silhouette. In this paper, we adopt the terminology of
shape variations in a curve. Our work addresses one part of a
Hertzmann et al.13, and call the approach “curve analogies.”
larger problem of learning line art styles, including learning
the placement of strokes, as well as the physical appearance
Our work is based on recent algorithms in image texture
of media (i.e. the buildup of graphite or paint). We focus on
synthesis, where a new texture is generated such that every
the speciﬁc subproblem of curve statistics, and assume that
pixel neighborhood in the new texture looks like some neigh-
the input and output curves are continuous.
borhood in the example texture. More recently, several meth-
ods have been developed that try to learn transformations be-
tween images, based on local neighborhood transformations.
2. Related Work
Our goal is to automatically learn how to generate an output
There are many different ways to design curve styles
curve from an input curve. This problem is surprisingly dif-
for illustration. One general strategy is to manually pro-
ﬁcult, because neighborhood sampling patterns depend on
gram curve shape procedures, an approach taken by sev-
the shape of the output curve. This means we must some-
eral authors6, 21, 23, 27. This approach gives great control to
how generate the output curve shape, its neighborhood pa-
the programmer, but can be somewhat unintuitive. An al-
rameterizations, and its correspondence to the source curve
ternative approach that we adopt is to create curve styles
simultaneously, even though each one depends on the oth-
from examples. This approach has the potential to provide
ers. Moreover, we must also account for possible rotations
an intuitive interface — since the user is required only to
and translations in drawings. Our approach is to ﬁx a curve
draw an example of the desired style — but requires an ef-
c The Eurographics Association 2002.

fective procedure for capturing style. One approach is to
p
p
2
p
7
simply warp an example stroke to the desired shape, as
p
p
1
3
f(t)
p
8
done by Hsu et al.16, or to copy displacements to the tar-
p
6
0
p
p
9
get shape, as done by Finkelstein and Salesin9. These ap-
4
p5
proaches give high quality results with relatively little effort,
but may give a distorted result when the target shape is very
different from the original shape, or when the base shapes
are noisy. Such noise can be ﬁltered out, at the cost of de-
t t t
t t
t
t
t
t t
tail in the source shape. Essentially, such approaches require
0
1 2
3 4
5
6
7
8
9
that all texture information resides in the high-frequencies
Figure 1: We represent each curve as a polyline: a list of
of the texture curve, and all “sweep” information is in the
control points p
low-frequencies of the target shape. These methods also pro-
i connected by straight line segments. Each
control point has a corresponding parameter value t
duce noticeable repetition, unless many example curves are
i and
speciﬁes that f (t
provided. Markosian et al.22 describe a system for placing
i) = pi. The curve is evaluated at a parame-
ter value t by linear interpolation between the adjacent con-
copies of small, hand-made strokes on a surface; in contrast,
trol points.
we focus on creating long, continuous curves. Computer vi-
sion and handwriting recognition researchers have put sub-
stantial effort into shape recognition2, 3, 29; to our knowledge,
these methods have not been applied to modifying shapes
or to line art synthesis. Some service bureaus will create a
Freeman and Hertzmann et al., such methods can be used
personalized typeface from a handwriting sample; our work
to produce line drawings by example. However, they can be
generalizes this approach. Freeman et al.11 present perhaps
quite slow, since they must process large amounts of train-
the ﬁrst system for synthesizing line art based on example
ing data (i.e. the set of image neighborhoods), often blur and
styles, and provide an inspiration for our work. This method
blend curves, and produce resolution-dependent bitmaps (as
produces high-quality results, but requires a large training
opposed to line art in vector form). Furthermore, they do not
corpus (e.g. over a hundred strokes), only handles strokes
capture rotationally invariant aspects of style. We show that
of roughly the same lengths as the example strokes, and
these methods can be transferred to the curve domain, and
does not incorporate other constraints (e.g. the method may
that this representation is better suited to learning line styles.
change the topology of the drawing in undesirable ways).
However, it should be pointed out that our method cannot yet
Our work addresses all of these problems. In work concur-
work directly from source images but, rather, requires input
rent to our own, Kalnins et al.18 describe a user interface for
curves in vector form.
designing rendering styles for 3D models. Our work is very
much in this spirit. Their method also includes a curve tex-
Wei and Levoy26 and Ying et al.28 describe methods for
ture synthesis procedure, but may suffer the same problems
synthesizing 3D shape by synthesizing displacement tex-
as copying displacements (as described above). Our method
tures. These methods model statistics of only the displace-
could be used as a more-powerful curve synthesis procedure
ments, whereas our method models statistics of the resulting
within their framework.
surface shape. Furthermore, we show how to learn transfor-
mations between shapes.
A few authors have described ways to learn stroke ar-
rangements. Salisbury et al.23 create hatchings interactively
At a high level, our constrained texture synthesis algo-
by copying from example hatchings. Chen et al.5 describe
rithms are similar to variational curve modeling. For exam-
a method for learning facial illustration styles from exam-
ple, Kobbelt and Schröder19 describe a subdivision scheme
ple. Their method does not produce very stylized lines; our
for fair interpolation of a set of constraints. Instead of using a
method could be used together with theirs to produce styl-
predetermined smoothness energy function, we use a learned
ized line drawings directly from images. Jodoin et al.17 de-
energy function; our single-resolution synthesis is analogous
scribe a system for learning hatchings as arrangements of
to a blurring ﬁlter, and our multiresolution scheme is analo-
strokes; our work is complementary, since we focus on the
gous to variational subdivison.
appearance of individual strokes.
This paper builds on recent work in example-based im-
3. Algorithms
age texture synthesis. Freeman et al.10 describe an algorithm
for learning scene interpretations from examples, and, sub-
In this section, we introduce a family of curve analogies
sequently, Efros and Leung8 demonstrated that simple recur-
problems and corresponding algorithms. The most general
sive sampling can produce new textures similar to examples.
problem statement is as follows: given an example “unﬁl-
Wei and Levoy25 demonstrated several accelerations to this
tered” curve A and example “ﬁltered” curve A , we would
technique, including multiresolution synthesis. Efros and
like to “learn” the transformation from A to A , and apply it
Freeman7, Harrison12 and Hertzmann et al.13 demonstrated
to a new curve B to generate an output curve B . In short, we
ways to generalize these synthesis algorithms by constrain-
wish to generate B to satisfy the relation
ing them with source images. As demonstrated by Efros and
A : A :: B : B
(1)
2

A'
A'(s )
j
A'
a' +t'
k
b'
B'(t )
t'
k
B'
i
B'
(a)
(b)
(c)
Figure 2: Curve synthesis problem statement. (a) Our goal is to produce a curve B such that the neighborhood around every
point B (ti) “looks like” some neighborhood in A . (b) Each neighborhood can be sampled with arc length spacing to get
samples ak and bk. In order to compare neighborhoods, we estimate a rigid transformation that aligns the samples. Here we
show a translation t that maps from the coordinate system of the A neighborhood to that of the B neighborhood. (c) Once
the transformation is determined, the two neighborhoods can be compared in the same coordinate frame. Each neighborhood
is sampled with arc length spacing, and the corresponding samples ak + t and bk are compared.
We begin by describing the simplest problem statement and
algorithm8 for processing shape. We use the primed symbols
algorithm, and then show how to add more features to create
A and B in anticipation of Section 3.3, when the unprimed
more general algorithms. Speciﬁcally, we begin with a curve
symbols will be introduced in a more general version of the
synthesis algorithm that randomly generates a B curve in
algorithm.
the style of an A curve. In this case, there are no A or B
In particular, we use A to deﬁne a cost function over pos-
curves and we are just randomly sampling curves from the
sible curves, and then attempt to generate a new curve B
distribution implied by A . Next, we show how to synthesize
a shape according to some constraints. We then show how to
with minimal cost:‡
generalize this method to create curve analogies, and, ﬁnally,
E(B ) = ∑mind(B ,ti,A ,sj)
(2)
how multiresolution synthesis can give greater speed, quality
j
i
and ﬂexibility.
This cost function states that we desire curves for which the
A 2D curve can be written in parametric form as a map-
neighborhood around each ti location in B “looks like” the
ping from a segment on the real line to the plane f :
neighborhood around some s j in A (Figure 2(a)). In other
[t
words, this cost function measures the “difference” between
min, tmax] → R2. For convenience, we use piecewise linear
curve representations, a.k.a. polylines (Figure 1). Each curve
the local shape of B around ti and the local shape of A
is represented internally by an ordered list of control points
around s j. We restrict the ti and s j samples to be taken from
(t
a ﬁnite set of samples in their respective domains.
i, pi), where each parameter value ti is unique in the list,
and the maximum and minimum parameter values in the list
Of critical importance is the deﬁnition of the neigh-
deﬁne the range [tmin,tmax]. A curve is evaluated at a speciﬁc
borhood distance metric d(B ,ti, A , s j). We represent each
t value by linear interpolation†, and rendered by drawing line
neighborhood as a set of K samples ak ≡ A (sk) and
segments between each control point. The curve is modiﬁed
bk = B (tk), k = {1..K}, taken in unit arc length incre-
by adding and deleting control points. A curve computed in
ments around A (s j) and B (ti), respectively. In addition, we
this representation can easily be converted to other repre-
use tangent features ∆b
sentations, e.g. by resampling or by creating a non-uniform
k ≡ (bk − bk−1)/ bk − bk−1 , and
∆a
B-spline from the control points.
k ≡ (ak − ak−1)/ ak − ak−1 , in order to better capture
the smoothness properties of the curve. We cannot directly
compare the 2D positions of points on the curves B and
3.1. Curve synthesis
A (as is typically done for image texture synthesis), since
translating a curve in 2D should not affect this measurement.
Problem statement. The curve texture synthesis problem is
More generally, we would like the neighborhood compari-
as follows: given an example curve A , generate a new curve
son to be invariant to rigid transformations. Hence, we use
B of a speciﬁed arc length LB that appears “similar” to the
the distance metric
example curve (Figure 2(a)).
d(B ,t
2
i, A , s j ) = min ∑ wk( R ak + t − bk
+
Our single-scale curve texture synthesis algorithm is an
R ,t
k
adaptation of Efros and Leung’s image texture synthesis
w
2
∆ R ∆ak − ∆bk
)
p
if (t, p) is a control point
‡
†
More formally, we infer a density of curves p(B ) from A , and
f (t) =
(t−tlower)pupper+(tupper−t)plower
otherwise
then sample from this density. The density is given by p(B ) =
tupper−tlower
where “upper” and “lower” denote the parameter values of the
e−E(B )/Z, where Z is a normalization constant. In this sense, our
control points immediately before and after t.
algorithm models the statistics of curves as samples from p(B ).
3

A'(S(i))
B'(ti)
A'
B'
(a)
(b)
A'(S(i-1))
B'(ti-1)
Figure 3: Curve coherence. (a) We would like to copy continuous segments of the example curve A to B . We say that a curve
segment is coherent if it was copied from a curve segment in the example curve. (b) For the polyline representation, we can test
for coherence by testing if the line segment between A (S(i)) and A (S(i − 1)) is identical to the line segment between B (ti) and
B (ti−1). We approximate this test by measuring the arc lengths A (S(i)) − A (S(i − 1)) and B (ti) − B (ti−1) .
A'
A'(s )
j
A'
a' +t'
k
p =A'(s )+t'
j
j
b'
B'(t
)
k
B'
max
B'(t
)
t'
new
B'
(a)
(b)
(c)
Figure 4: First pass of the single-scale curve synthesis algorithm. (a) Our goal is to pick a value of B (tnew) such that the
resulting neighborhood around tnew matches some neighborhood in A . (b) For each example neighborhood, we ﬁnd the rigid
translation that best aligns the two neighborhoods. Here we show a translation t . (c) The neighborhoods are compared in the
transformed coordinate frame, and the candidate p j for B (tnew) is a copy of A (s j) in the local coordinate frame.
where the rotation R and translation t together deﬁne a rigid
ment of three consecutive control points from the A curve
transformation, and wk and w∆ are user-deﬁned parameter
to the B curve.
settings. If only translation invariance is desired, then R can
The remainder of the B curve is generated by nearest-
be ﬁxed to be the identity matrix; this choice of whether to
neighbor sampling. We generate a new sample after the end
use rotation invariance depends on the texture being mod-
of the curve, at t
eled.
new = tmax + ∆t, where tmax is the current
maximum t value in B , and ∆t is a predetermined spacing
Several authors have observed that, for image textures,
(Figure 4(a)). Our goal now is to choose the value B (tnew)
copying patches of texture gives improved quality1, 7, 13, 20;
that minimizes equation (2), while holding the rest of the
we ﬁnd the same to be true for curves. Hence, in a similar
curve shape ﬁxed. As in the texture synthesis algorithms, the
manner to Hertzmann et al.13, we encourage copying coher-
algorithm searches for the value p∗ for B (tnew) that best
ent segments from A to B . Recall that B is constructed
matches the resulting neighborhood of some neighborhood
from control points (ti, pi) so that B (ti) = pi. Let S(i) be
j∗ in A , while neglecting the other terms in Equation 2 that
the source index for each control point in B , deﬁned as
depend on B (tnew):
S(i) = arg min j d(B ,ti, A , s j). We say that a control point
(
(
t
p∗, j∗) = arg
min
d(B ,tnew, A , s j)
(3)
i, pi) in B
is coherent with the preceding control point
(B (tnew), j)
(ti−1, pi−1) if S(i−1) < S(i) and the curve segment between
A (S(i − 1)) and A (S(i)) is identical to the curve segment
In order to estimate these values, the algorithm iterates over
between B (t
each neighborhood (indexed by j) in the example curve
i−1) and B (ti) (Figure 3). We can test for this
condition approximately by testing whether the correspond-
A . For each j, we compute the new value p j for B (tnew)
ing arc lengths are the same, which, in this case, reduces to
that minimizes d(B ,tnew, A , s j), and the corresponding cost
A (S(i)) − A (S(i − 1)) = B (ti) − B (ti−1) . We penal-
d j = arg minB (tnew) d(B ,tnew,A ,s j). We set j∗ to the index
ize non-coherent control points by multiplying the distance
of the neighborhood with the smallest cost d j, and p∗ to the
d(·) by (1 + κ/D) where D is the approximate control point
corresponding position p j. Finally, we set B (tnew) to the
spacing in the curve and κ is a parameter setting.
value p∗ corresponding to that sample, and the source in-
dex S(i) = j∗. We introduce randomness in the choice of the
source location j∗ in the same way as Efros and Leung8,
i.e. we uniformly sample from the set of s values for which
Algorithm. We now describe a single-scale curve synthesis
the neighborhood measurement is within a ﬁxed proportion
procedure; pseudocode for the algorithm is given at the end
(1 + ε) of the optimal choice.
of this section. Since, at the outset, the start of the curve is
entirely unconstrained, the algorithm ﬁrst initializes B (t) to
We now describe the computation of p j and d j for each
be an empty curve, and then copies a randomly-chosen seg-
neighborhood j. This search is very similar to the one de-
4

scribed by Efros and Leung8, except that neighborhoods
Once an initial shape for this curve has been generated, it
must be compared under rigid transformations. First, the
can be reﬁned by making successive passes over the curve.
method computes the optimal rigid transformation R j, t j
The algorithm for this is very similar to the ﬁrst pass: it iter-
that aligns the existing neighborhood samples {ak, ∆ak}
ates over every control point in B , collects the full neighbor-
to {bk, ∆bk} (Figure 4(b)), using standard point-matching
hood around that control point, and ﬁnds the best matching
techniques14, 15, 24. This 2D transformation is computed in
full neighborhood in the example data.
closed form. A causal neighborhood is used: the neighbor-
hood points in bk after ti are omitted (since they have not yet
3.2. Curve synthesis with constraints
been computed), along with their corresponding points in ak.
Note that if only translation invariance is desired, then R j is
Problem statement. We now describe a way to add con-
set to be an identity matrix. With these values for j, R j, and
straints to curve synthesis. We allow two kinds of con-
t
straints:
j , the cost function can be written
• Soft constraints specify that the curve should pass near
E(B ) = wnew p j − (R jA (s j) + t j) 2 +C
(4)
speciﬁc positions. Each soft constraint adds a term of the
form w
2
c B (tc) − qc
to the cost function in Equation (2),
where wnew is the weight of B (tnew) in equation (2), and C
where c is the index of the constraint.
contains terms that do not depend on p j. Hence, the optimal
• Hard constraints specify that the resulting curve must pass
choice for p j is given by transforming the position of A (s j)
through speciﬁc positions qc. A hard constraint is speciﬁed
to the neighborhood in B : this gives p j = R jA (s j) + t j
as B (tc) = qc and can be viewed as the limit of a soft con-
(Figure 4(c)). Finally, we evaluate the cost d j for this value
straint as wc → ∞.
p j, and apply the coherence penalty, if any. We determine
We represent a constraint as a tuple
t
whether a candidate is coherent by the approximate test
c, qc, wc , where
w
S(i − 1) < s
c = ∞ for hard constraints.
j and
A (S(i)) − A (S(i − 1)) < 3 B (t
2
i) −
B (ti−1) .
Note that a pair of hard constraints
t1, p1, ∞ and
t
This process can be accelerated by precomputing the set
2, p2, ∞
speciﬁes only that the curve must pass through
p
of neighborhoods in A . This sped up our system by a factor
1 and p2 in that order; the arc length of the curve between
those two points may vary widely depending on the style of
of ﬁve.
the training data. In practice, the arc length will also depend
This method is not strictly optimal because the sampling
on the search procedures used.
pattern and the rigid transformation all depend on the po-
sition of p j, and vice versa. For curve synthesis, we have
Algorithm. The algorithm for curve synthesis with con-
not found it necessary to directly account for this interde-
straints is quite similar to that without. In addition to the
pendence — we ﬁrst estimate the sampling pattern, then the
example curve A , a set of constraints tc, pc, wc is also pro-
rigid transformation, then the value of p j. In later sections,
vided. For this version of the algorithm, we assume that
we will describe cases where multiple passes are necessary.
the ﬁrst point on B is constrained; we use the constraint
to initialize the curve. The more general case can be han-
The entire algorithm can be summarized in pseudocode as
dled by synthesizing forward from the ﬁrst constraint, and
follows:
then synthesizing backwards from the last constraint. The
single-scale synthesis is not robust to arbitrary arrangements
function SYNTHESIZECURVE(A ):
of constraints (for example, it has little chance of reaching
initialize B with a randomly-chosen sequence of
a constraint that is far from the previous constraint), so we
three control points from A
use the multiscale synthesis algorithm (to be described in
while ARCLENGTH(B ) < LB
Section 3.4) whenever constraints are present. We only in-
tnew ← tmax + ∆t
troduce synthesis with constraints in isolation to make the
compute the curve samples bk around B (tmax)
presentation cleaner.
for each s j in a discrete set in the domain of A , do:
compute the curve samples a
For each pass of the algorithm, we iterate over the same t
k around A (s j )
compute R
values as before, but also insert control points (and perform
j , t j that align {ak, ∆ak} to {bk, ∆bk}
p
synthesis) at the constraints. Unconstrained control points
j ← R jA (s j) + t j
are synthesized exactly as in the previous algorithm. Con-
compute d j for the values of p j, R j, and t j
trol points speciﬁed by hard constraints are immediately re-
if p j is not coherent, d j ← d j(1 + κ/D)
placed with the position of the hard constraint.
dmin ← min j d j
compute the set of candidates { j : d j(1 + ε) ≤ dmin}
For control points with soft constraints, we must modify
choose j∗ by uniform sampling from the candidates
the search procedure according to the modiﬁed cost function.
insert control point (tnew, p j∗ ) into B
In the case of one soft constraint tc, qc, wc when tc = tnew,
S(i) ← s j∗
equation (4) becomes
return B
E(B ) = w
2
new p j − (R jA (s j) + t j) 2 + wc p j − qc
+C
5

R',t'
B
B'
B
B'
A'(s )
j
A
A
A'
B(t )
B'(t )
A'
A(s )
i
i
j
R',t
(a)
(b)
Figure 5: Curve analogies problem statement. Here we show an example using rotation invariance, where the B curve is rotated
90 degrees from the A curve. (a) Our goal is to generate a B curve, such that, for every point ti, there is some point s j such
that: the neighborhood around B (ti) looks like the neighborhood around A (s j), the neighborhood around B(ti) looks like the
neighborhood around A(s j), and the offset between B and B matches the offset between A and A . The offsets are measured
between the neighborhood centers-of-mass. (c) These comparisons are performed with respect to the rigid transformation that
best matches the curves. We use the same rotation matrix R to map between the different curves, but allow different translations
t (between A and B) and t (between A and B ). A different rigid transformation will be estimated for each candidate.
where wnew is the weight of the target point; we always use
B to satisfy the relation
wnew = 1 in practice. The overall outline of the algorithm is
the same: we compute the set of candidates p j and their asso-
A : A :: B : B
ciated costs d j, and choose the candidate with the least cost.
However, the candidates are now computed by optimizing
Of course, there are inﬁnitely many plausible ways to inter-
the above quadratic equation:
pret a training pair; in this paper, we pursue a “textural” form
w
of analogy, similar to the image methods7, 12, 13. However, the
new(R jA (s j) + t j) + wcqc
p j =
(5)
shape and position relationships of 2D curves present some
wnew + wc
additional complications.
The position of p j is no longer a direct copy of the example
First, we need some way of representing the correspon-
position as in Efros and Leung’s method — the constraint
dence between A and A , and between B and B . We do this
acts like a spring that “tugs” the value of B (tnew) toward q j.
by requiring the A, A curves to have the same parameter-
For a speciﬁc neighborhood and transformation, the cost
ization, and generating B with the same parameterization
in equation (4) acts like a soft constraint of the form
as B. Hence, a point A(s j) corresponds to the point A (s j),
tnew, R jA (s j) + t j, wnew . It is convenient to view this as
and B(ti) corresponds to B (ti). Note that this places no con-
a “virtual constraint;” in subsequent sections, we will add
straint on the relative arc lengths of the curves. If the exam-
more terms that can be written as virtual constraints. In
ple curves do not have the same parameterization, then they
our implementation, we convert all energy terms into vir-
are ﬁrst reparameterized.
tual constraints, and then compute p j as a weighted com-
bination p
More importantly, we need some way to measure shape
j = ∑c wcqc/ ∑c wc, for each constraint c, real or
relationships between input and output curves. For example,
virtual. The cost terms ∑
2
c wc p j − qc
are then added to d j.
one approach would be to represent the A and B curves as
This uniﬁed treatment of constraints makes it easier for one
displacement vectors from the A and B curves respectively.
piece of code to handle the different possible cost functions.
However, this representation would be very sensitive to noise
Of course, our technique could be generalized to allow non-
in the shape of the A and B curves — noise in the B curve
quadratic cost, for example, by computing p j by numerical
would distort the resulting B .
optimization.
Intuitively, the problem can be viewed as simultaneously
optimizing for three goals (Figure 5(a)):
3.3. Curve analogies
• The shape of B should “look like” the shape of A ,
Problem statement. The curve analogies problem is as fol-
• The shape of B should have the same relationship to the
lows: given an example “unﬁltered” curve A and example
shape of B as the shape of A has to A, and
“ﬁltered” curve A , we would like to “learn” the transforma-
• The relative positions and orientations of the B curve
tion from A to A , and apply it to a new curve B to generate an
should have the same relationship to B as the positions and
output curve B (Figure 5(a)). In short, we wish to generate
orientations of A do to A.
6

The ﬁrst two goals directly correspond to the matching terms
generating the value of B (tnew), if the neighborhood of
in the Image Analogies algorithms; we optimize these by
B(tnew) has a sharp concave curve, then we prefer to choose
using the same point and tangent features for B and B as in
an example where A(s j) also has a sharp curve. The optimal
Section 3.1. However, the third goal is new, and speciﬁc to
translations t and t are computed for each pair of curves
shape modeling.
separately, and then a single R is computed that aligns the
pairs of curves.
We measure the positional relationships by measuring the
offsets between curves. We use neighborhood centers-of-
Second, the curve offset term affects both the matching
mass as a sort of “summary” of the local neighborhood posi-
penalty d j and the computation of the optimal p j. Since p j
tion; we experimented with several alternatives (such as us-
is a candidate for B (tnew), we can rewrite the offset term as
ing offsets between corresponding curve points) and found
centers-of-mass to be the most effective. Let ak, a
d
k, bk, bk
OFFSET (B, B , ti, A, A , s j ) =
be neighborhood samples in the four curves around some t
∑
i
m2 p
k wk bk
2
j − 1
b + R (a − a) −
and s
m
w j+∑k wk
j sample points, and deﬁne the centers-of-mass as a =
∑wkak/∑wk,a = ∑wkak/∑wk,b = ∑wkbk/∑wk,b =
where we have deﬁned m = w j/(w j + ∑k wk) is the propor-
∑wkb
∑
k/ ∑ wk, and let R and t be the rigid transformation that
tion of mass of b due to p
k wk bk
j , so that b = mp j +
.
aligns the neighborhood of A to B. As before, these samples
w j+∑k wk
This is another quadratic error term, and can be converted
are taken in unit arc length increments around the neighbor-
into a virtual constraint
hood centers; the bk values are resampled from the portion
of the B curve that has already been generated. We want the
1
∑k wkbk
offset from A to A to match the offset from B to B (Figure
ti,
b + R (a − a) −
, wOFFSET m2
m
w j + ∑k wk
5(b)), which we express as
and included when solving for p j. This also means that we
dOFFSET (B, B ,ti, A, A , s j) = R (a − a) − (b − b) 2
do not begin the curve with random initialization. For the
ﬁrst sample in B , there are no b values, and the virtual con-
For example, when R is the identity matrix and A is above
straint reduces to t
A, B should be above B.
i, b + R (a − a), wOFFSET .
The analogy algorithm terminates when the maximum t
The new cost function for the output curve becomes
value in B is the same as in A . Multiple passes may be
E(B ) = ∑min d(B ,ti,A ,sj) + wBd(B,ti,A,sj)+
added to the curve analogy procedure as before, by consid-
j
i
ering full neighborhoods in B , and iterating over all control
w
points in B .
OFFSET dOFFSET (B, B , ti, A, A , s j )
(6)
The single-scale analogies algorithm can be summarized
in pseudocode as follows:
where the three terms here correspond to the three goals de-
scribed above, with weights wB and wOFFSET . The sum of
the three desires is also multiplied by a coherence term. The
function SYNTHESIZEANALOGY(A, A , B, constraints):
weight w
initialize B as an empty curve
B controls the relative importance of the B and B
curves in the matching.
for each constraint {tc, qc, wc}, do:
insert (tc, B(tc)) as a control point into B
As stated, this cost function allows different rotations and
for each control point (ti, pi) ∈ B, do:
translations for the input and output curves. We have found it
compute samples bk and bk around B(ti) and B (ti)
works best to use separate translations t and t for matching
for each s j in a discrete set in the domain of A , do:
A to B and A to B , respectively, but to use the same rotation
compute samples ak and ak around A(s j), A (s j)
from A to B as from A to B .
compute the R j, t j, t j that align {ak, ∆ak, ak, ∆ak}
Constraints may be included in the analogy just as in
to {bk, ∆bk, bk, ∆bk}
curve synthesis; note that the constraint t values correspond
get all constraints { ti, wc, qc } for ti
to the t values on the B curve.
p j ← ∑c wcqc/ ∑c wc
d
2
j ← ∑
In many cases, we are not interested in the shape of the
c wc p j − qc
if p j is not coherent, d j ← d j(1 + κ/D)
B curve, only in its overall sweep. In this case, we set wB to
dmin ← min j d j
be zero, which corresponds to disabling the “second goal”
compute the set of candidates { j : d j(1 + ε) ≤ dmin}
described above.
choose j∗ by uniform sampling from the candidates
insert control point (ti, p j∗ ) into B
Algorithm. The analogy algorithm is the same as curve
S(i) ← s j∗
synthesis, but with two important changes to the cost func-
return B
tion in equation (2), due to the two new cost terms.
First, the A/B curve matching term (6) affects the curve
In
addition,
we
deﬁne
a
function
alignment and the computation of d j. For example, when
REFINEANALOGY(A, A , B, ˆ
B ) that performs the same
7

task as the above procedure, but using ˆ
B as the initialization
Algorithm. These observations lead to the following mul-
of B . This means that full, non-causal neighborhoods are
tiresolution curve synthesis algorithm. First, compute the
computed instead of causual neighborhoods. Multipass
coarsest level of the B pyramid by curve synthesis from the
reﬁnement can be performed by ﬁrst synthesizing with
coarsest A . In order to remove the dependence on the ∆t
SYNTHESIZEANALOGY and then repeatedly applying
for the remaining levels, we resample the curve at unit arc
REFINEANALOGY.
length intervals after synthesis. We then compute the remain-
ing levels of the pyramid using the analogy algorithm from
Section 3.3. Each ﬁner level of the pyramid is initialized with
3.4. Multiresolution curve synthesis
the coarser level, and resampled to have a control point den-
There are a number of problems with the single-scale al-
sity proportional to the neighborhood sampling density be-
gorithms just described. First, we would like to be able to
fore applying the analogy reﬁnement. If any constraints are
use large neighborhood sizes, which can be computationally
present, the corresponding blurred versions are used at each
expensive. Second, we would like to propagate constraints
level.
from the end of the output curve to the beginning, which
The multiresolution synthesis algorithm can be summa-
may require several passes over the curve. Finally, there is
rized in pseudocode as follows:
a dependence on the ∆t step size parameter that can pro-
duce poor results when processing analogies, since the step
function SYNTHESIZECURVEMULTIRES(A ):
size is not directly related to arc length. We address these
create Gaussian pyramid of A
problems with a multiresolution synthesis procedure, based
ˆ
B
on the work of Wei and Levoy25. The idea is to generate
0 ← SYNTHESIZECURVE(A0)
B
a Gaussian pyramid of curves with the same statistics as a
0 ← RESAMPLE( ˆ
B0, 1)
for each remaining level , from coarest to ﬁnest, do:
Gaussian pyramid of curves corresponding to A . The Gaus-
ˆ
sian pyramid of curves is a list of L curves A such that A
B ← RESAMPLE(B
−1
−1, .5(1.2) −L)
is a blurred and subsampled version of A , and A
B ← REFINEANALOGY(A −1, A , B −1, ˆ
B )
L = A . This
deﬁnition of curve pyramid is analogous to Gaussian pyra-
return BL
mids used for images4.
The function RESAMPLE resamples a curve to have the
Problem statement. In multiresolution curve synthesis, we
speciﬁed arc length between control points.
are given an example A curve, and synthesize an output
curve B . This is done by synthesizing a Gaussian pyramid
of curves from coarse to ﬁne, by sampling from the distribu-
3.5. Multiresolution curve analogies
tion implied by the Gaussian pyramid of the A curve. The
ﬁnal output B is the ﬁnest level of the output pyramid.
The most general version of our algorithm is a direct exten-
sion of the preceding algorithms. In this case, synthesis at the
The coarsest level of the output pyramid is a single-scale
ﬁner levels of the pyramid can be viewed as a “generalized
curve synthesis problem based on the coarsest level of the
analogy:”
input pyramid. For the remaining levels of the pyramid, the
relationships can be expressed as an analogy:
A , A −1 : A :: B , B −1 : B
A −1 : A :: B −1 : B
The new curves simply add the corresponding terms to the
where indicates the pyramid level, − 1 is the next coarser
cost function — each curve pair has a matching term, and
level above .
there are two offset terms, to ensure that B has consistent
relationships with B and B −1. Each of these terms creates
Any constraints placed on the output curve must be ap-
corresponding terms in the neighborhood matching and the
plied to the higher levels as well. However, because the
virtual constraints; the resulting algorithm is a trivial gener-
coarser levels are smoothed versions of the ﬁner levels,
alization of the previous algorithms. This involves generaliz-
maintaining the constraints precisely may distort the texture
ing REFINEANALOGY to take multiple pairs of A, B curves
in undesirable ways. For example, if the coarsest level of A
that samples neighborhoods for alignment and comparison
is a straight line, then this texture cannot meet constraints
in each corresponding A, B pair.
that the ﬁner level can. The texture would have to bend to
meet non-colinear constraints, which would not be a valid
The multiresolution analogies algorithm can be summa-
sample of the texture. We address this problem by softening
rized in pseudocode as follows:
the constraints: the weight wc of each constraint tc, qc, wc
is replaced with (w−1
c
+ (L − )−1)−1 at level , where we
deﬁne ∞−1 = 0 for hard constraints.§
σ2c = w−1
c
. Softening corresponds to convolving the prior with an-
other Gaussian of variance (L − )−1; the resulting Gaussian density
§ This softening comes from viewing a constraint as a Gaussian
has variance (w−1
c
+ (L − )−1)−1. Hard constraints have zero vari-
prior probability density over the position of B (tc), with variance
ance.
8

in an analogy to apply this texture to the sweep of another
function SYNTHESIZEANALOGYMULTIRES(A, A , B):
curve. Since we are only interested in the sweep of the curve,
create Gaussian pyramids of A, A and B
ˆ
we set w
B
B to zero.
0 ← SYNTHESIZEANALOGY(A0, A0, B0)
B0 ← RESAMPLE( ˆ
B0, 1)
Figure 9 demonstrates a situation where the curve of A
for each remaining level , from coarest to ﬁnest, do:
is used, and thus wB is set to 2. The algorithm generates a
ˆ
B ← RESAMPLE(B
hand-drawn style based on the example, with different fea-
−1, .5(1.2) −L)
B ← REFINEANALOGY({A , A
tures at corners, vertical edges, and so on. Translation invari-
−1}, A ,
{B , B
ance is also used to capture orientation-dependent behavior.
−1}, ˆ
B )
return B
This means that, for example, the portions of B that were
L
drawn right-to-left are treated differently from those that
were drawn left-to-right. Both A and B were drawn clock-
wise. This explains the indentations that appear below the
3.6. Parameter settings
horizontal part of B : there are no examples of long horizon-
In practice, the parameters must be chosen carefully to yield
tal, right-to-left lines in A, so the long horizontal, right-to-
good results. However, we found that a single set of param-
left lines are replaced with a series of corners.
eters is effective for a large class of styles — we were able
Figure 10 illustrates the use of constraints to compose
to use almost the same settings for all of the results shown
drawings. Processing each edge separately produces curves
in this paper, as described below. All distance measurements
that do not intersect at the endpoints, and thus changes the
are in units of image-space distance.
topology of the drawing (Figure 10(b)). We created endpoint
When creating Gaussian pyramids, we resample each
constraints for each set of nearby endpoints; processing with
level to have a spacing of D = .5(1.2) −L between sam-
these constraints forces the corners to meet and preserves the
ples. We set the neighborhood sampling spacing (i.e. the
source topology.
distance between adjacent ak or bk samples) to be D/4
Figure 11 shows a larger composition, made by applying
for the level being synthesized. Note that the neighbor-
various curve style to elements of a line-drawn composition.
hood size should be large enough to capture the details
Similar methods could be used to for rendering silhouettes
of the style in A . We typically use neighborhoods of
extracted from 3D models in hand-drawn styles.
n = 81 or n = 101 samples wide. The sample weights wk
are approximations to a nonnormalized Gaussian density:
Figure 12 shows another example where the shape of the
n
n
source curve is used, and thus wB is set to 2. Here, we create
wk =
/
. Since this weighting de-
k
(n − 1)/2
a ﬁlter that makes the shape less rigid and more ﬂowing. For
pends only on k, the Gaussian has a larger spread at coarser
example, the bend at the beginning of the “M” is copied to
pyramid levels. The weight of the coarse scale w
the beginning of the “B.”
B −1 = 2.
Typical parameter values include wOFFSET = 3.8 for analo-
gies and wOFFSET = 19 for curve synthesis; w∆ = 100 or
5. Discussion and Future Work
w∆ = 100, 000; κ = 1/2, ∆t = 2D or ∆t = 4D.
We have presented a technique for learning transformations
The weight of source B curve is set to wB = 2 or wB = 0
of 2D shapes from examples. This method can be applied
depending on whether the shape of the B curve is important.
to a wide variety of input shapes, and thus we expect this
If wB = 0, then only the position and orientation information
method to be useful whenever hand-made stroke styles are
of B is used.
desired, such as in hatching and 3D silhouette rendering.
There are many potential avenues for improvement and en-
4. Applications and Experiments
hancement:
We now show several applications of the curve analogies
Learning parameters. Our current system includes an un-
framework. We used a simple 2D drawing interface to create
fortunate number of user-tuned parameters. It would be more
the curves and specify styles and constraints, although our
desirable to ﬁt a parametric model that automatically learns
methods can be applied to any 2D curves (such as object sil-
all parameters from the data. However, it is not at all clear
houettes rendered from a 3D model). Each of the individual
what model would be appropriate for representing arbi-
curves here took from a few seconds to about 1 minute to
trary styles; in particular, traditional time-series models of-
generate using a 1 GHz Pentium 4, depending on the length
ten have difﬁculty representing sharp features without using
of the B curve. This time could probably be substantially
a very large number of parameters (and thus requiring a vast
reduced by hand-optimizing the code. We used multireso-
corpus of training data).
lution synthesis in all experiments. Rotational invariance is
used for all ﬁgures except where noted.
Alternative representations. The polyline representation
Figure 6 illustrates curve synthesis using translation in-
does not represent smooth curves very efﬁciently; it should
variance — a long, loopy curve is automatically generated
be straightforward to translate our methods to a smooth
in the style of a short one. Figure 7 shows this texture used
curve representation (e.g. B-splines). Furthermore, such a
9

representation would be much faster, since fewer control
2.
Serge Belongie, Jitendra Malik, and Jan Puzicha.
Shape
points would be required to represent a curve. However,
Matching and Object Recognition Using Shape Contexts.
such a system would be more complicated to implement,
IEEE Trans. Pattern Anal. Machine Intell., 24(2), April 2002.
and might have difﬁculty capturing sharp features, unless
To appear. 2
they are explicitly included in the curve representation. It
3.
Andrew Blake and Michael Isard. Active Contours. Springer,
would also be interesting to experiment with displacement
1998. 2
vectors9, 18. We speculate that they may give good results in
4.
P. J. Burt and E. H. Adelson. A multiresolution spline with
some cases.
application to image mosaics. ACM Transactions on Graphics,
2(4):217–236, October 1983. 8
Drawing systems. Our system can be viewed as a ﬁrst step
5.
Hong Chen, Ying-Qing Xu, Heung-Yeung Shum, Song-Chun
toward a larger system that generates many curves or begins
Zhu, and Nan-Ning Zheng.
Example-based Facial Sketch
with images or 3D geometry as input, such as in the work of
Generation with Non-parametric Sampling.
Proc. Interna-
Kalnins et al.18 and Jodoin et al.17. For example, our system
tional Conference on Computer Vision, 2001. 2
currently assumes that every input curve produces exactly
one output, and thus cannot learn styles with broken curves,
6.
Cassidy Curtis. Loose and Sketchy Animation. In SIGGRAPH
such as the overdrawing common in sketchy styles. Addi-
98: Conference Abstracts and Applications, page 317, 1998. 1
tional types of constraints could be added, in the form of
7.
Alexei A. Efros and William T. Freeman. Image Quilting for
new cost terms. For example, we could specify constraints
Texture Synthesis and Transfer. In Proceedings of SIGGRAPH
on interactions between curves or between a curve and an
2001, pages 341–346, 2001. 2, 4, 6
image.
8.
Alexei A. Efros and Thomas K. Leung. Texture Synthesis by
Non-parametric Sampling. In IEEE International Conference
Additional features. Additional features of the curve —
on Computer Vision, pages 1033–1038, September 1999. 2,
such as pressure, thickness, and brush texture — could be
3, 4, 5
easily added to the synthesis as additional features in the
9.
Adam Finkelstein and David H. Salesin.
Multiresolution
curves A , B (although more training data may be required in
Curves. Proceedings of SIGGRAPH 94, pages 261–268, July
order to capture greater variability). Closed curves can eas-
1994. 2, 10
ily be modeled by allowing neighborhood sampling to wrap
10.
W. T. Freeman, E. C. Pasztor, and O. T. Carmichael. Learn-
around the curve.
ing Low-Level Vision. Intl. J. Computer Vision, 40(1):25–47,
2000. 2
Animation. Curve analogies can be used for deﬁning line
11.
William T. Freeman, Joshua B. Tenenbaum, and Egon Pasz-
styles for non-photorealistic animation. If stroke temporal
tor. An example-based approach to style translation for line
coherence is desired, then corresponding terms could be
drawings. Technical Report TR99-11, MERL, February 1999.
added to the energy functions. However, whereas ﬂickering
2
in full-frame painting styles may be very distracting, ﬂicker-
ing in line-drawn animation is often acceptable or even desir-
12.
Paul Harrison.
A Non-Hierarchical Procedure for Re-
able, since it gives a sense of life to otherwise still characters.
Synthesis of Complex Textures. Proc. WCSG, 2001. 2, 6
Many classic line-drawn animations did not use temporal co-
13.
Aaron Hertzmann, Charles E. Jacobs, Nuria Oliver, Brian Cur-
herence of the form that has recently become a goal of NPR
less, and David H. Salesin. Image Analogies. Proceedings of
research.¶
SIGGRAPH 2001, pages 327–340, 2001. 1, 2, 4, 6
14.
Berthold K. P. Horn. Closed Form Solution of Absolute Orien-
3D signal processing. We believe that the curve analogies
tation using Unit Quaternions. Journal of the Optical Society
approach could be extended to learning probability distribu-
of America, 4(4):629–642, April 1987. 5
tions over 3D surfaces, for recognition and 3D signal pro-
15.
Berthold K. P. Horn, Hugh M. Hilden, and Shahriar Negah-
cessing applications.
daripour. Closed Form Solution of Absolute Orientation us-
ing Orthonormal Matrices. Journal of the Optical Society of
America, 5(7):1127–1135, July 1988. 5
References
16.
Siu Chi Hsu and Irene H. H. Lee.
Drawing and Anima-
1.
Michael Ashikhmin. Synthesizing Natural Textures. 2001
tion Using Skeletal Strokes. In Proceedings of SIGGRAPH
ACM Symposium on Interactive 3D Graphics, pages 217–226,
’94 (Orlando, Florida, July 24–29, 1994), Computer Graphics
March 2001. 4
Proceedings, Annual Conference Series, pages 109–118, July
1994. 2
17.
Pierre-Marc Jodoin, Emric Epstein, Martin Granger-Pichi, and
¶
Victor Ostromoukhov. Hatching by Example: a Statistical Ap-
To name a few: “The Big Snit” by Richard Condie, “Special De-
proach. NPAR 2002 : Second International Symposium on Non
livery” by John Weldon and Eunice Macaulay, “Billy’s Balloon” by
Photorealistic Animation and Rendering, June 2002. To ap-
Don Hertzfeldt, the animations of Bill Plympton, “Crac” by Frédéric
pear. 2, 10
Back, and the feature ﬁlm “My Neighbors, the Yamadas” by Isao
Takahata.
18.
Robert D. Kalnins, Lee Markosian, Barbara J. Meier,
10

Michael A. Kowalski, Joseph C. Lee, Philip L. Davidson,
Matthew Webb, John F. Hughes, and Adam Finkelstein.
WYSIWYG NPR: Drawing Strokes Directly on 3D Models.
In Proceedings of SIGGRAPH 2002, July 2002. To appear. 2,
10
19.
Leif Kobbelt and Peter Schröder. A multiresolution framework
for variational subdivision. ACM Transactions on Graphics,
17(4):209–237, October 1998. 2
20.
Lin Liang, Ce Liu, Ying-Qing Xu, Baining Guo, and Heung-
Yeung Shum. Real-time texture synthesis by patch-based sam-
pling. ACM Transactions on Graphics, 20(3):127–150, July
2001. 4
21.
Lee Markosian, Michael A. Kowalski, Samuel J. Trychin,
Lubomir D. Bourdev, Daniel Goldstein, and John F. Hughes.
Real-Time Nonphotorealistic Rendering. In SIGGRAPH 97
Conference Proceedings, pages 415–420, August 1997. 1
22.
Lee Markosian, Barbara J. Meier, Michael A. Kowalski, Lor-
ing S. Holden, J. D. Northrup, and John F. Hughes. Art-based
Rendering with Continuous Levels of Detail. In NPAR 2000 :
First International Symposium on Non Photorealistic Anima-
tion and Rendering, pages 59–66, June 2000. 2
23.
Michael P. Salisbury, Sean E. Anderson, Ronen Barzel, and
David H. Salesin. Interactive Pen–And–Ink Illustration. In
Proceedings of SIGGRAPH ’94 (Orlando, Florida, July 24–
29, 1994), pages 101–108, July 1994. 1, 2
24.
Shinji Umeyama. Least-squares estimation of transformation
parameters between two point patterns. IEEE Trans. Pattern
Anal. Machine Intell., 13(4):376–380, 1991. 5
25.
Li-Yi Wei and Marc Levoy.
Fast Texture Synthesis Using
Tree-Structured Vector Quantization.
Proceedings of SIG-
GRAPH 2000, pages 479–488, July 2000. 2, 8
26.
Li-Yi Wei and Marc Levoy. Texture Synthesis Over Arbitrary
Manifold Surfaces. Proceedings of SIGGRAPH 2001, pages
355–360, August 2001. 2
27.
C. I. Yessios. Computer drafting of stones, wood, plant and
ground materials. In Computer Graphics (Proceedings of SIG-
GRAPH 79), volume 13, pages 190–198, August 1979. 1
28.
Lexing Ying, Aaron Hertzmann, Henning Biermann, and De-
nis Zorin. Texture and Shape Synthesis on Surfaces. Proc.
12th Eurographics Workshop on Rendering, pages 301–312,
June 2001. 2
29.
Song-Chun Zhu. Embedding Gestalt Laws in Markov Random
Fields. IEEE Trans. Pattern Anal. Machine Intell., 21(11),
November 1999. 2
11

Figure 6: Shape synthesis with translation invariance. The
example shape is shown on the left, and the output shape on
:
::
the right.
A
A'
:
::
:
:
B
B'
Figure 9: Shape analogy with translation invariance. The B
curve is generated by analogy to the other curves.
(a)
Figure 7: Shape analogy with rotational invariance. The B
curve is generated by analogy to the other curves.
(b)
(c)
Figure 10: Shape analogy with constraints. (a) Source A
and A pair. (b,Left) Source B curves. (b,Right) Applying
Figure 8: Source A and A curves used in Figure 11.
the style to each curve separately produces curves that do
not intersect at their endpoints. (c,Left) Source B curves
with hard constraints. (c,Right) Applying the style with con-
straints forces the curves to meet at their endpoints.
12

Figure 11: Combining shape analogies. The curves in the left image were processed with styles learned from the examples in
Figure 8, and composited manually to produce the right image. The character and texture of the example strokes are transfered
to the source drawing. Rotational invariance was used for all examples.
:
::
:
Figure 12: Shape analogy with rotational invariance, applied to handwriting. The output curve is generated by analogy to the
other curves. Here, we create a ﬁlter that makes the shape less rigid and more ﬂowing. For example, the bend at the beginning
of the “M” is copied to the beginning of the “B.”
13