#!/usr/local/bin/perl
# Copyright © 2001 by James F. Carter
# 2001-01-15, Perl-5.006.
# Produces a PostScript file of a Penrose tiling (fat and thin rhombi).
# Permission to use, copy, modify, distribute, and sell this
# software and its documentation for any purpose is hereby granted
# without fee, provided that the above copyright notice appear in all
# copies and that both that copyright notice and this permission notice
# appear in supporting documentation. James F. Carter makes no
# representations about the suitability of this software for any purpose.
# It is provided "as is" without express or implied warranty.
# See http://www.law.gwu.edu/facweb/claw/penrose.htm for the story of the
# notorious Kimberly-Clark affair. Sir Roger Penrose asserts copyright on the
# design, which has been assigned to Pentaplex Ltd., Brighouse, UK
# . The Penrose tiling produced by this
# program is intended as a background on a web page. The terms described to
# me for thus using it are:
# 1. Your site is not and will not be used for commercial gain.
# 2. You acknowledge on the page (and any others showing such a tiling) that
# the tiling is illustrated with the permission of Pentaplex Ltd, Brighouse, UK.
# 3. A link is made available to the pentaplex.com site.
# The fill color of the rhombi is hardwired in @colors (the first instance),
# about 90% of the way through. Edges can also be colorized, though I decided
# to deactivate that feature.
# How to use:
# penrose.pl -r -s 12.5 -X 0 -Y 0 -x 800 -y 600 -n 50000 > penrose.ps
# gs -q -g800x600 -sDEVICE=png16m -sOutputFile=penrose.png - < penrose.ps
# mogrify -format gif penrose.png
# Reference: http://www.math.ubc.ca/~robles/tiling/penrose.html
# This is a tutorial on tiling in general, including Penrose tiling. It has
# the original 2D rules, as well as discussion de Bruijn's integer lattice
# formulation.
# Rules for Penrose tiling:
# A. The narrow rhombus has angles of 36 and 144 degrees. Put a dot on one
# 144 deg vertex and on the two opposite edges put arrows pointing toward
# the other 144 deg vertex.
# B. The wide rhombus has angles of 72 and 108 degrees. Put a dot on one 72
# deg vertex and on the two opposite edges put arrows pointing away from
# the other 72 deg vertex.
# C. Around any point either all vertices must be dotted or must be
# undotted.
# D. Across any edge the two tiles must have no arrows or must have arrows
# pointing in the same direction.
# Vertex star patterns, i.e. patterns of rhombi around any common vertex:
# (The letters are apparently standard nomenclature.) * indicates center dot.
# A*. 5 wide ones with dots at the center.
# B. 5 wide ones with the undotted 72deg vertex at the center.
# C. Like type B with one wide rhombus replaced by 2 narrow ones.
# D. Like type B with 2 (non-adjacent) wide rhombi replaced by pairs of
# narrow ones.
# E. Open box: in the 144deg space between 2 wide rhombi, put (in order)
# a narrow one, a wide one, and a narrow one. Center vertex cannot
# be dotted.
# F*. Like type A with 2 adjacent wide rhombi replaced by one narrow one
# (sideways). Center vertices are dotted only.
# G. Necker cube: 2 wide rhombi with the 144deg space filled by a
# narrow one. Center cannot be dotted, and "together" vertices of the
# wide rhombi must be dotted.
# H*. Converse of G: 2 narrow rhombi with the 72deg space filled by
# a wide one. Center vertices are dotted only.
# At most of the vertex stars, most of the edges have the same polarity.
# Several have exceptions of polarity 1 (dot at opposite end): C,D,G.
# Type E has 2 1's, 1 2 (arrow in) and 2 3's (arrow out). This precludes
# summarizing the edge polarity at each vertex.
# Timing (on Pentium III 500 MHz) as a function of the scale factor (-s) using
# the default page dimensions.
# Sub-optimal choice of which point to surround next: global FIFO.
# Scale(-s) Rhombi Added Rhombi Kept Time (secs)
# 200 54 40 0.37
# 100 108 80 0.70
# 50 508 205 2.93
# 25 10366 687 58.27
# Better choice: most complete first; FIFO within completeness classes.
# 100 74 71 0.33
# 50 334 205 1.14
# 25 1352 700 4.25
# 12.5 7005 2564 20.76
# Time is roughly proportional to rhombi added. The ratio of rhombi added
# vs rhombi kept is an increasing function perhaps proportional to
# log(rhombi kept). Rhombi kept is approximately
# proportional to scale^(-2), but there are many rhombi along the perimeter.
package main;
use Carp;
use Getopt::Std;
our $opt_s = 100; # Scale factor, length of all the lines, in PS points
# The geometry is created with scale = 1 initially.
our $opt_x = 468; # Total width of pattern in points: 6.5in wide
our $opt_y = 648; # Total height of pattern in points: 9in high
our $opt_X = 72; # Origin of pattern in points
our $opt_Y = 72; # Origin of pattern in points
our $opt_n = 10000; # Max number of points allowed
our $opt_d; # Debug switch
our $opt_v; # Verbosity switch
# opt_r Print rhombi (filled polygons)
# opt_e Print edges (lines)
&getopts("s:x:y:X:Y:n:ervd");
die "You need one or both of -e (print edges) or -r (print rhombi).\n"
unless $opt_e || $opt_r;
$Xw = $opt_x / $opt_s;
$Yw = $opt_y / $opt_s;
our $halt; # Set to \%Edge or \%Point or 1 to halt processing.
# Wide rhombus: acute vertex is 72 deg, obtuse is 108 deg.
# Narrow rhombus: obtuse vertex is 144 deg, acute is 36 deg.
our $N2 = 5;
our $N = 2*$N2; # Number of rotational divisions
our $pi5 = 8 * atan2(1, 1) / $N; # 2pi/10 = 36 degrees.
# The terminology "spoke from $a" refers to one edge of a rhombus one of whose
# vertices is point $a. Each point has several spokes radiating from it, all
# of the same length.
# "Canonical" coordinates are a single integer: 10^$c * $x + 10*$y
# Polarity codes used in several places:
# 0 Dot at this vertex, or at the start point of an edge
# 1 Non-dot vertex (corner), or dot at the endpoint of an edge.
# 2 Arrows toward start point
# 3 Arrows away from start point, toward end point
our %polvtx = ('', qw(undef 0 dotted 1 corner 2 arrow-in 3 arrow-out 4 bogus));
# Generic arguments to $this->print($level):
# $level 0 or undef: just identifies the item (no \n)
# 1: brief debug info (1 line, no \n)
# 2: more debug info (multiple lines with ending \n)
# It returns a string; you print it yourself.
# ===================================
package Point;
# A Point object has all the necessary info about one point. Data members:
# x Its X coordinate (exact)
# y Its Y coordinate (exact)
# can Canonical integer representation of the coordinates
# dest A hash whose key is the canonical coord. of another point,
# to which an edge extends from this one. Value is a ref. to
# an Edge object. (Being phased out?)
# azs A hash whose key is azimuths (0..9). Value is a ref. to an
# Edge object.
# tiles A string of 10 bytes, 1 per azimuth. The bytes are:
# . unoccupied
# W start of wide rhombus
# w wide rhombus continued
# N Start of narrow rhombus
# n narrow rhombus continued
use Carp;
our %real; #Key = canonical coords; value = ref. to point with that coord.
# Creates a point. Args: x, y (to go in the blessed object).
sub Point::new {
my($class) = shift;
my($this) = bless({
dest => { },
azs => { },
tiles => '.' x $N
}, $class);
@{$this}{qw(x y)} = @_;
$this->{can} = $this->canon();
$real{$this->{can}} = $this;
$this;
}
# Converts exact coordinates to canonical integer coords.
# The same point, when approached from different directions, can have slightly
# different coords due to roundoff error, and when rounded could produce a
# different canonical integer. The subroutine checks if @a is in %real,
# looking at neighboring points if not found. Arguments:
# $this Only x and y need to be filled in.
# Returns the canonical integer coordinate, e.g. (8,2) -> 80020
# "Canonical" coordinates are a single integer: 10^$c * $x + 10*$y
$log10 = log(10);
$Clg = int(2.55 + log(($main::Xw > $main::Yw) ? $main::Xw : $main::Yw)/$log10);
$C = int(0.5 + exp($log10 * $Clg));
$Cfmt = sprintf "%%%dd", 2*$Clg; #Format for printing canonical coords
sub Point::canon {
my($this) = @_;
my($xi, @xi, @dC);
my($dC) = $C;
foreach $x (@{$this}{qw(x y)}) {
push(@xi, ($xi = int(10*$x)));
if ($x >= 0) {
if ($x < 0.1*$xi + 0.000001) {
push(@dC, -$dC);
} elsif ($x > 0.1*$xi + 0.099999) {
push(@dC, $dC);
}
} else {
if ($x > 0.1*$xi - 0.000001) {
push(@dC, $dC);
} elsif ($x < 0.1*$xi - 0.099999) {
push(@dC, -$dC);
}
}
$dC = 1;
}
push(@dC, $dC[0] + $dC[1]) if @dC > 1;
my($ac) = $C*$xi[0] + $xi[1];
unless (exists($real{$ac})) {
foreach $dC (@dC) {
next unless exists($real{$ac + $dC});
$ac += $dC;
last;
}
}
$ac;
}
# Convert a vector to its direction.
# \%aa Base point of an edge ($this)
# \%bb Endpoint of an edge
# Returns: Integer in 0..9 such that $aa->endpt($bb, $result, xx) would
# reproduce the vector.
sub Point::myatan {
my($aa, $bb) = @_;
my($ing) = atan2($bb->{y} - $aa->{y}, $bb->{x} - $aa->{x})/$pi5;
$ing += $N if $ing < 0;
int($ing + 0.5) % $N;
}
# Like myatan, find the azimuth to a neighbor of $this that *does* have an
# edge going to it.
# $this Base point of the edge
# $bc Canonical coordinate of the neighbor
# Returns: Azimuth to that neighbor, looked up from saved data.
sub Point::az {
my($this, $bc) = @_;
defined(${$this}{dest}{$bc}) ? ${$this}{dest}{$bc}->az($this) :
$this->myatan($real{$bc});
}
# Is this a dotted vertex? Returns 1 if so. All edges out of a dotted vertex
# have polarity 0. The vertex must have at least one edge for this to work.
sub Point::isdot {
!(values %{$_[0]{dest}})[0]->pol($_[0]);
}
# Returns the neighbor (via an existing edge) in a given direction.
# $this Base point
# $az Azimuth, integer from 0..9, to be multiplied by pi/5.
# Returns: The existing point in that direction 1 unit away from $this.
sub Point::spoke {
my($this, $az) = @_;
defined($this->{azs}{$az}) ?
$this->{azs}{$az}->to($this) : undef;
}
# Removes a point. You should remove all the edges first.
sub Point::remove {
my($this) = @_;
delete $real{$this->{can}};
}
# Produces the endpoint of an edge. The resulting point is not caused to
# "exist" by being added to %real or $this->{dest}. Args:
# $this Base point
# $ing Azimuth, integer from 0..9, to be multiplied by pi/5.
# Returns: A new point 1 unit away from $this in the requested direction.
# The official representation is NOT returned whether or not
# it exists. It has x, y and can members only.
sub Point::endpt {
my($this, $ing) = @_;
my($b) = bless({
x => $this->{x} + cos($pi5 * $ing),
y => $this->{y} + sin($pi5 * $ing),
dest => { },
tiles => '.' x $N,
}, ref($this));
$b->{can} = $b->canon();
return $b;
}
# Produces the endpoint of an edge and adds it as an actual point. Args:
# $this Base point
# $b Endpoint produced by $this->endpt()
# $az Azimuth, integer from 0..9, to be multiplied by pi/5.
# $pol Polarity of the edge.
# Returns: A new point 1 unit away from $this in the requested direction.
# If the point already exists its official representation is
# returned. Whether or not existing, an edge is enrolled from
# the base point to the endpoint.
sub Point::addpt {
my($this, $b, $az, $pol) = @_;
$az = ($az + $N) % $N; #Make sure it's in 0..9
my($azopp) = ($az + $N2) % $N; #Azimuth in opposite direction
# Register the new endpoint, or use the registered copy.
my($bc) = $b->{can};
if (exists($real{$bc})) {
$b = $real{$bc};
} else {
$real{$bc} = $b;
}
# Check for placing an edge different from the existing one.
my($e);
if (exists($this->{azs}{$az})) {
$e = $this->{azs}{$az};
if ($pol != $e->pol($this)) {
print STDERR "addpt tried to replace an edge with different polarity: new = ",
$polvtx{$pol}, "\n\t", $e->print(1);
$halt = $e->flip($this);
}
}
# Add an edge to the new endpoint (if not there already).
my($ac) = $this->{can};
$e = Edge->new($this, $b, $pol, $az);
$b->{dest}{$ac} = $e unless exists($b->{dest}{$ac});
$b->{azs}{$azopp} = $e unless exists($b->{azs}{$azopp});
$this->{dest}{$bc} = $e unless exists($this->{dest}{$bc});
$this->{azs}{$az} = $e unless exists($this->{azs}{$az});
return $b;
}
# Adds or deletes a tile at a vertex, in the {tiles} string.
# $this Point (vertex) which is one corner of the tile.
# $az Azimuth of the A-B edge of the new tile. It must exist.
# $w Type of tile (w or n or '.'). Dot is for deleting an
# existing tile.
# Returns A string, which the caller could assign to $this->{tiles}.
sub addtile {
my($this, $az, $w) = @_;
my($e) = $this->{azs}{$az};
my($type) = $w . $e->pol($this);
my($daz) = $Rhomb::rhangle{$type};
my($tile) = uc($w) . ($w x ($daz-1));
my($res) = $this->{tiles};
substr($res, $az, $daz) = $tile;
my($j) = length($res) - $N;
if ($j > 0) {
substr($res, 0, $j) = substr($res, $N);
substr($res, $N) = '';
}
$res;
}
# Vertex star patterns. Key is star code letter, value last byte is 'D' for
# a dotted vertex, 'U' for undotted, preceeded by the potential content of
# Point->{tiles} (q.v. for codes) repeated twice.
%pickstar_OBSOLETE = qw(
A dWwWwWwWwWwWwWwWwWwWw
B uWwWwWwWwWwWwWwWwWwWw
C uWwWwWwWwNNWwWwWwWwNN
D uWwNNWwNNWwWwNNWwNNWw
E uWwwWwwNWwNWwwWwwNWwN
F dWwWwWwNnnnWwWwWwNnnn
G uWwwWwwNnnnWwwWwwNnnn
H dNnnnNnnnWwNnnnNnnnWw
);
# $vertexstar{$point->{tiles}} contains the letters of the vertex stars which
# are possible for that configuration of tiles, in all 10 cyclic shifts, and
# with each tile replaced by dots in a binary pattern.
{
my(@ptns) = ( #Number of variants ($N * (1< [qw(Ww Ww Ww Ww Ww)], # 320
B => [qw(Ww Ww Ww Ww Ww)], # 320 doesn't distinguish dot-in or out
C => [qw(Ww Ww Ww Ww N N)], # 640
D => [qw(Ww N N Ww N N Ww)], # 1280
E => [qw(Www Www N Ww N)], # 320
F => [qw(Ww Ww Ww Nnnn)], # 160
G => [qw(Www Www Nnnn)], # 80
H => [qw(Nnnn Nnnn Ww)], # 80
); # 3200 total members of %vertexstar
my($type, $ptn, @dots, $imax, $i, $j, $vtx0, $vtx);
while (($type, $ptn) = splice(@ptns, 0, 2)) {
@dots = map {'.' x length($_)} @$ptn;
$imax = 1<<@$ptn;
foreach $i (0..($imax-1)) {
$vtx0 = join('', map {($i & (1<<$_)) ? $ptn->[$_] : $dots[$_]}
0..(@$ptn-1));
foreach $j (0..($N-1)) {
$vtx = substr($vtx0, $j) . substr($vtx0, 0, $j);
$vertexstar{$vtx} .= $type
unless index($vertexstar{$vtx}, $type) >= 0;
}
}
}
}
# Checks if a particular kind of rhombus can be added at the point. (Only
# checks this one vertex.)
# $this Point where the rhombus might go.
# $az Azimuth of its A-B edge (doesn't have to exist).
# $w Type of rhombus (w or n).
# Returns: Legal vertex star codes; undef if vertex is illegal; or 1 if
# the $az edge doesn't exist (yet).
sub Point::check {
exists($_[0]->{azs}{$_[1]}) ? $vertexstar{$_[0]->addtile($_[1], $_[2])} : 1;
}
# Picks which kind of rhombus is feasible at a given point.
# $this Point at which a rhombus is wanted.
# Returns a list:
# [0] Azimuth where rhombus goes, or undef if $this is finished, i.e.
# it has no place for a rhombus.
# [1] Type of rhombus (n or w), or undef (with a defined azimuth) if
# a rhombus is needed but none is feasible.
# [2] Alternative type, undef if only one type is feasible.
sub Point::pickpoly {
my($this) = @_;
# The next polygon to be placed will start at this azimuth,
# corresponding to the first dot (unoccupied azimuth), except
# wrapping backward if azimuth 0 is unoccupied. Quick exit if
# the point is actually already finished.
my($az);
$az = index($this->{tiles}, '.');
return () unless $az >= 0; #This vertex is finished
if ($az <= 0 && substr($this->{tiles}, -1, 1) eq '.') {
$this->{tiles} =~ /(\.*)$/;
$az = $N - length($1);
}
if (!exists($this->{azs}{$az})) { #{tiles} eq '....'
#There could only be one edge. Rescue it. This only happens
#on the 1st point.
$az = (sort {$a <=> $b} keys %{$this->{azs}})[0];
if (!defined($az)) {
#No edges at all. All the rhombi got deleted first.
#Caller will take care of tossing the point.
return (0); #Signal failure to place rhombus
}
}
my(@nxtile, $w, $ntile);
foreach $w (qw(w n)) {
next unless Rhomb->check($this, $az, $w);
push(@nxtile, $w);
}
# If we've reached an illegal configuration, there will be
# no keys in %nxtile. The caller will take care of it.
($az, @nxtile);
}
# File a Point on a list of points segregated by the number of un-done
# azimuths. Finished points (without dots) are skipped. A point could be on
# several lists without hurting anything.
# $this Point to go on the list.
# \%ptlist Keys are the number of dots in $this->{tiles} and values are
# lists of Point refs.
sub ptlist_OBSOLETE {
my($this, $ptlist) = @_;
my $j = $this->{tiles} =~ tr/././;
push(@{$ptlist->{$j}}, $this) if $j > 0;
}
# Packs the unoccupied space around $this with rhombi
# $this Point to be surrounded with rhombi. Points outside the clip
# area are bypassed. OK to call on a point which is already
# surrounded by rhombi.
# \@rhlist The new rhombi are appended to this list. All listed rhombi
# are legal.
# \%ptlist The new points are appended to these lists. "New" means is a
# vertex of an added rhombus; the program may not be completely
# efficient at excluding finished points or points already on
# the list. See method ptlist() for details.
# Returns: 1 normally (including bypassing a point outside the clip area,
# or if the given point is actually finished already),
# or undef if the point could not be finished (e.g. it ended up
# with an illegal vertex star configuration).
sub Point::surround {
my($this, $rhlist, $ptlist) = @_;
# Quick exit for a point that's already finished or is deleted.
return 1 unless exists($real{$this->{can}}) &&
index($this->{tiles}, '.') >= 0;
# Bypass points outside the clip area. "Inside" means actually
# inside the clip rectangle, or in a 1-unit square at each
# corner. The reason for the squares is, there was an edge
# that ran from one point outside the rectangle to another,
# cutting across a corner. There should have been a tile
# placed covering the corner, but there would not have been
# because both endpoints were "outside" the actual rectangle.
unless ((-0.05 < $this->{x} && $this->{x} < $main::Xw &&
-0.05 < $this->{y} && $this->{y} < $main::Yw) ||
((abs($this->{x}) <= 1 || abs($this->{x} - $main::Xw) <= 1) &&
(abs($this->{y}) <= 1 || abs($this->{y} - $main::Yw) <= 1))) {
return 1; #Outside clip area counts as success
}
# Start placing polygons.
my($b, $bc, $j, %newpt);
my($res) = 1; #Assume correct placement
POLY: {
my($az, @ppy) = $this->pickpoly();
undef $res if defined($az) && !@ppy; #Oops, needs rhombus but illegal
last unless @ppy > 0; #Exit if rhombus can't be placed
my($rh) = Rhomb->new($this, $az, @ppy);
push(@$rhlist, $rh);
foreach $b (@{$rh->{vtx}}) {
$newpt{$b->{can}} = $b;
}
redo;
}
foreach $b (values %newpt) {
$ptlist->add($b);
}
return $res;
}
# "Prints" a point. Args:
# $this Point to be printed.
# $full How much to print (0=ID, 1= 1 line debug, 2 = full with \n's).
# Full format shows each edge.
# $fmt Format for X and Y (there's a default if you omit it).
# Returns: A string.
sub Point::print {
my($this, $full, $fmt) = @_;
return defined($this) ? "$this" : "undef" unless ref($this) eq 'Point';
$fmt = '%7.2f' unless defined($fmt);
my($res) = sprintf "$Cfmt$fmt$fmt", @{$this}{qw(can x y)};
$res .= " `${$this}{tiles}'" if $full >= 1;
if ($full >= 2) {
my($d, $e, $bc, $fl);
$d = $this->{dest};
my(@dests) = sort {$d->{$a}{az} <=> $d->{$b}{az}} keys %$d;
foreach $bc (@dests) {
$e = $d->{$bc};
if (ref($e) eq 'Edge') {
$fl = ($e->{from}{can} eq $this->{can}) ? 'or' : 'fl';
$e = $e->flip($this);
$res .= "\n\t" . $e->print(1);
} else {
$res .= "\n\t$bc -> `$e'";
}
}
$res .= "\n";
} elsif ($full == 1) {
my($dest) = join(' ', keys %{$this->{dest}});
$dest = "(none)" unless $dest;
$res .= " -> $dest";
}
$res;
}
# ===================================
package Edge;
# An Edge is necessarily polarized, but its member functions yield values
# appropriate to a specified endpoint. Data members:
# from The base point (origin when it was originally created)
# to End point when originally created
# pol Polarity of the line (see codes above)
# az Azimuth, integer in 0..9.
# flip Whether flipped end for end (codes: fl = flipped, un = not)
# rhombi Array pointing to 1 or 2 rhombi which this edge is part of
use Carp;
# Constructs an Edge. Args in order: class, from, to, $pol, az. If az is
# undef, the program will compute it. Polarity is inferred from $from and $to.
sub Edge::new {
my($this) = bless({flip => 'un'}, shift);
@{$this}{qw(from to pol az)} = @_;
$this->{az} = $this->{from}->myatan($to) unless defined($this->{az});
$this;
}
# If $a == "from" point, returns $this. If not, $this is copied in reverse
# order and a ref. to the copy is returned.
sub Edge::flip {
my($this, $a) = @_;
bless(($a->{can} == $this->{from}{can}) ? $this : {
from => $this->{to},
to => $this->{from},
pol => $this->{pol} ^ 1,
az => ($this->{az} + $N2) % $N,
flip => (($this->{flip} eq 'un') ? 'fl' : 'un'),
rhombi => $this->{rhombi},
}, ref($this));
}
# Returns the polarity of the Edge, assuming $a is the "from" point.
sub Edge::pol {
my($this, $a) = @_;
($a->{can} == $this->{from}{can}) ? $this->{pol} : ($this->{pol} ^ 1);
}
# Returns the azimuth of the Edge, assuming $a is the "from" point.
sub Edge::az {
my($this, $a) = @_;
($a->{can} == $this->{from}{can}) ? $this->{az} : ($this->{az} + $N2) % $N;
}
# Returns the opposite end of the edge from $a.
sub Edge::to {
my($this, $a) = @_;
($a->{can} == $this->{from}{can}) ? $this->{to} : $this->{from};
}
# Returns a canonical symbol for an edge, independent of the direction.
sub Edge::can {
my($this) = @_;
my($ac) = $this->{from}{can};
my($bc) = $this->{to}{can};
($ac < $bc) ? "$ac-$bc" : "$bc-$ac";
}
# Removes an edge from whatever needs removing. The edge itself is not
# destroyed, but an endpoint left without edges is destroyed.
sub Edge::remove {
my($this) = @_;
my($e, $a);
foreach $e (($this, $this->flip($this->{to}))) {
$a = $e->{from};
delete $a->{azs}{$e->{az}};
delete $a->{dest}{$e->{to}{can}};
$a->remove() if scalar(%{$a->{azs}}) <= 0;
}
}
# Produces (on STDOUT) an arbitrary path oriented on an edge. Args:
# $this Ref. to the edge
# $a Origin point (this subroutine flips the edge if needed)
# $x Coord. along the length of the edge, 0=from, 1=to.
# $y Transverse coordinate
# $move 1 means this is the 1st point in a sub-line, 0 = continuation.
# (The above 3 args repeat arbitrarily many times.)
# Returns: Nothing.
sub Edge::pspath {
my($this, $a) = splice(@_, 0, 2);
my($x0) = $opt_X + $opt_s * $a->{x};
my($y0) = $opt_Y + $opt_s * $a->{y};
my($az) = $this->az($a) * $pi5;
my($cs) = cos($az) * $opt_s;
my($si) = sin($az) * $opt_s;
my($x, $y, $mv, @path);
push(@path, 'newpath');
while (($x, $y, $mv) = splice(@_, 0, 3)) {
push(@path, sprintf("%5.1f %5.1f",
$x0 + $cs*$x - $si*$y,
$y0 + $si*$x + $cs*$y),
($mv ? 'moveto' : 'lineto'));
}
push(@path, "stroke\n");
print join(' ', @path);
}
# Produces (on STDOUT) a representation of one edge, with a marker for arrows
# or dots if $opt_d. Args: the color (if present,
# additional arg(s) will be joined with ' ' and printed as the argument
# to setcolor.)
%psarrow = (
n => [ qw(0 0 1 1 0 0) ], #Normal display, no arrow
0 => [ qw(0 0 1 1 0 0 0.1 0.0 1 0.2 0.1 0 0.3 0.0 0 0.2 -0.1 0
0.1 0.0 0) ],
1 => [ qw(0 0 1 1 0 0 0.9 0.0 1 0.8 0.1 0 0.7 0.0 0 0.8 -0.1 0
0.9 0.0 0) ],
2 => [ qw(0 0 1 1 0 0 0.55 0.1 1 0.45 0.0 0 0.55 -0.1 0) ],
3 => [ qw(0 0 1 1 0 0 0.45 0.1 1 0.55 0.0 0 0.45 -0.1 0) ],
);
sub Edge::postscript {
my($this) = shift;
print join(' ', @_, "setcolor ") if @_;
my($pol) = $opt_d ? $this->{pol} : 'n';
$this->pspath($this->{from}, @{$psarrow{$pol}});
}
# Dumps an Edge (returns a string). Argument 0 or undef -> report endpoints
# as canonical, 1 -> report in full.
sub Edge::print {
my($this, $full) = @_;
return "undef" unless ref($this) eq 'Edge';
my($btn) = ' btn ' . join(',', map {$_->print()} @{$this->{rhombi}});
if ($full == 2) {
return sprintf "%-11s az %d pol %-9s %s%s\n\tfrom %s\n\tto %s\n",
$this->can(), $this->{az}, $polvtx{$this->{pol}},
$this->{flip}, $btn,
$this->{from}->print(1), $this->{to}->print(1);
} elsif ($full == 1 ) {
return sprintf "%-11s az %d pol %-9s %s from %-5s to %-5s%s",
$this->can(), $this->{az}, $polvtx{$this->{pol}},
$this->{flip}, $this->{from}{can}, $this->{to}{can}, $btn;
} else {
return $this->{can};
}
}
# ===================================
package Rhomb;
# The Rhomb object represents a specific rhombus. Members:
# vtx List of vertices
# edges List of edges (edge[i] runs from vtx[i] to vtx[(i+1)%4])
# type Type of the rhombus (w or n)
# alttype If another choice is possible for the type, this is it.
# can Canonical coordinate of the center of the rhombus
# color Color index to use when drawing it
use Carp;
# Vertex letters:
# Wide: D/<--/C Narrow:D\-->\C
# / v \ ^
# A/*--/B A\*--\B (A is the dotted vertex)
# Vertex order. w=wide, n=narrow rhombus. Values are the polarity to be
# given to &addpts() when adding the next edge, and the digit in the key is the
# corresponding polarity when the previous edge was added. Rhombi are usually
# traversed in the positive direction, i.e. ABCD in the figures above, but
# everything is supposed to work if mirror image flipped.
# A B C D
%rhpol = qw( w1 0 w0 2 w2 3 w3 1
n1 0 n0 3 n3 2 n2 1 );
%rhangle = qw( w1 3 w0 2 w2 3 w3 2
n1 1 n0 4 n3 1 n2 4 );
# Creates a new Rhomb. Args:
# $class Class of object ('Rhomb'). This is a static member function.
# $a First vertex (A).
# $az Azimuth of the A-B edge (which must exist).
# $w Type of rhombus ('w' or 'n'). The program infers the
# orientation from the polarity of $e.
# $altw Alternate type (if any; undef if not).
# Returns: A rhombus object, completely filled in, and all vertices and
# edges are referenced as-is if existing, or created if not.
sub Rhomb::new {
my($this) = bless({ }, shift);
my($a, $az) = splice(@_, 0, 2);
my($a0) = $a;
@{$this}{qw(type alttype)} = @_;
my($w) = $this->{type};
my($e) = $a->{azs}{$az};
my($ef) = $e->flip($a);
my($b, $pol) = @{$ef}{qw(to pol)};
my($xav, $yav, $daz, $i, $j);
foreach $i (0..3) {
#Add the previous vertex and edge to the rhombus.
push(@{$this->{vtx}}, $a);
push(@{$this->{edges}}, $e);
$xav += $a->{x};
$yav += $a->{y};
#Indicate at the previous vertex what kind of tile was placed.
$daz = $rhangle{"$w$pol"};
substr($a->{tiles}, $az, $daz) = uc($w) . ($w x ($daz-1));
$j = length($a->{tiles}) - $N;
if ($j > 0) {
substr($a->{tiles}, 0, $j) = substr($a->{tiles}, $N, $j);
substr($a->{tiles}, $N, $j) = '';
}
#Tie the rhombus to the edge.
push(@{$e->{rhombi}}, $this);
last if $i >= 3; #Don't need to "create" initial edge
#Advance to the next vertex, which may or may not exist and
#which may or may not have an edge going to it from prev vtx.
$a = $b;
$az = ($az + $daz) % $N;
$pol = $rhpol{"$w$pol"};
$b = $a->endpt($az);
#If the next vertex has an edge, check if it's the right kind.
if (exists($a->{azs}{$az})) {
$e = $a->{azs}{$az};
$b = $e->to($a);
if ($pol != $e->pol($a)) {
$halt = $e;
$j = $i+1;
print STDERR "Incompatible polarity from vertex $j: ",
$a0->print(1),
"\n\tAdding ", $polvtx{$pol}, " edge ",
(qw(A-B B-C C-D D-A))[$j],
", existing edge:\n\t", $e->print(2);
}
} else {
#Create a new edge to the next vertex.
$b = $a->addpt($b, $az, $pol);
$e = $a->{dest}{$b->{can}};
}
}
$this->{can} = Point->new(0.25*$xav, 0.25*$yav)->{can};
$this;
}
# Removes a rhombus. Disconnects it from its vertex points, and removes
# edges and entire points if they are not part of a neighboring rhombus.
sub Rhomb::remove {
my($this) = @_;
my($daz, $i, $j);
my($w) = $this->{type};
my($a) = $this->{vtx}[0];
my($e, $ef, $b, $az, $pol);
foreach $e (@{$this->{edges}}) {
$ef = $e->flip($a);
($b, $az, $pol) = @{$ef}{qw(to az pol)};
#At the previous vertex, invalidate the rhombus.
$daz = $rhangle{"$w$pol"};
substr($a->{tiles}, $az, $daz) = ('.' x $daz);
$j = length($a->{tiles}) - $N;
if ($j > 0) {
substr($a->{tiles}, 0, $j) = substr($a->{tiles}, $N, $j);
substr($a->{tiles}, $N, $j) = '';
}
#Disconnect the rhombus from the edge.
splice(@{$e->{rhombi}}, ($e->{rhombi}[0]{can} != $this->{can}), 1);
#Lose the edge if there's no rhombus on the other side.
#Includes deleting the endpoint if appropriate.
$e->remove() if substr($a->{tiles}, ($az + $N - 1) % $N, 1) eq '.';
#Advance to the next vertex.
$a = $b;
}
}
# Checks if a new Rhomb will fit in the given space. Same args as new().
# $class Class of object ('Rhomb'). This is a static member function.
# $a First vertex.
# $az Azimuth of the A-B edge (which must exist).
# $w Type of rhombus ('w' or 'n'). The program infers the
# orientation from the polarity of the $az edge.
# Returns: 1 if OK.
# Bad configurations are possible that this subroutine cannot detect. For
# example, suppose you're checking a narrow rhombus that goes upward, but a
# narrow rhombus projects horizontally through that space but doesn't share
# a vertex. That won't be recognized.
sub Rhomb::check {
my($class, $a, $az, $w) = @_;
my($res); #The eventual result (undef -> OK)
my($a0) = $a;
my($e) = $a->{azs}{$az};
my($ef) = $e->flip($a);
my($b, $pol) = @{$ef}{qw(to pol)};
my($i, $j, $daz, $dots);
foreach $i (0..3) {
#At the previous vertex, is the azimuth range unoccupied
#by any tile?
$daz = $rhangle{"$w$pol"};
$dots = substr($a->{tiles}, $az, $daz);
$j = $az + $daz - $N;
$dots .= substr($a->{tiles}, 0, $j) if $j > 0;
if ($dots ne '.' x $daz) {
$res = "Rhombus overlaps another at vertex " .
substr('ABCD', $i, 1);
last;
}
#Does the new tile make a legal vertex star at this vertex?
unless ($a->check($az, $w)) {
$res = "Illegal star at vertex " . substr('ABCD', $i, 1);
last;
}
last if $i >= 3; #A-B edge is compatible by definition
#Advance to the next vertex, which may or may not exist and
#which may or may not have an edge going to it from prev vtx.
$a = $b;
$az = ($az + $daz) % $N;
$pol = $rhpol{"$w$pol"};
$b = $a->endpt($az);
$b = $Point::real{$b->{can}} if exists($Point::real{$b->{can}});
#If the next vertex has an edge, check if it's the right kind.
if (exists($a->{azs}{$az})) {
$e = $a->{azs}{$az};
$b = $e->to($a);
$j = $e->pol($a);
if ($pol != $j) {
$res = sprintf "Edge %s would be %s but existing edge is %s\n",
(qw(A-B B-C C-D D-A))[$i+1], $polvtx{$pol}, $polvtx{$j};
last;
}
}
}
$res eq '';
}
# Prints (on STDOUT) a rhombus as PostScript. Args: the color (if present,
# additional arg(s) will be joined with ' ' and printed as the argument
# to setcolor.)
sub Rhomb::postscript {
my($this) = shift;
my($b, @path);
push(@path, @_, "setcolor") if @_;
push(@path, 'newpath');
my($op) = 'moveto';
foreach $b (@{$this->{vtx}}) {
push(@path, sprintf("%5.1f %5.1f",
$opt_X + $opt_s * $b->{x},
$opt_Y + $opt_s * $b->{y}),
$op);
$op = 'lineto';
}
push(@path, "closepath", "fill\n");
print join(' ', @path);
}
# Prints the Rhomb (returns a string). Generic arg: 0 = identify, 1 = 1-line
# debug, 2 = full debug dump.
sub Rhomb::print {
my($this, $full) = @_;
my($v, $j);
my($res) = sprintf "$Point::Cfmt(%s)", $this->{can},
(defined($this->{alttype}) ? "${$this}{type}/${$this}{alttype}" :
$this->{type});
if ($full <= 0) {
# That's all you get.
} elsif ($full == 1) {
$res .= ' [';
foreach $j (0..3) {
$v = $this->{vtx}[$j];
$res .= sprintf "$Point::Cfmt/%d ", $v->{can},
$this->{edges}[$j]->pol($v);
}
substr($res, -1) = ']';
} else { # $full == 2
foreach $j (qw(vtx edges)) {
foreach $v (@{$this->{$j}}) {
$res .= "\n\t" . $v->print(1);
}
}
$res .= "\n";
}
return $res;
}
# ===================================
package Ptlist;
# The Ptlist object contains lists of points keyed by the number of azimuths
# not filled by rhombi, that is, the number of dots in $Point->{tiles}.
# Creates a new Ptlist object. No arguments beyond the class (static
# member function).
sub new {
my($class) = @_;
bless({ }, $class);
}
# Adds a Point.
# $this Ptlist to add it to
# $a Ref. to the Point to be added.
# $end If 1, $a is added at the beginning of the respective list.
# Otherwise (normally) it's added at the end.
sub add {
my($this, $a, $end) = @_;
my($j) = $a->{tiles} =~ tr/././;
$end ? unshift(@{$this->{$j}}, $a) : push(@{$this->{$j}}, $a)
unless $j <= 0;
}
# Removes and returns a point, similar to shift(). The first point is
# returned from the list with the smallest number of dots. This strategy
# is supposed to minimize the number of rhombi that have to be thrown out
# because nearly-done points end up with an illegal vertex star pattern.
sub takeoff {
my($this) = @_;
my($j, $res);
foreach $j (sort {$a <=> $b} keys %$this) {
$res = shift @{$this->{$j}};
delete $this->{$j} unless @{$this->{$j}};
last if defined($res);
}
$res;
}
# Debug printout of a Ptlist.
sub print {
my($this, $FD) = @_;
my($j);
foreach $j (sort {$a <=> $b} keys %$this) {
print $FD "$j: ", join(' ', map {$_->{can}} @{$this->{$j}}), "\n";
}
}
# ===================================
package main;
# Prints (on stdout) a PostScript path. Arguments are alternating X,Y coords
# in real scaling. There may be arbitrarily many (at least 2 pairs).
sub main::pspath { #main
my(@path) = ('newpath');
my($op) = 'moveto';
while (@_ > 0) {
push(@path,
sprintf("%5.1f %5.1f",
($opt_X + $opt_s*$_[0]),
($opt_Y + $opt_s*$_[1])),
$op);
$op = "lineto";
splice(@_, 0, 2);
}
push(@path, "stroke\n");
print join(' ', @path);
}
our $incom = Ptlist->new(); # List of points that (may) need rhombi.
our @rhlist; # List of rhombi
our %done; # Key = canonical, value = 1 if point ever was on @incom. That
# means that its surrounding polygons either have already been
# added, or will be added when we get around to it.
#Main thread: Initialization
{
my($a) = Point->new(0, 0); #Ultimate origin point
my($b) = $a->addpt($a->endpt(0), 0, 0); #First edge
$incom->add($a); #It needs polygons around it
$done{$a->{can}}++; #This point (will be) done.
}
# Generate the Penrose tiling. (Point->surround knows to
# skip points that are finished or that are outside the clip
# rectangle.)
{
my($a, $ac, $b, $bc, $c, $ct, $t, $j);
my(@stats) = qw(0 0 0); #rhombi {added, tossed}, alternates tried
my($endmsg) = "finished normally";
my($go) = 1;
INCOM: { #Loop ends when !$go
# Find the next point from among points having the fewest dots
# in {tiles}, that is, points which are the nearest to being
# finished. This supposedly minimizes wasted work when nearly
# finished points end up with an illegal configuration.
$a = $incom->takeoff();
unless (defined($a)) {
undef $go;
next;
}
#If a point was created, removed, then created again,
#use the NEW version.
$a = $Point::real{$a->{can}};
next unless (defined($a)); #Skip if already removed.
$stats[0] -= @rhlist;
$t = $a->surround(\@rhlist, $incom);
$stats[0] += @rhlist;
if (!$t) {
if (scalar(%{$a->{azs}}) <= 0) {
$a->remove();
undef $a;
}
# Ended up with illegal vertex star. Toss rhombi
# until one is found where a different type could
# have been placed.
TOSS: while (@rhlist) {
last TOSS if $halt;
my($rh) = pop @rhlist;
$b = $rh->{vtx}[0];
my($w) = $rh->{alttype};
my($az) = $rh->{edges}[0]->az($b);
$rh->remove();
$stats[1]++;
if (!defined($w)) {
#No alternative rhombus type
next;
} elsif (!exists($b->{azs}{$az})) {
# Paranoia, this hasn't been seen to happen.
printf STDERR "At $Point::Cfmt could have placed %s but lost edge %d:\n\t%s",
$b->{can}, $w, $az, $b->print(2);
} elsif (! Rhomb->check($b, $az, $w)) {
#Alternative rhombus is hemmed in by others
next;
} else {
# Insert the alternative rhombus type.
$rh = Rhomb->new($b, $az, $w);
foreach $c (@{$rh->{vtx}}) {
$incom->add($c);
}
push(@rhlist, $rh);
$stats[0]++;
$stats[2]++;
last TOSS; #Resume normal operation.
}
}
unless (@rhlist) {
$endmsg = "terminated, every rhombus failed";
undef $go;
}
$incom->add($a, 1) if defined($a);
}
if ($halt) {
$endmsg = "terminated because an error was found";
undef $go;
}
if (++$ct >= $opt_n) {
$endmsg = "terminated, exceeded max $opt_n steps";
undef $go;
}
} continue {
if (!$go || ($ct % 10000) == 0) {
$stats[3] = scalar(@rhlist);
printf STDERR
"Statistics: %d rhombi added, %d removed, %d alternates, %d kept\n",
@stats;
}
redo if $go;
}
print STDERR "Penrose tiling generation ", $endmsg, ".\n";
print STDERR " Used $ct steps of max $opt_n; give -n to increase\n";
}
# The 4-color map problem. Not exactly optimized.
our $maxcolor = 0;
{
my(@colinit); # $color 1 $color 1, To initialize %colors
my(@cstats);
my($rh, $e, $rnbr, $c);
foreach $rh (@rhlist) {
my(%colors) = @colinit;
foreach $e (@{$rh->{edges}}) {
$rnbr = $e->{rhombi}[$e->{rhombi}[0]{can} == $rh->{can}];
delete $colors{$rnbr->{color}} if exists($rnbr->{color});
}
$c = (sort keys %colors)[0];
unless (defined($c)) {
$c = $maxcolor++;
push(@colinit, $c, 1);
push(@cstats, 0);
}
$rh->{color} = $c;
$cstats[$c]++;
}
print STDERR "Color counts: @cstats ($maxcolor colors in map)\n";
}
# Print the rhombi in PostScript.
{
last unless $opt_r;
# So far I haven't needed more than 5 colors. The first and
# second colors are the most common, with the 3rd close behind.
# 4 is on less than 10% of the tiles. 5 is very rare.
my(@colors) = (
# '.5 1 1', '1 .5 1', '1 1 .5', '.5 .5 1', '.5 1 .5', '1 .5 .5', #Pastel
# '.9 .9 .9', '.8 .8 .8', '1 1 1', '.7 .7 .7', '.6 .6 .6', #Gray
# '1 1 1', '.9 .9 .7', '.85 .85 .6', '1 1 .8', '.7 .7 .5', #Yellow
#These are CIELSH hue 17 = Pantone 11-0508 etc.
# '1 1 1', '1 .85 .55', '1 .8 .4', '.92 .73 .37', '.8 .65 .33', #Desert
# 100/0/0 95/19/17 90/25/17 85/25/17 80/25/17 (CIELSH)
#These came out pretty good, but the effect breaks up the
#letters so it's hard to read. Have to try something else.
# '1 1 1', '1 .89 .66', '.92 .78 .52', '.80 .68 .45', '.69 .58 .38', #Lighter
# All the same brightness and hue, different saturations.
# CIELSH 88/*/17 where * (saturation) is 48, 36, 24, 12, 0.
# '1 .76 .30', '.93 .75 .39', '.87 .74 .50', '.79 .73 .61', '.72 .72 .72',
# Hue 17 is too orange. Let's try CIELSH 87/*/20 with * in
# 84, 42, 0, 63, 21, so the common boundaries have more
# distinction.
'1 .81 0', '.86 .76 .31', '.7 .7 .7', '.93 .78 .14', '.79 .73 .5',
);
printf STDERR "WARNING, map colored with %d colors, but only %d preset.\n",
$maxcolors, scalar(@colors) if $maxcolors > @colors;
my($rh, $c);
print "%!\n[ (DeviceRGB) ] setcolorspace\n";
foreach $rh (@rhlist) {
$c = $rh->{color};
$c = ($c < @colors) ? $colors[$c] : '1 1 1';
$rh->postscript($c);
}
}
# Print the edges in PostScript.
# Colorizing the edges was a failure; they're printed in black.
{
last unless $opt_e; #Edges can be omitted
my($k, @col, %colors); #Keys are in order of commonness
my(@colors) = qw(
0-1 1 .6 .6
0-2 .6 1 .6
1-2 .6 .6 1
0-3 .8 .8 .4
1-3 .4 .8 .8
2-3 .8 .4 .8
other .7 .4 .4
);
while (($k, @col) = splice(@colors, 0, 4)) {
$colors{$k} = $colors{reverse($k)} = join(' ', @col);
}
my($a, $ac, $b, $bc, $e, $ec, $rh);
undef %done;
if ($opt_d) {
print "1 0 0 setcolor ";
&pspath(0, 0, $Xw, 0, $Xw, $Yw, 0, $Yw, 0, 0);
}
print "0 0 0 setcolor\n";
while (($ac, $a) = each(%real)) {
while (($bc, $e) = each(%{$a->{dest}})) {
$ec = $e->can();
next if exists($done{$ec});
undef @col;
foreach $rh (@{$e->{rhombi}}) {
push(@col, $rh->{color})
}
$k = join('-', @col);
$k = 'other' unless exists($colors{$k});
#OBSOLETE $e->postscript($colors{$k});
$e->postscript();
$done{$ec}++;
}
}
if (ref($halt) eq 'Edge') {
print "100 0 0 setcolor\n";
$halt->pspath($halt->{from}, qw(
0.08 0.1 1 0.12 0.0 0 0.08 -0.1 0
0.28 0.1 1 0.32 0.0 0 0.28 -0.1 0
0.48 0.1 1 0.52 0.0 0 0.48 -0.1 0
0.68 0.1 1 0.72 0.0 0 0.68 -0.1 0
0.88 0.1 1 0.92 0.0 0 0.88 -0.1 0));
} elsif (ref($halt) eq 'Point') {
print "100 0 0 setcolor\n";
my($e) = values %{$halt->{dest}};
$e->pspath($halt, qw(0.1 0.0 1 0.0 -0.1 0 -0.1 0.0 0
0.0 0.1 0 0.1 0.0 0));
}
}
print "showpage\n";