This is a blog post I made to the "mollyrocket.com" forums. Those forums were rendered unusable during a software changeover, so I attempted to reconstruct the post here from the slightly broken version still accessible on the forum.

L-Systems Considered Harmful

Sean Barrett, 2007-2009??
email: sean at this domain
twitter: nothings

I was looking at the paper Procedural Modeling of Cities by Parish and Mueller (which was referenced by the guy making Introversion’s game Subversion). They use an L-system to describe their road-building algorithm. I had looked at the paper before and my eyes had glazed over when L-systems were mentioned and I mostly blew off that part of the paper. But this time I decided to wade through and figure out what that part of the paper was describing.

L-systems are apparently de rigeur in procedural modelling, especially if what you’re modelling has branches. Here is what they say in their introduction:

Grammar-based generation of models (especially L-systems) are employed in computer graphics mainly to create plant geometry [26]. L-systems have evolved into very sophisticated and powerful tools [20, 27] and have been extensively used in the modeling of plant ecosystems described in [8].
and here is what they say in their abstract:
We propose a system using a procedural approach based on L-systems to model cities. [. . .] For the creation of a city street map, L-systems have been extended with methods that allow the consideration of global goals and local constraints and reduce the complexity of the production rules.
So, that’s admirable, to try to augment L-systems to simplify the L-system-y part of the thing, right? Well, only if you should be using L-systems in the first place.

Here’s the extended L-system for road-building, directly from the paper (’w’ is omega, and ‘e’ is epsilon):

w: R(0, initialRuleAttr) ?I(initRoadAttr, UNASSIGNED)
p1: R(del, ruleAttr) : del<0 -> e
p2: R(del, ruleAttr) > ?I(roadAttr,state) : state==SUCCEED
    { globalGoals(ruleAttr,roadAttr) creates the parameters for: 
          pDel[0-2], pRuleAttr[0-2], pRoadAttr[0-2] }
    -> +(roadAttr.angle)F(roadAttr.length)
      B(pDel[1],pRuleAttr[1],pRoadAttr[1]),
      B(pDel[2],pRuleAttr[2],pRoadAttr[2]),
      R(pDel[0],pRuleAttr[0]) ?I(pRoadAttr[0],UNASSIGNED)[i]
p3: R(del, ruleAttr) > ?I(roadAttr, state) : state==FAILED -> e
p4: B(del, ruleAttr, roadAttr) : del>0 -> B(del-1, ruleAttr, roadAttr)
p5: B(del, ruleAttr, roadAttr) : del==0 -> [R(del, ruleAttr)?I(roadAttr, UNASSIGNED)]
p6: B(del, ruleAttr, roadAttr) : del<0 -> e
p7: R(del, ruleAttr) < ?I(roadAttr,state) : del<0 -> e
p8: ?I(roadAttr,state) : state==UNASSIGNED
    { localConstraints(roadAttr) adjusts the parameters for:
        state, roadAttr}
    -> ?I(roadAttr, state)
p9: ?I(roadAttr,state) : state!=UNASSIGNED -> e
(Don’t bother closely reading this; I have a clearer rewrite down a bit.)

That’s pretty complicated, so it seems like it must be a good idea that they moved some of the complexity into those globalGoals and localConstraints functions, since otherwise it would be that much worse, right?

Now, if you don’t know what an L-system is, here’s a quick explanation. It’s a formalism defined around string rewriting. You begin with the string defined by w. Then, you rewrite the whole string according to the productions, in parallel, producing a new string. Then you repeat on the new string. L-systems support probabilistic substitution, and they can make use of left-and-right contexts around given tokens of the part of the string they’re rewriting.

But, there’s an issue here: string rewriting for something huge like a city would be slow. The formalism is based around string rewriting, but you might not want to actually use string rewriting. So what should you actually do? Or, to get more directly to my point, what do authors describing systems like the above expect readers to do if they implement it themselves?

