%!PS-Adobe- % escher_grid.ps -- using the right math % Steve Witham % 2005.02.13 % derived from spiral45.ps, which was derived from polar_grid.ps % % Attempting to emulate the grid-sketch Escher did for his % picture "Print Gallery." Both the print and his grid are shown % at my page: % http://www.tiac.net/~sw/2005/05/escher_grid_2 % % I owe the deciphering of the underlying math to the paper % "The Mathematical Structure of Escher's Print Gallery," % by Hendrik Lenstra and Bart de Smit: % http://escherdroste.math.leidenuniv.nl/notices_desmit-lenstra.pdf % % Their work is illustrated with text, graphics and animations at % "Escher and the Droste Effect": % http://escherdroste.math.leidenuniv.nl/ % % ----- about this program ----- % % I count five coordinate systems being used in this program. % If you order them according to which are transformed into % which, you get this sequence of coordinate systems and the mappings % between them: % 1) "transformed coordinates". Here we're always working % relative to a right triangle whose hypotenuse goes from (0,0) to (1,0) % and whose right angle is above the x axis. In this system, the grid % lines are straight horizontal and vertical lines, and it's easy to do % the calculations to space the grid lines and clip them against the % triangle's sides. % % Sometimes we work below the x axis in the area bounded by extending % the triangle's other two sides, plus other vertical lines, but % the points (0,0) and (1,0) still provide the grid spacing reference. % % 1-2) The map between the "transformed coordinates" and the "rectangular % world" is affine, and the peak of the triangle maps to the origin % in the rectagular world. % The map from rectangular world to "transformed coordinates" is % interp_map_xform % and the function back to rectangular world is % interp_map_invx % % 2) The "rectangular world" is the straight-line world that we can imagine % that "Print Gallery" is a distorted version of. The "Droste Effect" % site mentioned above shows one of Escher's sketches in these % coordinates. In the rectangular world, the grid lines are straight. % % Lenstra and de Smit % realized that the simplest or most consistent inference from Escher's % picture is that in the rectangular world, the town contains the gallery, % which contains the print, which portrays the town...*etcetera* in an % infinite sequence of sizes, like the picture on the Droste % cocoa box which shows a woman holding a Drost cocoa box. In fact, % Escher's picture has a white blur in the center that, to be consistent, % would be filled in with a twisty version of this sequence. % % 2-3) Lenstra and de Smit worked out that if you treat 2D locations as % complex numbers and let % h = location in the rectangular back-inferred world % and % w = location in "Print Gallery"'s curvy world, % then to map between them you raise a complex number to a constant % complex power: % w = h ^ ( ( 2 pi i ) / ( 2 pi i + ln 256 ) ) % and the inverse mapping is raising to the inverse power: % h = w ^ ( ( 2 pi i + ln 256 ) / ( 2 pi i ) ) % % There is one point that both functions map to itself: the origin % ( 0, 0 ). This is in the center of Escher's picture. % % In this program these mappings are done by the functions % Escher_map and Escher_inv_map, respectively. % % 3) There is a problem that these "functions" are really relationships % that are not one-to-one. Because the procedures to calculate them % return single results, discontinuities show up in the mappings they % do. To avoid these, the program stays within an area where the % mapping is continuous both ways. This safe area is a right triangle % in the rectangular world, with the right angle at the origin. % Two of its corners are found by mapping backwards from important % points in Escher's grid. Here are the corners in both systems: % rectangular world Escher's grid % ( 0, 0 ) ( 0, 0 ) % Escher_inv_map( ( 1, -1 ) ) ( 1, -1 ) dot he labeled "A" % Escher_inv_map( ( 1, +1 ) ) ( 1, +1 ) dot he labeled "D" % % The "safe area" of the curvy world of "Print Gallery" is a spiral % quarter-slice of the plane. % % 3-4) The procedure "four_triangles" draws this spiral slice of the grid % four times, rotating the picture -90, 0, 90 and 180 degrees. % % 4) Welcome to the curvy world of "Print Gallery". % % 4-5) Coordinates in the curvy world are multiplied by about 2.5 inches % to get coordinates on the rendered page. % In the program this is the function "scalept". % % 3-5) The program has a function called "Escher_map_and_scale" which is % defined as { Escher_map scalept }. % % 5) The rendered page has (0,0) in the center. % /inch { 72 mul } def /swap { exch } def % This works best if set to the printer resolution % (which is settable in my printer). /resolution 1 600 div inch def /minwidth resolution 6 mul def false setstrokeadjust % ----- Line "weight" (thickness & darkness) ---- /setlineweight { % "thickness" -- dup dup minwidth lt { pop minwidth } { % else resolution div ceiling resolution mul } ifelse dup setlinewidth div 1 exch sub setgray } def % Old weight-adjusting formula: /fadeweight minwidth 1 div def /sqrtfade minwidth fadeweight div sqrt def /setlineweight_old { % "thickness" -- dup fadeweight gt { % Above this weight thicken like sqrt: minwidth div sqrt minwidth mul } { % else Below it, fade like x^2: fadeweight div dup mul fadeweight mul } ifelse dup dup minwidth lt { pop minwidth } { % else resolution div ceiling resolution mul } ifelse dup setlinewidth div 1 exch sub setgray } def % These adjust the weight based on the magnification (abs value % of the derivative) of a mapping function that's assumed to be % conformal. /adjust_lineweight_setpt0 { % x0 y0 x0' y0' -- /adjust_lineweight_y0' exch def /adjust_lineweight_x0' exch def /adjust_lineweight_y0 exch def /adjust_lineweight_x0 exch def } def /adjust_lineweight_moveto { % x0 y0 -- 2 copy adjust_lineweight_mapfn 2 copy moveto adjust_lineweight_setpt0 } def /adjust_lineweight_setpt1_getmag { % x1 y1 -- magnification 2 copy /adjust_lineweight_y1 exch def /adjust_lineweight_x1 exch def adjust_lineweight_mapfn /adjust_lineweight_y1' exch def /adjust_lineweight_x1' exch def adjust_lineweight_x1' adjust_lineweight_y1' adjust_lineweight_x0' adjust_lineweight_y0' c_sub c_abs adjust_lineweight_x1 adjust_lineweight_y1 adjust_lineweight_x0 adjust_lineweight_y0 c_sub c_abs div } def % adjust_lineweight_calibrate sets the mapping function and % calibrates the weight for near two nearby points. % /adjust_lineweight_calibrate { % weight x1 y1 x0 y0 mapfn -- /adjust_lineweight_mapfn exch def adjust_lineweight_moveto newpath adjust_lineweight_setpt1_getmag div /adjust_lineweight_weight exch def /adjust_lineweight_ref_r' adjust_lineweight_x0' adjust_lineweight_y0' c_abs def } def % NOTE: THIS STROKES, CLEARS PATH, THEN DOES MOVETO. /adjust_lineweight_lineto { % x1 y1 -- final_weight adjust_lineweight_setpt1_getmag adjust_lineweight_weight mul % Up to this point the weight must be proportional to the magnification % so a continuation of the same line has the same weight as it crosses % between two of the quarter slices. % But now that that's true we can do what we want to the weight. % % The total weight of four lighter lines between two heavier lines should % be less that that of one of the thick lines, not just to make the thick % lines stand out more, but because otherwise, given the infinite succession % of resolutions, the total weight would be infinite. So, let's use a % power law. With 3/2 power, 4 * ( (1/4)^(3/2) ) = 1/2. % 3 2 div exp % This is adjusted for my printer rather than the screen (ghostview): 1.25 exp % 1.5 exp % for screen % % Now do some stuff based on how close to the center: /adjust_lineweight_r' adjust_lineweight_x1' adjust_lineweight_y1' c_abs adjust_lineweight_ref_r' div dup 1 lt { .35 exp } { pop 1 } ifelse def adjust_lineweight_r' div dup adjust_lineweight_r' gt { pop adjust_lineweight_r' } if .02 inch mul % "Sudden death" rule: % /foo .005 def % dup foo lt { foo div .5 sub 2 mul foo mul } if % dup foo lt { foo div .25 sub 1.5 mul foo mul } if dup 0 gt { dup setlineweight adjust_lineweight_x1' adjust_lineweight_y1' lineto stroke } if adjust_lineweight_x1 adjust_lineweight_y1 adjust_lineweight_x1' adjust_lineweight_y1' 2 copy moveto adjust_lineweight_setpt0 } def % ------ ----- /circle { % x y r -- 0 360 arc stroke } def /circle_through_point { % center_x center_y pt_x pt_y -- 4 copy pop pop % cx cy cx cy px py c_sub c_abs circle } def /dot { % x y r -- 0 360 arc fill } def 8.5 inch 2 div 11 inch 2 div translate /center { % -- x y 0 0 } def /tan { dup sin exch cos div } def /pow { % base exponent -- real exp } def /M_PI 3.1415926 def /M_E 2.7182818 def /deg_to_rad { 180 div M_PI mul } def /rad_to_deg { 180 mul M_PI div } def % ----- Complex Math ----- /c_add { % a b c d -- a+c b+d 3 -1 roll add 3 1 roll add exch } def /c_sub { % a b c d -- a-c b-d neg 3 -1 roll add 3 1 roll sub exch } def /c_neg { neg exch neg exch } def /c_mul { % a b c d -- (ac-bd) (ad+bc) 4 copy exch 4 -1 roll mul 3 1 roll mul sub 5 1 roll 4 -1 roll mul 3 1 roll mul add } def /c_div { % num_real num_imag denom_real denom_imag -- quot_real quot_imag % No divide-by-zero test. c_inv c_mul % Gee, that's something of an anticlimax. } def /c_inv { % real imag -- real' imag' 2 copy c_abs^2 dup % real imag r^2 r^2 4 -1 roll swap div % imag r^2 real/r^2 3 1 roll div neg % real/r^2 -imag/r^2 } def /c_abs^2 { % real imag -- magnitude^2 dup mul exch dup mul add } def /c_abs { % real imag -- magnitude c_abs^2 sqrt } def /polar { % real imag -- radius angle_in_radians /polar_y swap def /polar_x swap def polar_x polar_y c_abs dup 0 eq { 0 } % else { polar_y polar_x atan dup 180 gt { 360 sub } if deg_to_rad } ifelse } def /c_ln { % real imag -- real imag % no zero test polar swap ln swap } def /cis { % radius angle -- real imag /cis_angle exch def dup cis_angle rad_to_deg cos mul swap cis_angle rad_to_deg sin mul } def /c_pow { % base_real base_imag power_real power_imag -- real imag /c_pow_d exch def % power is c + i d /c_pow_c exch def c_ln % base is a cis b /c_pow_b exch def /c_pow_ln_a exch def M_E c_pow_c c_pow_ln_a mul c_pow_b c_pow_d mul sub pow dup /c_pow_E swap def c_pow_d c_pow_ln_a mul c_pow_b c_pow_c mul add dup /c_pow_PHI swap def cis } def % ----- The power constant for the Escher map ----- % It's ( 2 pi i ) / ( 2 pi i + ln 256 ) % % = ( 4 pi^2 + 2 pi ( ln 256 ) i ) / ( 4 pi^2 + ( ln 256 )^2 ) /ln256 256 ln def /four_pi^2 M_PI dup mul 4 mul def /Escher_pwr_D four_pi^2 ln256 dup mul add def /Escher_pwr_real four_pi^2 Escher_pwr_D div def /Escher_pwr_imag 2 M_PI mul ln256 mul Escher_pwr_D div def /Escher_pwr { % -- real imag Escher_pwr_real Escher_pwr_imag } def % ----- Interpolation ----- /interp { % x0 y0 x1 y1 fraction -- x y /interp_fraction swap def /interp_y1 swap def /interp_x1 swap def 1 interp_fraction sub mul interp_y1 interp_fraction mul add swap 1 interp_fraction sub mul interp_x1 interp_fraction mul add swap } def /interp_lineto { % x0 y0 x1 y1 n_steps -- % Assumes you already did moveto or lineto x0 y0. /interp_lineto_n_steps exch def 4 copy 1 1 interp_lineto_n_steps { interp_lineto_n_steps div interp lineto 4 copy } for pop pop pop pop } def /interp_map_lineto { % x0 y0 x1 y1 n_steps mapfn -- % Assumes you already did moveto or lineto x0 y0 mapfn. /interp_map_lineto_mapfn exch def /interp_map_lineto_n_steps exch def 4 copy 1 1 interp_map_lineto_n_steps { interp_map_lineto_n_steps div interp interp_map_lineto_mapfn lineto 4 copy } for pop pop pop pop } def /interp_map_adjust_lineto { % x0 y0 x1 y1 n_steps mapfn -- % Assumes you already did adjust_lineweight_moveto or _lineto x0 y0 mapfn. /interp_map_lineto_mapfn exch def /interp_map_lineto_n_steps exch def 4 copy 1 1 interp_map_lineto_n_steps { interp_map_lineto_n_steps div interp % interp_map_lineto_mapfn adjust_lineweight_lineto pop % Pop the final line weight for now. 4 copy } for pop pop pop pop } def /interp_map_triangle { % x0 y0 x1 y1 n_steps mapfn -- /interp_map_triangle_mapfn exch def /interp_map_triangle_n_steps exch def /interp_map_triangle_y1 exch def /interp_map_triangle_x1 exch def /interp_map_triangle_y0 exch def /interp_map_triangle_x0 exch def % Say P is ( x0, y0 ) and Q is ( x1, y1 ). interp_map_triangle_x0 interp_map_triangle_y0 interp_map_triangle_x1 interp_map_triangle_y1 interp_map_xform_setup % --- Lines parallel to line PQ: --- /interp_map_triangle_ystep 1 interp_map_triangle_n_steps div def % "Quantize" the min y (find the next even step): interp_map_triangle_miny interp_map_triangle_ystep div ceiling interp_map_triangle_ystep mul % for-loop starting value interp_map_triangle_ystep % step % Don't go higher than the origin. interp_map_triangle_maxy dup interp_map_xform_b gt { pop interp_map_xform_b } if % limit { % for /interp_map_triangle_y exch def interp_map_triangle_y interp_map_xform_b div interp_map_xform_a mul interp_map_triangle_y interp_map_invx 2 copy adjust_lineweight_moveto 1 1 interp_map_xform_a sub interp_map_triangle_y mul interp_map_xform_b div sub interp_map_triangle_y interp_map_invx % interp_map_triangle_n_steps 4 mul 48 /interp_map_triangle_mapfn load interp_map_adjust_lineto % stroke } for % --- Lines perpendicular to line PQ: --- % Naive min and max x at slice at miny: /interp_map_triangle_minx interp_map_triangle_miny interp_map_xform_a mul interp_map_xform_b div def /interp_map_triangle_maxx interp_map_triangle_miny 1 interp_map_xform_a sub mul interp_map_xform_b div neg 1 add def % Now convert these to step numbers, but % avoid drawing zero-height lines at the corners, by % pushing them inwards... interp_map_triangle_minx interp_map_triangle_n_steps mul .5 add ceiling % start 1 % step interp_map_triangle_maxx interp_map_triangle_n_steps mul .5 sub floor % limit { /interp_map_triangle_x exch interp_map_triangle_n_steps div def interp_map_triangle_x interp_map_triangle_miny interp_map_invx 2 copy adjust_lineweight_moveto interp_map_triangle_x dup dup interp_map_xform_a lt { % left of origin interp_map_xform_b mul interp_map_xform_a div } { % else right of origin neg 1 add interp_map_xform_b mul interp_map_xform_a neg 1 add div } ifelse dup interp_map_triangle_maxy gt { pop interp_map_triangle_maxy } if interp_map_invx 48 /interp_map_triangle_mapfn load interp_map_adjust_lineto % stroke } for } def % interp_map_xform is an affine transform that maps % ( x0, y0 ) => ( 0, 0 ) % ( x1, y1 ) => ( 1, 0 ) % ( 0, 0 ) => ( interp_map_xform_a, interp_map_xform_b > 0 ) /interp_map_xform { % x y -- x' y' /interp_map_xform_y exch def /interp_map_xform_x exch def interp_map_xform_x interp_map_xform_c mul interp_map_xform_x interp_map_xform_s mul exch interp_map_xform_y interp_map_xform_s mul sub interp_map_xform_a add exch interp_map_xform_y interp_map_xform_c mul add interp_map_xform_b add } def /interp_map_xform_zero { % -- % Set interp_map_xform_a and _b so ( _x0, _y0 ) => ( 0, 0 ). /interp_map_xform_a 0 def /interp_map_xform_b 0 def interp_map_xform_x0 interp_map_xform_y0 interp_map_xform /interp_map_xform_b exch neg def /interp_map_xform_a exch neg def } def % interp_map_invx is the inverse of interp_map_xform. /interp_map_invx { % x' y' -- x y /interp_map_invx_y exch def /interp_map_invx_x exch def interp_map_invx_x interp_map_invx_c mul interp_map_invx_x interp_map_invx_s mul exch interp_map_invx_y interp_map_invx_s mul sub interp_map_invx_a add exch interp_map_invx_y interp_map_invx_c mul add interp_map_invx_b add } def /interp_map_invx_setup { % -- interp_map_xform_c interp_map_xform_s c_inv /interp_map_invx_s exch def /interp_map_invx_c exch def % Set interp_map_invx_a and _b so ( 0, 0 ) => ( xform_x0, _y0 ): /interp_map_invx_a interp_map_xform_x0 def /interp_map_invx_b interp_map_xform_y0 def } def /interp_map_xform_setup { % x0 y0 x1 y1 -- % Stores xform params in globals ..._c, _s, _a, _b. /interp_map_xform_y1 exch def /interp_map_xform_x1 exch def /interp_map_xform_y0 exch def /interp_map_xform_x0 exch def /interp_map_xform_dx interp_map_xform_x1 interp_map_xform_x0 sub def /interp_map_xform_dy interp_map_xform_y1 interp_map_xform_y0 sub def /interp_map_xform_dx2 interp_map_xform_dx dup mul def /interp_map_xform_dy2 interp_map_xform_dy dup mul def interp_map_xform_dx2 interp_map_xform_dy2 ge { /interp_map_xform_c interp_map_xform_dx interp_map_xform_dx2 interp_map_xform_dy2 add div def /interp_map_xform_s interp_map_xform_c interp_map_xform_dy mul neg interp_map_xform_dx div def } if interp_map_xform_dx2 interp_map_xform_dy2 lt { /interp_map_xform_s interp_map_xform_dy neg interp_map_xform_dx2 interp_map_xform_dy2 add div def /interp_map_xform_c interp_map_xform_s interp_map_xform_dx mul neg interp_map_xform_dy div def } if interp_map_xform_zero % Flip so ( 0, 0 ) => above the x axis: 0 0 interp_map_xform exch pop 0 lt { /interp_map_xform_s /interp_map_xform_s neg def interp_map_xform_zero } if interp_map_invx_setup } def % ***** Because of how setlineweight works, I can't change the global scale! ***** /the_scale 2.5 inch def /scalept { the_scale mul swap the_scale mul swap } def /edge 4.125 inch def /-edge edge neg def /draw_axes { % -- .0025 inch setlineweight -edge 0 moveto edge 0 lineto stroke 0 -edge moveto 0 edge lineto stroke -edge -edge moveto edge -edge lineto edge edge lineto -edge edge lineto -edge -edge lineto closepath stroke } def % draw_axes /Escher_map { % real imag -- x y Escher_pwr c_pow } def /Escher_map_and_scale { % real imag -- x y Escher_map scalept } def /Escher_inv_map { % x y -- real imag Escher_pwr c_inv c_pow } def /four_triangles { % n_steps -- << exch /n_steps exch % def >> begin -90 90 180 { % for gsave rotate 1 -1 Escher_inv_map 1 1 Escher_inv_map n_steps /Escher_map_and_scale load interp_map_triangle grestore } for end } def % /Times-Roman findfont 24 scalefont setfont % .35 setgray center moveto (hello?) show /the_scale 2.5 inch def % ----- Interpolation on grid lines determines radii of circles ----- % grid_intersection locates a point on Escher's grid at the intersection of % a "vertical" line numbered from 0 at the bottom-right dot (point A) % to 6 at the bottom-left dot (point B) % but only lines that pass to the right of the center, % and a "horizontal" line numbered 0 at the top-right dot (point D) % to 6 at the bottom-right dot (point A) % The result point is in page coordinates. /grid_intersection { % leftness*6 downness*6 -- x y 6 div exch 6 div dup 16 div neg 1 add -1 1 Escher_inv_map 1 1 Escher_inv_map 5 -1 roll interp % downness leftness x0 y0 3 -1 roll 1 -1 Escher_inv_map -1 -1 Escher_inv_map 5 -1 roll interp % downness x0 y0 x1 y1 5 -1 roll interp Escher_map_and_scale } def /set_corner_weight { % weight -- 1 1 Escher_inv_map 1 47 48 div Escher_inv_map /Escher_map_and_scale load adjust_lineweight_calibrate } def % ========== NOW ACTUALLY DRAWING ========== % MINIMUM AND MAXIMUM Y % The idea here is that we are only drawing stuff in a horizontal band, % i.e. between two vertical lines, in the transformed space. This % can extend below and/or above the base of the triangle (i.e. % miny and maxy can be < 0 or > 0). /interp_map_triangle_miny 5 24 div def /interp_map_triangle_maxy 1 3 div def 1 256 div set_corner_weight 1536 four_triangles 1 64 div set_corner_weight 384 four_triangles % ----- smallish circle around the middle ----- 5.5 5.5 grid_intersection 2 copy % small circle point .25 setgray .025 inch dot c_abs .005 inch setlineweight center 3 -1 roll circle 1 16 div set_corner_weight 96 four_triangles 1 4 div set_corner_weight 24 four_triangles /interp_map_triangle_miny 1 6 div def % /interp_map_triangle_maxy 5.01 24 div def /interp_map_triangle_maxy 5 24 div def 1 64 div set_corner_weight 384 four_triangles 1 16 div set_corner_weight 96 four_triangles 1 4 div set_corner_weight 24 four_triangles 1 1 div set_corner_weight 6 four_triangles % ----- medium-sized circle around the middle ----- 4 5 4.5 16 div add grid_intersection 2 copy % big circle point .25 setgray .025 inch dot c_abs .005 inch setlineweight center 3 -1 roll circle /interp_map_triangle_miny -1 2 div def /interp_map_triangle_maxy 1 6 div def 1 16 div set_corner_weight 96 four_triangles 1 4 div set_corner_weight 24 four_triangles 1 1 div set_corner_weight 6 four_triangles % ----- little circles around points A, B, C & D ----- /little_circle_r 0 6 grid_intersection 1 12 div 6 grid_intersection c_sub c_abs def /little_circle { little_circle_r circle } def .01 inch setlineweight -1 the_scale mul -1 the_scale mul little_circle 1 the_scale mul -1 the_scale mul little_circle 1 the_scale mul 1 the_scale mul little_circle -1 the_scale mul 1 the_scale mul little_circle showpage