I see three possibilities:

  • implement it as string rewriting
  • implement it as module rewriting (a list or tree of heterogenous nodes, which are still rewritten in parallel as with the strings)
  • the L-system is just a formalism; flatten it out to the simplest code possible, not L-systemy code

    I have no idea what the authors’ intent is. I assume it’s one of the first two, presumably the second. I assume it’s not the third, because the third is a bunch of work, and if you think that’s the way to do it, why wouldn’t you do that work in the paper, rather than make everyone do it themselves?

    Of course, a paper can only be so long, and the L-system formulation is small and precise, so maybe you assume #3 but don’t spell it out because that’s busywork that wastes space in the paper. But, I’d argue, at least for the above formula, I don’t think the authors ever spelled it out. In fact, it turns out that the L-system above isn’t even really an L-system at all, as is obvious if you try to reduce it. Maybe implementing it as an L-system made it easier for them to prototype or something, but if there was any value there, it’s gone now.

    To demonstrate this, I am going to apply mode #3, and we will see where it takes us. (Or you can skip to the end of this post and see the punchline.)

    Let me first rewrite the rules in something a little less compact, more readable, but largely unchanged. I’m going to re-order the production rules and put the most complicated ones last just to make it a tad more readable (since all the production rules are applied in parallel, the order doesn’t matter, but I’ll keep the original numbering for clarity).

    Initial state:
    
    [ road(0, ra) query(qa, UNASSIGNED) ]
    
    rules:
    
    1. delete road(x,ra) if x < 0
    3. delete road(x,ra) if followed by query(qa, FAILED)
    4. decrement x in branch(x, qa, ra) if x > 0
    5. replace branch(x,ra,qa) with subtree [road(x, ra) query(qa, UNASSIGNED)] if x == 0
    6. delete branch(x,ra,qa) if x < 0
    7. delete query(qa,state) if preceded by road(x,ra) and x < 0
    9. delete query(qa,state) if state != UNASSIGNED
    8. use local constraints(qa) to update query(qa,UNASSIGNED) to query(nqa,state)
    2. given road(x,ra) query(qa,SUCCEED):
         run globalGoals(qa,ra) to compute pdel[3], nqa[3], nra[3]
         replace road(x,ra) with:
             segment(ra.angle, ra.length)
             branch(pdel[1], nra[1], nqa[1])
             branch(pdel[2], nra[2], nqa[2])
             road(pdel[0], nra[0])
             query(nqa[0], UNASSIGNED)
    
    Now, if you work through what’s going on here, you’ll find that a lot of the stuff here is doing various kinds of bookkeeping necessary to make things work in an L-system. But even within the L-system form, there’s room to simplify this. For example, we can short-circuit the query-deletion rule from 9 into 3 (when it’s FAILED) and 2 (when it’s SUCCEED):
    3. delete road(x,ra) query(ini_query, FAILED)
    2. given road(x,ra) query(qa,SUCCEED):
          run globalGoals(qa,ra) to compute pdel[3], nqa[3], nra[3]
          replace both tokens
              segment(ra.angle, ra.length)
              branch(pdel[1], nra[1], nqa[1])
              branch(pdel[2], nra[2], nqa[2])
              road(pdel[0], nra[0])
              query(nqa[0], UNASSIGNED)
    
    Note that in the original #2, it looked for a road/query pair and only replaced the road, leaving rule #9 to delete the query. So now I’ve just made it replace both, and rule #9 goes away.

    There are similar things that can be done, but we quickly notice now that road&query are always a pair anyway; for example production #1 and production #7 make it look like more flexible updating is happening, but in fact they’re just always updated as a pair. So we replace them with roadq(x,ra,qa,state):

    Initial state:
    
    roadq(0, ra, qa, UNASSIGNED)
    
    rules:
    
    1. delete roadq(x,ra,qa,state) if x < 0
    3. delete roadq(x,ra,qa,FAILED)
    4. decrement x in branch(x, qa, ra) if x > 0
    5. replace branch(x,ra,qa) with subtree [roadq(x, ra, qa, UNASSIGNED)] if x == 0
    6. delete branch(x,ra,qa) if x < 0
    7. [implied by 1]
    9. [folded into 2,3]
    8. use local constraints(qa) to update roadq(x,ra,qa,UNASSIGNED) to roadq(x,ra,nqa,state)
    2. given roadq(x,ra,qa,SUCCEED):
           run globalGoals(qa,ra) to compute pdel[3], nqa[3], nra[3]
           replace roadq(x,ra, qa, SUCCEED) with:
               segment(ra.angle, ra.length)
               branch(pdel[1], nra[1], nqa[1])
               branch(pdel[2], nra[2], nqa[2])
               roadq(pdel[0], nra[0], nqa[0], UNASSIGNED)
    
    Note that many of these simplifications are all still occuring in parallel on the same l-system update round as before. However, because no rules have any contexts anymore, we don’t actually have to do this; we can delete things a round earlier without affecting the results.

    Since branch(x,. . .) only decrements x if it’s greater than 0, the only way it can become negative and invoke rule 6 is if it’s created negative. Since it’s never used for context, we delete rule 6 and replace rule 2:

    2. given roadq(x,ra,qa,SUCCEED):
           run globalGoals(qa,ra) to compute pdel[3], nqa[3], nra[3]
           replace road(x,ra) with:
               segment(ra.angle, ra.length)
               if (pdel[1] >= 0) branch(pdel[1], nra[1], nqa[1])
               if (pdel[2] >= 0) branch(pdel[2], nra[2], nqa[2])
               roadq(pdel[0], nra[0], nqa[0], UNASSIGNED)
    
    (I assume the reason they do this the other way is that you can’t express this conditionally-add-something in regular L-systems, that it’s normal to use a rewrite rule instead. But it wouldn’t be hard to add this feature to the L-system formalism, not if you’re adding this function call stuff anyway!)

    Similarly, the roadq deletion rule in #1 can only occur if the roadq is installed negative. Because we’ll never run branch with a negative, rule #5 will never install a negative roadq, so the only place it can appear is in rule 2, so we remove rule #1 and update the last line of rule 2:

    if (pdel[0] >= 0) roadq(pdel[0], nra[0], nqa[0], UNASSIGNED)
    

    Here’s the current rule set:

    3. delete roadq(x,ra,qa,FAILED)
    4. decrement x in branch(x, qa, ra) if x > 0
    5. replace branch(x,ra,qa) with subtree [roadq(x, ra, qa, UNASSIGNED)] if x == 0
    8. use local constraints(qa) to update roadq(x,ra,qa,UNASSIGNED) to roadq(x,ra,nqa,state)
    2. given roadq(x,ra,qa,SUCCEED):
           run globalGoals(qa,ra) to compute pdel[3], nqa[3], nra[3]
           replace road(x,ra) with:
             segment(ra.angle, ra.length)
             if (pdel[1] >= 0) branch(pdel[1], nra[1], nqa[1])
             if (pdel[2] >= 0) branch(pdel[2], nra[2], nqa[2])
             if (pdel[0] >= 0) insert roadq(pdel[0], nra[0], nqa[0], UNASSIGNED)
    
    Rule #3 can only be invoked by rule #8, so we can delete rule 3 and embed it in 8:
    8. given roadq(x,ra,qa,UNASSIGNED):
           run localConstraints(ra,qa) to compute nqa,state
           replace roadq(x,ra,qa,UNASSIGNED) with
              if (state != FAILED) roadq(x,ra,nqa,state)
    
    At this point, the end of rule #8 is the only place that can invoke rule #2, so we move it in:
    4. decrement x in branch(x, qa, ra) if x > 0
    5. replace branch(x,ra,qa) with subtree [roadq(x, ra, qa, UNASSIGNED)] if x == 0
    8. given roadq(x,ra,qa,UNASSIGNED), delete it and:
           run localConstraints(ra,qa) to compute nqa,state
              if (state == SUCCEED) 
                 { run globalGoals(qa,ra) to compute pdel[3], nqa[3], nra[3]
                   insert segment(ra.angle, ra.length)
                   if (pdel[1] >= 0) insert branch(pdel[1], nra[1], nqa[1])
                   if (pdel[2] >= 0) insert branch(pdel[2], nra[2], nqa[2])
                   if (pdel[0] >= 0) insert roadq(pdel[0], nra[0], nqa[0], UNASSIGNED)
    

    Now, if we change from turtle-style motion with nesting at [] (i.e. [] means save/restore a turtle state) to assuming explicit coordinates are stored in the roadq()s and segments(), we can get rid of the [] calls entirely. In that case, we can migrate the decrement from branch to roadq, and get rid of branch().

    4. decrement x in roadq(x, qa, ra, state) if x > 0
    8. given roadq(x,ra,qa,UNASSIGNED) with x = 0, delete it and:
           run localConstraints(ra,qa) to compute nqa,state
              if (state == SUCCEED) 
                 run globalGoals(qa,ra) to compute pdel[3], nqa[3], nra[3]
                 insert segment(ra.angle, ra.length)
                 if (pdel[1] >= 0) insert roadq(pdel[1], nra[1], nqa[1], UNASSIGNED)
                 if (pdel[2] >= 0) insert roadq(pdel[2], nra[2], nqa[2], UNASSIGNED)
                 if (pdel[0] >= 0) insert roadq(pdel[0], nra[0], nqa[0], UNASSIGNED)
    
    At this point we can delete the fourth state parameter to roadq().

    Now, at this point, we still have an L-system, just one that’s been augmented with this support that when you write out the replacement for an existing module, you can conditionally choose whether certain things are written, and one where we’ve removed the turtle-graphics style information and assumed it’s all tracked explicitly by globalGoals().

    But what should be clear is that, by adding these small L-system features and mechanically simplifying the rules, we’ve simplified it down to the point where it’s hardly an L-system at all.

    In fact, we can now express it as an actual pseudocode algorithm, rather than an L-system, which has the advantage that it’s clear what exactly it does (rather than the behavior being buried in interactions of productions), and understanding it doesn’t require understanding L-systems.

    So here is their actual algorithm:

    initialize list of roadqueries L to a single entry:
        roadq(0,ra,qa)
    
    initialize segment list S to empty
    
    until L is empty
       for each entry R = roadq(x,ra,qa) in L
          if x > 0
             decrement x
          else
             L.delete( R )
             compute (nqa,state) = localConstraints(qa)
             if (state == SUCCEED)
                compute (pdel[3],nqa[3],nra[3]) = globalGoals(qa,ra)
                S.add ( segment(ra) );
                if (pdel[1] >= 0) L.add( roadq(pdel[1], nra[1], nqa[1]) )
                if (pdel[2] >= 0) L.add( roadq(pdel[2], nra[2], nqa[2]) )
                if (pdel[0] >= 0) L.add( roadq(pdel[0], nra[0], nqa[0]) )
    
    Switching to this algorithm in their paper would have the aforementioned benefits of clarity, but would also allow them to skip all the explanation in their paper where they explain how the various productions implement the various hacks to allow emulating a regular algorithm, so it would actually be shorter even if it wasn’t already shorter than the L-system (which it is).

    We can now consider changing the algorithm. It currently spends a lot of time iterating the list and decrementing the x’s without doing any work (depending on the actual values of the xs). We can fix this. It’s possible that computing globalGoals() and localConstraints() will dominate the performance anyway (the list L may stay reasonably small, depending on the pdel values -- the performance of the original system, which left the segments() in the same list, was O(N^2)-ish though) so it doesn’t need fixing. But fixing it lets us express it more simply anyway. What we want to do is go ahead and add a timestep t that we increment every time through the loop, then instead of storing pdel and decrementing it, we store t+pdel and wait until t reaches that value.

    Once we’ve got that, we can get rid of incrementing t and just use a priority queue:

    initialize priority queue Q with a single entry:
        r(0,ra,qa)
    
    initialize segment list S to empty
    
    until Q is empty
        pop smallest r(t,ra,qa) from Q
        compute (nqa, state) = localConstraints(qa)
        if (state == SUCCEED) {
            add segment(ra) to S
            compute (pdel[3], nqa[3],nra[3]) = globalGoals(qa,ra)
            if (pdel[0] >= 0) add r(t+1+pdel[0], nra[0], nqa[0]) to Q
            if (pdel[1] >= 0) add r(t+1+pdel[1], nra[1], nqa[1]) to Q
            if (pdel[2] >= 0) add r(t+1+pdel[2], nra[2], nqa[2]) to Q
    
    Now we can simplify the algorithm by having globalGoals directly add elements to Q, instead of returning the parameters and leaving the caller to create them, thus giving us the actual algorithm underlying the road-building system:
    initialize priority queue Q with a single entry:
         r(0,ra,qa)
    
    initialize segment list S to empty
    
    until Q is empty
        pop smallest r(t,ra,qa) from Q
            compute (nqa, state) = localConstraints(qa)
            if (state == SUCCEED) {
                add segment(ra) to S
                addZeroToThreeRoadsUsingGlobalGoals(Q, t+1, qa,ra)
    
    So, yeah. Nice L-system you got there, Mueller and Parish.
    Here are some of the additional comments I posted in the thread:
    perhaps an inefficient L-system implementation was just fast enough for their purposes, so they stopped thinking about it after they built it?
    That’s why I called out explicitly that the priority queue implementation isn’t necessarily for speed, but for clarity of description. A priority queue is a formalism, just like an L-system is, so it allows you to use a concise, precise description rather writing out a system that involves putting a delay on everything which you manually decrement.

    It shoulld be easy to guess from my last version of the algorithm exactly what they’re doing (they can use the ‘pdel’ field to decide what order to expand roads in, since presumably roads are blocked/stopped by other roads so the order matters).

    Also, it turns out that there’s apparently no magic in this L-system addressing the closedness itself, other than by doing that very blocking. None of that is easy to guess in the L-system version.

    In fact, let me quote that bit you said above:

    perhaps [. . .] they stopped thinking
    That is exactly what I meant by the title. I think L-systems rot your brain; or, if not that, at least they make people not actually invest effort into thinking or communicating.

    Maybe they make more sense for plant synthesis. But I have a suspicion that most uses out there would turn out to be using similarly degenerate subsets that would be better expressed in other ways.