ping luser- -droog: SVG or PS project for you?

[tlvp]

I'm not sure whether you're up to taking on others' tasks, but here's a
geometrical one my own PS and/or SVG capabilities just aren't up to at all.

I'd like to find code for getting 2D-projected views of the 3D object
whose x- and y- and z-coordinates (after suitable scaling, say) satisfy:

: Max [ x^2+y^2, x^2+z^2, y^2+z^2 ] = 4 (-2 <= x&y&z <= +2) .

You can visualize this as the result of starting with a 4-inch length of
two-inch radius wooden doweling (that helps assure, say, x^2+y^2 <= 4),
spinning it against a belt sander with spinning axis perpendicular to --
and bisecting -- the central axis from the one disc-shaped end to the other
until you've removed enough material to assure x^2+z^2 <= 4 (if it were a
stubby D cell battery, it would be spinning end over end, grinding the +
and - poles down until they're well rounded, sort of), and then using a
third axis perpendicular to the two axes already mentioned as a spinning
axis for grinding down the last bits (still at the + and - poles), to
assure y^2+z^2 <= 4.

[If you have Mathematica, tell it you want it to give you

: RegionPlot3D[Max[x^2+y^2,x^2+z^2,y^2+z^2]<=4,{x,-2,2},{y,-2,2},{z,-2,2}]] .

(Or try with the last 6 instances of the values "2", in the braces,
replaced by "2.01".)
With Maple, I suspect some very similar formulation will work as well. But
I have neither of them, and Wolfram|Alpha balks at giving the help I want.]

Anyway, viewed "head-on" from the perspective of any of those three axes,
the resulting object has circular silhouette. Viewed from virtually any
other perspective, not properly circular at all.

Any way to come up with a SVG animation showing views of this beast as
one's viewing angle varies from elevations of 0 up through 45 degrees,
and/or one's lateral angle varies between -45 and +45 degrees (leftward
and/or rightward), would be welcome. Or even just a short slide-show of PS
views :-) . (I think of your snowfkake iterations as a "slide-show" model.)

TIA, if possible; and no hard feelings if not. Cheers, -- tlvp
--
Avant de repondre, jeter la poubelle, SVP.

[luser- -droog]

On Tuesday, July 17, 2012 3:48:33 AM UTC-5, tlvp wrote:
> I&#39;m not sure whether you&#39;re up to taking on others&#39; tasks, but here&#39;s a
> geometrical one my own PS and/or SVG capabilities just aren&#39;t up to at all.
>
> I&#39;d like to find code for getting 2D-projected views of the 3D object

> whose x- and y- and z-coordinates (after suitable scaling, say) satisfy:
>

> : Max [ x^2+y^2, x^2+z^2, y^2+z^2 ] = 4 (-2 &lt;= x&amp;y&amp;z &lt;= +2) .

>
> You can visualize this as the result of starting with a 4-inch length of

> two-inch radius wooden doweling (that helps assure, say, x^2+y^2 &lt;= 4),

> spinning it against a belt sander with spinning axis perpendicular to --
> and bisecting -- the central axis from the one disc-shaped end to the other

> until you&#39;ve removed enough material to assure x^2+z^2 &lt;= 4 (if it were a

> stubby D cell battery, it would be spinning end over end, grinding the +

> and - poles down until they&#39;re well rounded, sort of), and then using a

> third axis perpendicular to the two axes already mentioned as a spinning
> axis for grinding down the last bits (still at the + and - poles), to

> assure y^2+z^2 &lt;= 4.

>
> [If you have Mathematica, tell it you want it to give you
>

> : RegionPlot3D[Max[x^2+y^2,x^2+z^2,y^2+z^2]&lt;=4,{x,-2,2},{y,-2,2},{z,-2,2}]] .
>
> (Or try with the last 6 instances of the values &quot;2&quot;, in the braces,
> replaced by &quot;2.01&quot;.)

> With Maple, I suspect some very similar formulation will work as well. But
> I have neither of them, and Wolfram|Alpha balks at giving the help I want.]
>

> Anyway, viewed &quot;head-on&quot; from the perspective of any of those three axes,

> the resulting object has circular silhouette. Viewed from virtually any
> other perspective, not properly circular at all.
>
> Any way to come up with a SVG animation showing views of this beast as

> one&#39;s viewing angle varies from elevations of 0 up through 45 degrees,
> and/or one&#39;s lateral angle varies between -45 and +45 degrees (leftward

> and/or rightward), would be welcome. Or even just a short slide-show of PS

> views :-) . (I think of your snowfkake iterations as a &quot;slide-show&quot; model.)

>
> TIA, if possible; and no hard feelings if not. Cheers, -- tlvp
> --
> Avant de repondre, jeter la poubelle, SVP.

I'll give it a go. Might have to reread the second half
of Mathematical Illustrations. It would be a bridge to
one of my ideas: an animated stereogram of a rotating
tesseract.

[tlvp]

On Tue, 17 Jul 2012 08:16:41 -0700 (PDT), luser- -droog wrote,
and I am astonished at all the substitutions of HTML entity names for the
actual LessThan or Apostrophe or Ampersand or QuoteMark characters I used:

> Path: n102.xanadu-bbs.net!optima5.xanadu-bbs.net!optima11.xanadu-bbs.net!optima2.xanadu-bbs.net!v102.xanadu-bbs.net!xanadu-bbs.net!news.glorb.com!npeer03.iad.highwinds-media.com!news.highwinds-media.com!feed-me.highwinds-media.com!postnews.google.com!glegroupsg2000goo.googlegroups.com!not-for-mail
> From: luser- -droog <mij...@yahoo.com>
> Newsgroups: comp.lang.postscript
> Subject: Re: ping luser- -droog: SVG or PS project for you?
> Date: Tue, 17 Jul 2012 08:16:41 -0700 (PDT)
> Organization: http://groups.google.com
> Lines: 47
> Message-ID: <aff1ff46-4811-4bfe...@googlegroups.com>
> References: <1ca6lvkhx25dt.14pkcuyecuct2$.d...@40tude.net>
> NNTP-Posting-Host: 99.191.124.223
> Mime-Version: 1.0
> Content-Type: text/plain; charset=ISO-8859-1
> X-Trace: posting.google.com 1342538202 8793 127.0.0.1 (17 Jul 2012 15:16:42 GMT)
> X-Complaints-To: groups...@google.com
> NNTP-Posting-Date: Tue, 17 Jul 2012 15:16:42 +0000 (UTC)
> Cc: mPiOsUcB...@att.net
> In-Reply-To: <1ca6lvkhx25dt.14pkcuyecuct2$.d...@40tude.net>
> Complaints-To: groups...@google.com
> Injection-Info: glegroupsg2000goo.googlegroups.com; posting-host=99.191.124.223;
> posting-account=G1KGwgkAAAAyw4z0LxHH0fja6wAbo7Cz
> User-Agent: G2/1.0
> X-Received-Bytes: 3613
> Xref: n102.xanadu-bbs.net comp.lang.postscript:2330

Think it's Google Groups doing that? Or your browser? Or other? (OtOH, your
own Apostrophe, in your "I'll give it a go" phrase, is just fine :-) .

Anyway, I'm heartened to see that you'll "give it a go" -- I won't hold my
breath for it, but I will be interested to see what you come up with.

Cheers, and TIA, -- tlvp

[luser- -droog]

On Tuesday, July 17, 2012 3:48:33 AM UTC-5, tlvp wrote:

Alright, version 1 has some obvious flaws.
But it looks really strange.
I think I should solve for z instead of sampling.
And there's no projection here, just dropping the z
(interpreting it as a grey value).

I borrowed heavily from my mandelbrot set program
in the use of `image`.

%!

/max {
2 copy lt { exch } if pop
} def

%x y z f bool
/f {
dup mul 3 1 roll %z^2 x y
dup mul 3 1 roll %y^2 z^2 x
dup mul 3 1 roll %x^2 y^2 z^2
2 index 2 index add %x^2 y^2 z^2 x^2+y^2
3 index 2 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2
3 index 3 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2 y^2+z^2
max %x^2 y^2 z^2 x^2+y^2 max(x^2+z^2,y^2+z^2)
max %x^2 y^2 z^2 max(x^2+y^2, max(x^2+z^2,y^2+z^2))
4 1 roll pop pop pop %max(...)
4 eq
} def

/tlvp {
/res exch def %resolution
/str res string def %1-row buffer
/++ { dup load 1 add store } bind def %increment proc
/yc 0 def %y cursor
res res 8 [ res 0 0 res neg 0 res ]
{ %image
0 1 res 1 sub { %for each byte of the string, incr x
/xc exch def %x cursor
/x0 xc res div 4 mul 2 sub def %
/y0 yc res div 4 mul 2 sub def %

0 1 res 1 sub { %check all z values in range
/zc exch def
/z0 zc res div 4 mul 2 sub def
x0 y0 z0 f {
str xc
z0 2 add 4 div 255 mul cvi
put
} if
} for
} for
/yc ++ %incr y for next row
str %yield string buffer to `image`
} image
} def

612 612 scale
100 tlvp

showpage

[luser- -droog]

On Saturday, August 18, 2012 1:40:51 AM UTC-5, luser- -droog wrote:
> On Tuesday, July 17, 2012 3:48:33 AM UTC-5, tlvp wrote:
>
> > I'm not sure whether you're up to taking on others' tasks, but here's a
>
> >
>
> > geometrical one my own PS and/or SVG capabilities just aren't up to at all.
>
> >
>
> >
>
> >
>
> > I'd like to find code for getting 2D-projected views of the 3D object
>
> >
>
> > whose x- and y- and z-coordinates (after suitable scaling, say) satisfy:
>
> >
>
> >
>
> >
>
> > : Max [ x^2+y^2, x^2+z^2, y^2+z^2 ] = 4 (-2 <= x&y&z <= +2) .
>
> >
>
> >

[snip]

>
>
>
> Alright, version 1 has some obvious flaws.
>
> But it looks really strange.
>
> I think I should solve for z instead of sampling.
>
> And there's no projection here, just dropping the z
>
> (interpreting it as a grey value).
>
>
>
> I borrowed heavily from my mandelbrot set program
>
> in the use of `image`.
>

Forgot to clear the bytes, so it was smearing.
Now it uses 128 for a mid-grey background.
And a little fuzzy on the notion of equality
starts to bring it to life, I think.

%!

/oldeq /eq load def
/eq {
ceiling exch ceiling oldeq
} def

/max {
2 copy lt { exch } if pop
} def

%x y z f bool
/f {
dup mul 3 1 roll %z^2 x y
dup mul 3 1 roll %y^2 z^2 x
dup mul 3 1 roll %x^2 y^2 z^2
2 index 2 index add %x^2 y^2 z^2 x^2+y^2
3 index 2 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2
3 index 3 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2 y^2+z^2
max %x^2 y^2 z^2 x^2+y^2 max(x^2+z^2,y^2+z^2)
max %x^2 y^2 z^2 max(x^2+y^2, max(x^2+z^2,y^2+z^2))
4 1 roll pop pop pop %max(...)
4 eq
} def

/tlvp {

/dep exch def %depth

/res exch def %resolution
/str res string def %1-row buffer
/++ { dup load 1 add store } bind def %increment proc
/yc 0 def %y cursor
res res 8 [ res 0 0 res neg 0 res ]
{ %image
0 1 res 1 sub { %for each byte of the string, incr x
/xc exch def %x cursor
/x0 xc res div 4 mul 2 sub def %
/y0 yc res div 4 mul 2 sub def %

str xc 128 put

0 1 dep 1 sub { %check all z values in range
/zc exch def
/z0 zc dep div 4 mul 2 sub def

x0 y0 z0 f {
str xc
z0 2 add 4 div 255 mul cvi
put
} if
} for
} for
/yc ++ %incr y for next row
str %yield string buffer to `image`
} image
} def

0 150 translate
612 612 scale
100 100 tlvp %x&y z don't push it with these, it's O(n^3)!

showpage

[tlvp]

You've got it :-) ! ... But the "project z out of the picture" strategy
loses the sort of perspective that gives a better feel for the object.

It's a so-called "Steinmetz Solid" (q.v. -- midway down the first page of
Google hits, as they come up for me, are some neat small renditions of it).

[luser- -droog]

On Saturday, August 18, 2012 9:44:21 PM UTC-5, tlvp wrote:
> On Sat, 18 Aug 2012 00:14:56 -0700 (PDT), luser- -droog wrote:
>
>
> > Forgot to clear the bytes, so it was smearing.
>
> > Now it uses 128 for a mid-grey background.
>
> > And a little fuzzy on the notion of equality
>
> > starts to bring it to life, I think.
>
> >
>

[snip]

>
>
> You've got it :-) ! ... But the "project z out of the picture" strategy
>
> loses the sort of perspective that gives a better feel for the object.
>
>
>
> It's a so-called "Steinmetz Solid" (q.v. -- midway down the first page of
>
> Google hits, as they come up for me, are some neat small renditions of it).
>

I've made a little progress here. It still generates the points
with the O(n^3) brute force algorithm, but it saves the points
to a data file called "stein.pts" so subsequent pages are much
less agonizing.

Rotations and projections are applied to each point, and then
it draws a little circle to represent the projected point.

%!
% Point-field sampling
% with data caching (in a file),
% point-wise axial rotations,
% and perspective projection.

/fuzz 10000 def %the "grain" of eq

/oldeq /eq load def
/eq {

fuzz mul round exch fuzz mul round oldeq

} def

/max {
2 copy lt { exch } if pop
} def

%x y z f bool
/f {
dup mul 3 1 roll %z^2 x y
dup mul 3 1 roll %y^2 z^2 x
dup mul 3 1 roll %x^2 y^2 z^2
2 index 2 index add %x^2 y^2 z^2 x^2+y^2
3 index 2 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2
3 index 3 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2 y^2+z^2
max %x^2 y^2 z^2 x^2+y^2 max(x^2+z^2,y^2+z^2)
max %x^2 y^2 z^2 max(x^2+y^2, max(x^2+z^2,y^2+z^2))
4 1 roll pop pop pop %max(...)
4 eq
} def

/filename (stein.pts) def

/low -2.2 def
/hi 2.2 def

%generate data by brute force, cache in file, plot
/pointfieldtocache {
/res exch def
/dt 1 res div def
/fuzz res .5 mul def
%fuzz affects the "closeness" of
%the equality test. a lower fuzz will allow more
%values to be equal.
% 'res' gives thin lines
% 'res .5 mul' gives wider "ribbons"

/outfile filename (w) file def
/outbuf 128 string def
low dt hi {
low dt hi {
low dt hi { % xW yW zW "world" coords
3 copy f {
3 copy 3 -1 1 { -1 roll
outbuf cvs outfile exch writestring
outfile (\n) writestring
} for %dump points to file
3 copy project
%2 copy exch = =
%2 copy transform exch = = ()=
2 copy 2 copy moveto dt .5 mul 0 360 arc moveto
fill
/flushpage where {pop flushpage} if
} if
pop
} for
pop
} for
pop
} for
outfile closefile
} def

%plot cached data from file
/pointfieldfromcache {
/res exch def
/dt 1 res div def
/infile filename (r) file def
/it 1 def
%cvx exec
%count 3 idiv {

{
{
/it it 1 add def
infile token not {stop} if % bail-out
infile token not {stop} if % on any datafile issues
infile token not {stop} if
%3 copy
project
2 copy 2 copy moveto dt .5 mul 0 360 arc moveto
it res mod 0 oldeq {
fill /flushpage where {pop flushpage} if
} if
} loop
} stopped pop
%} repeat
} def

/pointfield { %check if there's a readable data file
{ filename (r) file closefile } stopped not
//pointfieldfromcache %yes! read it.
{ pop pop pointfieldtocache } ifelse %no! make one.
} def

% x y z ang -> x y' z'
/rotx {
/theta exch def
/z exch def
/y exch def
y theta cos mul
z theta sin mul sub
y theta sin mul
z theta cos mul add
} def

% x y z ang -> x' y z'
/roty {
/theta exch def
/z exch def
/y exch def
/x exch def
x theta cos mul
z theta sin mul add
y
x theta sin mul neg
z theta cos mul add
} def

% x y z ang -> x' y' z
/rotz {
/theta exch def
/z exch def
/y exch def
/x exch def
x theta cos mul
y theta sin mul sub
x theta sin mul
y theta cos mul add
z
} def

% Eye coords
/ex .2 def %a little x-y skew adds "drama"
/ey .2 def
/ez 5 def

% x y z -> X Y
/project {
ang roty
ang .25 mul rotx
/z exch def
/y exch def
/x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
} def

10 10 360 {
/ang exch def
%matrix currentmatrix
300 400 translate
100 100 scale
60 pointfield
%setmatrix fill
/flushpage where {pop flushpage} if
showpage
} for

[tlvp]

On Thu, 23 Aug 2012 13:00:15 -0700 (PDT), luser- -droog wrote:

> I've made a little progress here. It still generates the points
> with the O(n^3) brute force algorithm, but it saves the points
> to a data file called "stein.pts" so subsequent pages are much
> less agonizing.

Ah, yes, why carry out the same calculations at every re-run-through?

> Rotations and projections are applied to each point, and then
> it draws a little circle to represent the projected point.
>
> %!
> % Point-field sampling
> % with data caching (in a file),
> % point-wise axial rotations,

> % and perspective projection. ...

Details suppressed in this reply, as I don't speak to them directly.

> ... showpage
> } for

That's neat, luser -drug, many thanks; and the rapid succession of images
helps fill in some details the imagination requires.

What are the chances, though, that, in each successive sweeping out of the
next image, the color used at each given moment can be somehow keyed to
that moment? By that I mean (dividing the time-parametrization interval for
each sweeping out into 4 main subintervals [0, 1/4], [1/4, .5], [.5, .3/4],
and [3/4, 1], to use the colors [rgb] = [FF, (4t)xFF, 0] for t in [0, 1/4],
the colors [rgb] = [(2 - 4t)xFF, FF, 0] for t in [1/4, .5],
the colors [rgb] = [0, FF, (4t - 2)xFF] for t in [.5, .3/4], and
thecolors [rgb] = [0, (4 - 4t]xFF, FF] for t in [3/4, 1],
thereby swooping through the visible spectrum from red = [FF 0 0] through
[FF FF 0] through green = [0 FF 0] through [0 FF FF] to blue = [0 0 FF] on
each "sweep".

I can't quite isolate that sweep-time parameter in your depiction of the
routine, or I'd try to see what that coloration-effect achieves, myself.

[tlvp]

On Fri, 24 Aug 2012 01:45:13 -0400, tlvp failed to clear out a stray
decimal point or two from where it didn't belong -- cleared out now:

> On Thu, 23 Aug 2012 13:00:15 -0700 (PDT), luser- -droog wrote:
>
>> I've made a little progress here. It still generates the points
>> with the O(n^3) brute force algorithm, but it saves the points
>> to a data file called "stein.pts" so subsequent pages are much
>> less agonizing.
>
> Ah, yes, why carry out the same calculations at every re-run-through?
>
>> Rotations and projections are applied to each point, and then
>> it draws a little circle to represent the projected point.
>>
>> %!
>> % Point-field sampling
>> % with data caching (in a file),
>> % point-wise axial rotations,
>> % and perspective projection. ...
>
> Details suppressed in this reply, as I don't speak to them directly.
>
>> ... showpage
>> } for
>
> That's neat, luser -drug, many thanks; and the rapid succession of images
> helps fill in some details the imagination requires.
>
> What are the chances, though, that, in each successive sweeping out of the
> next image, the color used at each given moment can be somehow keyed to
> that moment? By that I mean (dividing the time-parametrization interval for

> each sweeping out into 4 main subintervals [0, 1/4], [1/4, .5], [.5, 3/4],

> and [3/4, 1], to use the colors [rgb] = [FF, (4t)xFF, 0] for t in [0, 1/4],
> the colors [rgb] = [(2 - 4t)xFF, FF, 0] for t in [1/4, .5],

> the colors [rgb] = [0, FF, (4t - 2)xFF] for t in [.5, 3/4], and

> thecolors [rgb] = [0, (4 - 4t]xFF, FF] for t in [3/4, 1],
> thereby swooping through the visible spectrum from red = [FF 0 0] through
> [FF FF 0] through green = [0 FF 0] through [0 FF FF] to blue = [0 0 FF] on
> each "sweep".
>
> I can't quite isolate that sweep-time parameter in your depiction of the
> routine, or I'd try to see what that coloration-effect achieves, myself.

Ah, typos -- sorry about that. Cheers, -- tlvp

[luser- -droog]

Yeah, parametrization is the next step.
And real matrix transformations.
Should be able to radically improve the speed
and clarity. I might scale it back to just 2 cylinders
until I get a better handle on it.

[luser- -droog]

Well, here's one way to add color.
The generator is really just a triple loop

for x = -2.2 .. 2.2 step 1/res
for y = -2.2 .. 2.2 step 1/res
for z = -2.2 .. 2.2 step 1/res
if (f(x,y,z)) { plot(x,y,z) }

So the "sweep" is along the x-axis of the object
(ie. pre-rotation). So if we snag the x coord
and shift and scale it to 0..1, we can get a
rainbow with sethsbcolor (my go-to color trick).
The nice benefit is that it turns pale blue
right in the center where the cross section becomes
a circle. Everywhere else it's a truncated circle
or a square.

I've also reduced the number and size of the points
and tweaked the fuzz to get more moire.

49(1)01:43 AM:ps 0> cat 3d2.ps

%!
% Point-field sampling
% with data caching (in a file),
% point-wise axial rotations,
% and perspective projection

% and color.

/fuzz 10000 def %the "grain" of eq
/oldeq /eq load def
/eq {
fuzz mul round exch fuzz mul round oldeq
} def

/max {
2 copy lt { exch } if pop
} def

%x y z f bool
/f {
dup mul 3 1 roll %z^2 x y
dup mul 3 1 roll %y^2 z^2 x
dup mul 3 1 roll %x^2 y^2 z^2
2 index 2 index add %x^2 y^2 z^2 x^2+y^2
3 index 2 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2
3 index 3 index add %x^2 y^2 z^2 x^2+y^2 x^2+z^2 y^2+z^2
max %x^2 y^2 z^2 x^2+y^2 max(x^2+z^2,y^2+z^2)
max %x^2 y^2 z^2 max(x^2+y^2, max(x^2+z^2,y^2+z^2))
4 1 roll pop pop pop %max(...)
4 eq
} def


/filename (stein.pts) def

/low -2.2 def
/hi 2.2 def

/pointsize { dt .3 mul } def

%generate data by brute force, cache in file, plot
/pointfieldtocache {
/res exch def
/dt 1 res div def

/fuzz res .2 mul def

%fuzz affects the "closeness" of
%the equality test. a lower fuzz will allow more
%values to be equal.
% 'res' gives thin lines
% 'res .5 mul' gives wider "ribbons"

/outfile filename (w) file def
/outbuf 128 string def
low dt hi {
low dt hi {
low dt hi { % xW yW zW "world" coords
3 copy f {
3 copy 3 -1 1 { -1 roll
outbuf cvs outfile exch writestring
outfile (\n) writestring
} for %dump points to file

2 index 2.2 add 4.4 div .5 .6 sethsbcolor

3 copy project
%2 copy exch = =
%2 copy transform exch = = ()=

2 copy 2 copy moveto pointsize 0 360 arc moveto

fill
/flushpage where {pop flushpage} if
} if
pop
} for
pop
} for
pop
} for
outfile closefile
} def

%plot cached data from file
/pointfieldfromcache {
/res exch def
/dt 1 res div def
/infile filename (r) file def
/it 1 def
%cvx exec
%count 3 idiv {

{
{
/it it 1 add def
infile token not {stop} if % bail-out
infile token not {stop} if % on any datafile issues
infile token not {stop} if
%3 copy

2 index 2.2 add 4.4 div .5 .6 sethsbcolor
project
2 copy 2 copy moveto pointsize 0 360 arc moveto

30 pointfield

[tlvp]

On Thu, 23 Aug 2012 23:58:25 -0700 (PDT), luser- -droog wrote:

> ...

> Well, here's one way to add color.
> The generator is really just a triple loop
>
> for x = -2.2 .. 2.2 step 1/res
> for y = -2.2 .. 2.2 step 1/res
> for z = -2.2 .. 2.2 step 1/res
> if (f(x,y,z)) { plot(x,y,z) }
>
> So the "sweep" is along the x-axis of the object
> (ie. pre-rotation). So if we snag the x coord
> and shift and scale it to 0..1, we can get a
> rainbow with sethsbcolor (my go-to color trick).
> The nice benefit is that it turns pale blue
> right in the center where the cross section becomes
> a circle. Everywhere else it's a truncated circle
> or a square.
>
> I've also reduced the number and size of the points
> and tweaked the fuzz to get more moire.
>
> 49(1)01:43 AM:ps 0> cat 3d2.ps
> %!
> % Point-field sampling
> % with data caching (in a file),
> % point-wise axial rotations,
> % and perspective projection
> % and color.
>

> ... [details snipped] ...

>
> /flushpage where {pop flushpage} if
> showpage
> } for

The color adds very nicely to the effect, thank you, moy droog :-) .

Two things keep surprising me, though, in the graphics swept out: first,
the large number of times the eye beholds squares unfolding; and, second,
the total absence of any swept out figures whose silhouette is a disc.

After all, peering straight at the figure directly along any of the three
principal axes (x-axis, y-axis, z-axis) should reveal a circular outline.

The squares are the cartographic "contour level-lines" with respect to
heights measured along those axes, but small squares near the top and
bottom should grow as squares only up to a certain point, and then start to
balloon out to become more and more circle-like, finally becoming a true
circle at mid-height exactly. Why doesn't that strike the eye more
blatantly, at least occasionally?

I know, it sounds like "Mrrble, mrrble, grubs again!" but it's not, really,
I'm just curious :-) ... and hoping maybe to catch a glimpse of a tiny
variation in how to compute (or to color) the lines being drawn, so as
better to capture the spirit of that Steinmetz Solid. And, in the meantime,
I think what you have is great! Again, my thanks!

Cheers, -- tlvp

PS: something about Google Groups that keeps interspersing new blank lines
between each pair of posted lines being quoted? Makes everything you quote
get stretched to twice, or twice-twice, or even twice-twice-twice its
original line-count :-) . No chance of convincing you to use a *proper*
newsreader? like Thunderbird? Dialog? Pan? Agent? MesNews? (or other?)

[Luser droog]

Trying KNode with a third free server.


--- Posted via news://freenews.netfront.net/ - Complaints to ne...@netfront.net ---

[Luser droog]

Well, I think to make any real progress, we need to transform it into
an easier-to-calculate form. Casselman frequently uses parametrization,
so I thought that was the thing to do. So,

Max(x^2+y^2, y^2+z^2, x^2+z^2)=4

If we take x=r*cos t, y=r*sin t, then for one "third" of the shape,

r^2 * cos^2 t + r^2 * sin^2 t = 4
r^2 (sin^2 t + cos^2 t) = 4
r = +/- 2

and 0 < t < 2pi

That reduces it to 2 variables t and z, right?
But I'm not sure what to do next. Keep z as it is?
Or convert it to an angle, too?

I think I'll go back to the simple side and try to draw
the cylinders themselves as wireframes (circles one way,
rectangles the other two). Then I should be able to
draw each one constrained by the other two, I think.

[tlvp]

Here's another way to think of the geometric object itself.

Start by picturing a cube -- sides of length 4, say -- with spindle-shafts
drilled perpendicularly through it at each face.

Now chuck it (with the help of a spindle through one of those shafts), and
lathe-turn it down until everything at distance 2 or more from the center
of rotation is cut away. See it? A cylinder of height 4 and diameter 4.

Now pick a new shaft, chuck it with help of a spindle through that shaft,
and turn it down until everything at distance 2 or more from the center of
rotation is cut away. See it? A "roundly beveled former cylinder ... ."

Now, using a spindle through the last shaft, chuck it again, and proceed as
earlier, i.e., turn it down until everything at distance 2 or more from the
center of rotation is cut away. See it? The Steinmetz Solid.

Think there's any hope of showing how it's gradually revealed through that
last lathe-turning procedure?

Me, I'm just not up to keeping tabs on all the recursion/iteration needed.

Cheers, -- tlvp

--
Avant de repondre, jeter la poubelle, SVP.

[PS: how sweet to have no pointless blank lines quoted in follow-ups :-) .]

[Luser droog]

I can definitely imagine it better after reading this.
But drawing it? Not so sure. To simulate a lathe, we need
a material that can be cut and removed. That would lead
me back toward the point-field. Then we could maybe animate
a tool revolving and deactivate the corresponding points in
the field. But I think it would be painfully slow.

This one simplifies (perhaps too much) the shape to projected
polygon rings. Imagine a big hotel atrium with a suspended
sculpture of hula-hoops and you're ascending upon it in the
glass elevator.

This looks worse than the previous, but I still consider it
a step forward because the calculation is more elegant.
Most of the execution time of the other program is calculating
points that fail to show up. Because it has to make sure not
to draw them.

Now with this wireframe, I think I can evolve it to "fill-in"
the gaps and make rectangles out of points in corresponding
rings. Then we can get a normal vector for the "face" and
only draw the "front" of the shape.

The other -lemma is finally learning how those matrix projections
work. And to squash the rotation sequence into a single matrix.
That seems to call for a fully-general matrix multiplication
routine. And a utility to do transposes. Of course I know these
things exist, but to "own" it, I have to write my own. Even if it's
straight from some wikipedia pseudocode. It's like a Hunter Thompson
thing. Where he re-typed the works of Hemmingway to know what
it felt like to write those books. Once I've written my own version
of something, I'm on the inside of it rather than the outside.
Like gematria.

Anyhoo. On to the picture.

554(1)03:41 PM:ps 0> cat 3d3.ps
%!
/.setlanguagelevel where { pop 2 .setlanguagelevel } if

% x y z ang -> x y' z'

/rotx { 3 dict begin

/theta exch def
/z exch def
/y exch def
y theta cos mul
z theta sin mul sub
y theta sin mul
z theta cos mul add

end } def

% x y z ang -> x' y z'

/roty { 4 dict begin

/theta exch def
/z exch def
/y exch def
/x exch def
x theta cos mul
z theta sin mul add
y
x theta sin mul neg
z theta cos mul add

end } def

% x y z ang -> x' y' z

/rotz { 4 dict begin

/theta exch def
/z exch def
/y exch def
/x exch def
x theta cos mul
y theta sin mul sub
x theta sin mul
y theta cos mul add
z

end } def

% Eye coords
/ex .2 def

/ey .2 def
/ez 5 def

% x y z -> X Y
/project {

%ang roty
%ang .25 mul rotx
%alpha rotz
beta roty
gamma rotx
3 dict begin

/z exch def
/y exch def
/x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul

end } def


/moveto { z project //moveto } def
/lineto { z project //lineto } def

%h R N
/cylinder { 2 dict begin
/N exch def
/R exch def
/h exch def
/dz 1 N div def
/dt 360 dz mul def

0 dz 1 {
h mul h 2 div sub /z exch def
R 0 moveto
0 dt 360 {
dup cos R mul
exch sin R mul
lineto
} for
closepath
} for

end } def

%/alpha 0 def
/beta 0 def
/gamma 0 def
0 10 180 { /ang exch def
/ey ang 180 div 5 mul 2.5 sub def

/dep 30 def
300 400 translate
60 dup dup scale 1 exch div setlinewidth

%/alpha ang def
/beta 0 def
/gamma 0 def
4 2 dep cylinder

%/alpha ang def
/beta 90 def
/gamma 0 def
4 2 dep cylinder

%/alpha ang def
/beta 0 def
/gamma 90 def
4 2 dep cylinder

%gsave %false upath
%erasepage
%grestore %uappend
stroke

/flushpage where { pop flushpage } if

%copypage
showpage

} for

[Luser droog]

I've added the parallels by drawing each segment as a little box.
Sadly, this makes it about 4x slower.

I have an intuition that eofill can discern between the the front
and back rectangles since their orientations are quite literally
reversed. But that doesn't help me avoid stroking them.

566(0)12:22 PM:ps 0> cat 3d3.ps

/cylinder { 5 dict begin

/N exch def
/R exch def
/h exch def
/dz 1 N div def
/dt 360 dz mul def

/hdz h dz mul def

0 dz 1 dz sub {

h mul h 2 div sub /z exch def

%R 0 moveto

0 dt 360 {
dup
dup cos R mul
exch sin R mul

2 copy moveto % R t z
/z z hdz add def
lineto % R t z+dz
dt sub

dup cos R mul
exch sin R mul

2 copy lineto % R t-dt z+dz
/z z hdz sub def
lineto % R t-dt z
closepath
} for
%closepath

} for

end } def

%/alpha 0 def
/beta 0 def
/gamma 0 def
0 10 180 { /ang exch def
/ey ang 180 div 5 mul 2.5 sub def

/dep 20 def

[Luser droog]

Well, something's gone very wrong here. But it's the failures
that are most interesting, right?

The last thing I did was add the 'model' calls into the vectors
and out of the 'project' procedure. But that seems right.
I need to apply the model rotation before I take the normal vector
for visibility. Maybe I need to normalize the normal.
or the norml :)

I've inlined the matrix stuff, to save y'all some cut-and-paste
or directory hunting or usenet archiving or giving up before you
can laugh at this poor excuse for a cylinder. Somebody left it
on the railroad tracks.

%!
%The ol' cylinder, take 4

%(mat.ps) run
%!
%mat.ps
%Matrix and Vector math routines

/.error where { pop /signalerror { .error } def } if

/dot { % u v
2 copy length exch length ne {
/dot cvx /undefinedresult signalerror
} if
% u v
0 % u v sum
0 1 3 index length 1 sub { % u v sum i
3 index 1 index get exch % u v sum u_i i
3 index exch get % u v sum u_i v_i
mul add % u v sum
} for % u v sum

3 1 roll pop pop % sum
} bind def

% [ x1 x2 x3 ] [ y1 y2 y3 ] cross [ x2*y3-y2*x3 x3*y1-x1*y3 x1*y2-x2*y1 ]
/cross { % u v
dup length 3 ne
2 index length 3 ne or {
/cross cvx /undefinedresult signalerror
} if
% u v
exch aload pop 4 3 roll aload pop % x1 x2 x3 y1 y2 y3
[
5 index 2 index mul % ... [ x2*y3
3 index 6 index mul sub % ... [ x2*y3-y2*x3
5 index 5 index mul % ... [ x2*y3-y2*x3 x3*y1
8 index 4 index mul sub % ... [ x2*y3-y2*x3 x3*y1-x1*y3
8 index 5 index mul % ... [ x2*y3-y2*x3 x3*y1-x1*y3 x1*y2
8 index 7 index mul sub % ... [ x2*y3-y2*x3 x3*y1-x1*y3 x1*y2-x2*y1
]
7 1 roll 6 { pop } repeat
} bind def

/transpose { STATICDICT begin
/A exch def
/M A length def
/N A 0 get length def
[
0 1 N 1 sub { /n exch def
[
0 1 M 1 sub { /m exch def
A m get n get
} for
]
} for
]
end } dup 0 6 dict put def

/matmul { STATICDICT begin
/B exch def
B 0 get type /arraytype ne { /B [B] def } if
/A exch def
A 0 get type /arraytype ne { /A [A] def } if
/Q B length def
/R B 0 get length def
/P A length def
Q A 0 get length ne {
/A A transpose def
/P A length def
Q A 0 get length ne {
A B end /matmul cvx /undefinedresult signalerror
} if
} if

[
0 1 R 1 sub { /r exch def
[
0 1 P 1 sub { /p exch def
0
0 1 Q 1 sub { /q exch def
A p get q get
B q get r get mul
add
} for
} for
]
} for
]

end } dup 0 10 dict put def

%u v {operator} vop u(op)v
%apply a binary operator to corresponding elements
%in two vectors producing a third vector as result
/vop { 1 dict begin
/op exch def
2 copy length exch length ne {
/vop cvx end /undefinedresult signalerror
} if

[ 3 1 roll % [ u v
0 1 2 index length 1 sub { % [ ... u v i
3 copy exch pop get % u v i u_i
3 copy pop get % u v i u_i v_i
op exch pop % u v u_i(op)v_i
3 1 roll % u_i(op)v_i u v
} for % [ ... u v
pop pop ]

% x y z -> x' y' z'
/model {

%ang roty
%ang .25 mul rotx
%alpha rotz
beta roty

%gamma rotx
} def

% Eye coords
/ex 0 def
/ey 0 def

/ez 5 def

% x y z -> X Y
/project {

3 dict begin
/z exch def
/y exch def
/x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
end } def

%h R N
/cylinder { 20 dict begin

/N exch def
/R exch def
/h exch def
/dz 1 N div def
/dt 360 dz mul def
/hdz h dz mul def

0 dz 1 dz sub {
h mul h 2 div sub /z exch def

0 dt 360 { /t exch def
/v1 [ t cos R mul
t sin R mul
z model ] def
/v4 [ v1 aload pop pop
z hdz add model ] def
/t t dt add def
/v2 [ t cos R mul
t sin R mul
z model ] def
/v3 [ v2 aload pop pop
z hdz add model ] def
/normal v4 v1 {sub} vop
v2 v1 {sub} vop
cross def
[normal aload pop 1] [ex ey ez 1] dot 0 ge {
/action { moveto /action { lineto } def } def
[ v1 v2 v3 v4 ]
{ aload pop project action }
forall
closepath stroke
} if

} for
} for
end } def

300 400 translate
1 50 dup dup scale div setlinewidth
/beta 90 def
4 2 10 cylinder

[M. Joshua Ryan]

Here's the stupid. ^ the x and y of v4 and v3 get modelled twice.

> /t t dt add def
> /v2 [ t cos R mul
> t sin R mul
> z model ] def
> /v3 [ v2 aload pop pop
> z hdz add model ] def

So take the 'model's out of the above and add this here:

[ v1 v2 v3 v4 ] {
aload 4 1 roll model 4 3 roll astore pop
} forall

> /normal v4 v1 {sub} vop
> v2 v1 {sub} vop
> cross def

I also tried normalizing the normal, but, of course it doesn't
change anything. The dot product still has the same sign.

/nlen 0 normal { dup mul add } forall def
/normal normal [nlen nlen nlen] {div} vop def

> [normal aload pop 1] [ex ey ez 1] dot 0 ge {

Changing 'ge' to 'le' shows the front of the shape instead of the back.
Since the eye is at +z, the eye-vector is -z.

[luser- -droog]

I snagged a quick-and-dirty light model from Rogers' Procedural Elements.
Take off the wire-frame training wheels, and burn the background to
bring up the white edge and now it looks like a thing!

Still not the "right" thing yet, but we're getting there.

505(1)08:08 PM:ps 0> cat 3d4.ps
%!
%A shaded cylinder! Woohoo!

%length of a vector
/mag { 0 exch { dup mul add } forall } def

gamma rotx
} def

% Eye coords

/ex .1 def
/ey .1 def
/ez 5 def
/eyedir [ex ey ez]
dup mag [ exch dup dup ]{div} vop

def

% x y z -> X Y
/project {
3 dict begin
/z exch def
/y exch def
/x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
end } def

/light
[ 3 -7 -2 1 ]
dup mag [ exch dup dup dup ]{div} vop
def
/Ia .4 def % Incident Ambient Intensity
/Ka .4 def % Ambient Diffuse reflection constant
/Il .5 def % Incident intensity of Lightsource
/Kd .3 def % Diffuse reflection constant

%h R N
/cylinder { 20 dict begin
/N exch def
/R exch def
/h exch def
/dz 1 N div def
/dt 360 dz mul def
/hdz h dz mul def

0 dz 1 dz sub {
h mul h 2 div sub /z exch def

0 dt 360 { /t exch def
/v1 [ t cos R mul
t sin R mul

z ] def

/v4 [ v1 aload pop pop

z hdz add ] def

/t t dt add def
/v2 [ t cos R mul
t sin R mul

z ] def

/v3 [ v2 aload pop pop

z hdz add ] def
[ v1 v2 v3 v4 ] {

aload 4 1 roll model 4 3 roll astore pop
} forall
/normal v4 v1 {sub} vop
v2 v1 {sub} vop
cross def

/nlen normal mag def

/normal normal [nlen nlen nlen] {div} vop def

[normal aload pop 1] [eyedir aload pop 1] dot 0 lt {

/action { moveto /action { lineto } def } def
[ v1 v2 v3 v4 ]
{ aload pop project action }
forall
closepath

gsave
[normal aload pop 1]
light
%[ex ey ez neg 1] %"radiant"
dot
Il Kd mul mul
Ia Ka mul add
setgray
fill
grestore

stroke
} if

} for
} for
end } def

300 400 translate

280 dup dup moveto
dup neg dup neg lineto
dup neg dup lineto
dup neg lineto closepath .6 setgray fill
1 70 dup dup scale div setlinewidth

/beta 0 def
/gamma 0 def

4 2 50 cylinder

/beta 90 def
/gamma 0 def

4 2 50 cylinder

/beta 0 def
/gamma 90 def

4 2 50 cylinder

[luser- -droog]

Here's a small improvement over the previous.
Three cylinders each chopped by the other two.
Color selected by (hue,sat,bright) =
(model-theta, model-z, that-light-model-I-found)
One unfortunate effect is changing the resolution
really messes up your lighting.

I've got it set-up to grow "out of nowhere"
dramatically by setting the depth of the approximation
by the fibonacci sequence. dep=8 is curiously invisible.
And it starts to take shape around 89.


534(1)02:56 PM:ps 0> cat 3d4.ps
%!
%Steinmertz solid (?)

/prep {
transform
2 {
exch
%floor
round
%ceiling
%2 mul cvi 2 div %round
} repeat
itransform

/ex 3 def
/ey 3 def
/ez 10 def

/eyedir [ex ey ez]
dup mag [ exch dup dup ]{div} vop
def

% x y z -> X Y
/project {
3 dict begin
/z exch def
/y exch def
/x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
end } def

/light

[ 3 -7 -4 1 ]
%dup mag [ exch dup dup dup ]{div} vop

def
/Ia .4 def % Incident Ambient Intensity

/Ka .5 def % Ambient Diffuse reflection constant
/Il .3 def % Incident intensity of Lightsource
/Kd .4 def % Diffuse reflection constant

/front { 0 gt } def
/wirecolor
{ .5 .5 .5 sethsbcolor } def %medium-green wirecolor
%{} def %use fill-color as wirecolor
/wire? true def
/fill? true def
%/fudge 1.0 def
/fudge 1.03 def
%/fudge 1.06 def

%h R N
/cylinder { 20 dict begin

%/facen 0 def

/N exch def
/R exch def
/h exch def
/dz 1 N div def
/dt 360 dz mul def

/fdt dt fudge mul def

/hdz h dz mul def

/fhdz hdz fudge mul def

0 dz 1 dz sub {
h mul h 2 div sub /z exch def

0 dt 360 { /t exch def

% Generate a small face on x^2+y^2=R^2
% by x = cos, y = sin, z = z
% w = dz, h ~= k*dt
%/facen facen 1 add def

/v1 [ t cos R mul
t sin R mul
z ] def
/v4 [ v1 aload pop pop

z fhdz add ] def
/t t fdt add def

/v2 [ t cos R mul
t sin R mul
z ] def
/v3 [ v2 aload pop pop

z fhdz add ] def

% Check if face is inside
% x^2+z^2=R^2 and y^2+z^2=R^2
true [ v1 v2 v3 v4 ] {
aload pop
3 { 3 1 roll dup mul } repeat % x^2 y^2 z^2
dup 4 1 roll % z x y z
R dup mul dup 6 1 roll % r z x y z r
3 1 roll add exch lt % r z x yz<r
4 1 roll add exch lt % yz<r zx<r
and and
} forall % bool %within bounds of the other 2 cylinders

{ %if

% Perform model->world transform

[ v1 v2 v3 v4 ] {
aload 4 1 roll model 4 3 roll astore pop
} forall

% Calculate normal vector in World coords

/normal v4 v1 {sub} vop
v2 v1 {sub} vop
cross def
/nlen normal mag def
/normal normal [nlen nlen nlen] {div} vop def

% Check visibility of a point on the face
%/front { 0 gt } def
normal v1 [ex ey ez] {sub} vop dot
front % bool %facing the eye?

% Draw the face
{ %if
%[ v1 v2 v3 v4 ] ==
% Draw the projected path

/action { moveto /action { lineto } def } def

%facen 2 mod 0 eq {
[ v1 v2 v3 v4 ]
%} {[ v1 v4 v3 v2 ]} ifelse
{ aload pop project
%2 copy 2 array astore ==
prep
action }
forall
closepath
%()= flush


% Calculate illumination
t 360 div % hue
z h 2 div add h div % sat

[normal aload pop 1]
light
%[ex ey ez neg 1] %"radiant"
dot
Il

1 light 0 3 getinterval
v1 { sub } vop mag div mul

Kd mul mul
Ia Ka mul add

%.2 add
dup 1 gt { pop 1 } if
%setgray % bright
sethsbcolor

% Draw according to control bools
wire? { gsave wirecolor stroke grestore } if
fill? { gsave
fill
%/flushpage where{pop flushpage}if
grestore } if
newpath

} if % visible

} if % inside other two cylinders

} for % dt 0..360
} for % dz
end } def % cylinder

/depp 0 def
/dep 1 def
{
%/dep 1 add
/dep dep depp add /depp dep def def %fibonacci
(dep = )print dep =

300 400 translate
280 dup dup moveto
dup neg dup neg lineto
dup neg dup lineto
dup neg lineto closepath
.6 setgray
fill
1 70 dup dup scale div setlinewidth

%/front { 0 lt } def

/beta 0 def
/gamma 0 def

4 2 dep cylinder

%/front { 0 gt } def

/beta 90 def
/gamma 0 def

4 2 dep cylinder

/beta 0 def
/gamma 90 def

4 2 dep cylinder

showpage

} loop

[luser- -droog]

Next step is clipping the faces more closely.

%!
%Spinning Steimertz Solid (approx.)

/omega 0 def
/makespinmat {
/spinmat
[[ omega cos 0 omega sin ]
[ 0 1 0 ]
[ omega sin neg 0 omega cos ]]
def
} def
/spin { spinmat matmul } def

[ v1 v2 v3 v4 ] {
aload pop
3 { 3 1 roll dup mul } repeat % x^2 y^2 z^2
dup 4 1 roll % z x y z
R dup mul dup 6 1 roll % r z x y z r
3 1 roll add exch lt % r z x yz<r
4 1 roll add exch lt % yz<r zx<r
and

} forall
4 copy
or or or %true if any inside
5 1 roll
and and and %false if any outside
exch %use any-outside
pop

% bool %within bounds of the other 2 cylinders

{ %if

% Perform model->object transform

[ v1 v2 v3 v4 ] {
aload 4 1 roll model 4 3 roll astore pop
} forall

% Perform object->world transform
[ v1 v2 v3 v4 ] {
dup spin 0 get exch copy pop

v1 { sub } vop mag div mul %scale light by dist

Kd mul mul
Ia Ka mul add
%.2 add
dup 1 gt { pop 1 } if
%setgray % bright
sethsbcolor

% Draw according to control bools
wire? { gsave wirecolor stroke grestore } if
fill? { gsave
fill
%/flushpage where{pop flushpage}if
grestore } if
newpath

} if % visible

} if % inside other two cylinders

} for % dt 0..360
} for % dz
end } def % cylinder

%/depp 0 def
%/dep 1 def
/dep 77 def
{
%/dep 1 add
%/dep dep depp add /depp dep def def %fibonacci
%(dep = )print dep =
/omega omega 10 add def
makespinmat

300 400 translate
280 dup dup moveto
dup neg dup neg lineto
dup neg dup lineto
dup neg lineto closepath
.6 setgray
fill
1 70 dup dup scale div setlinewidth

%/front { 0 lt } def
/beta 0 def
/gamma 0 def
4 2 dep cylinder

%/front { 0 gt } def
/beta 90 def
/gamma 0 def
4 2 dep cylinder

/beta 0 def
/gamma 90 def
4 2 dep cylinder

showpage

pstack

} loop

[luser- -droog]

Trying to add clipping turned out to be a real mess.
So this is a step backwards, but up as well.
Just the wireframe for now, but a much better
separation of functions, and simpler coding overall, I think.

And it caches the faces in a file (self-executing, so it's a
"sprite" as well) that it executes three times for the three
cylinders.

And since I'm assuming file-access, I'm not inlining mat.ps
anymore. It's available in other messages to the group
(and a buggy one on SO).

[PS. I know I've gone back to google, but as long as don't
quote anybody, it's inoffensive, right? I've got to go
newsserver shopping again...]


593(1)09:02 PM:ps 0> cat 3d5.ps
%!
%Steinmertz, take 5: NEEDFORSPEED
(mat.ps)run

/cylxz <<
/inside { % x y z . bool
exch pop
dup mul exch dup mul add R^2 le } bind
>> def

/cylyz <<
/inside { % x y z . bool
dup mul exch dup mul add R^2 le
exch pop } bind
>> def

/checkface { % [ v1 v2 v3 v4 ] dict
begin
{ aload pop inside } forall
4 copy eq eq eq { % all-eq
pop pop pop
}{ % uneq in>=1 out>=1
pop pop pop %TODO: clipface
} ifelse
end } bind def

/writeface {
out ([) writestring
[ v1 v2 v3 v4 ]{ % forall vertices
out ([) writestring
{ % forall coords of vertex
out exch outbuf cvs writestring
out ( ) writestring
} forall
out ( ) writestring
out (]) writestring
} forall
out (] drawface\n) writestring
out flushfile
} def

/steinmertz-gen { 20 dict begin

/N exch def
/R exch def

/R^2 R dup mul def

/h exch def
/dz 1 N div def
/dt 360 dz mul def

/hdz h dz mul def

0 dz 1 dz sub {
h mul h 2 div sub /z exch def

0 dt 360 dt sub { /t exch def

/v1 [ t cos R mul
t sin R mul
z ] def

/v4 [ t cos R mul
t sin R mul
z hdz add ] def
/t t dt add def

/v2 [ t cos R mul
t sin R mul
z ] def

/v3 [ t cos R mul
t sin R mul
z hdz add ] def

[ v1 v2 v3 v4 ] cylxz checkface {
[ v1 v2 v3 v4 ] cylyz checkface {
writeface
} if
} if
} for
} for
end } bind def


% Generate Data File to cache the faces

/data (stein.fac) def
{
/in data (r) file def
} stopped {
clear
/out data (w) file def
/outbuf 128 string def

% Generate the faces
%h R N
4 2 21 steinmertz-gen

flush
out closefile
/out null def

/in data (r) file def
} if
/reopen { in closefile
/in data (r) file def } def




/drawface { DICT begin % [[ v1 v2 v3 v4 ]]

/action { moveto /action { lineto } def } def

{ % [ x y z ]
MO matmul %[[x' y' z']]
0 get aload pop %x' y' z'
proj %X Y
action
} forall
closepath stroke
end } dup 0 2 dict put bind def

/I3 [[1 0 0]
[0 1 0]
[0 0 1]] def
/MO I3 def %default model-object transform

/E [ 0 0 10 ] def %eye point
/crackE { % set pointers into E
/ex E 0 1 getinterval cvx def
/ey E 1 1 getinterval cvx def
/ez E 2 1 getinterval cvx def
} def crackE
/E^ E dup [ exch mag neg dup dup ] { div } vop def %eye center vec

/proj { DICT begin /z exch def /y exch def /x exch def

1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul

end } dup 0 10 dict put bind def

300 400 translate
1 90 dup dup scale div setlinewidth

in cvx exec reopen
/MO [[1 0 0]
[0 0 1]
[0 1 0]] def % 90 rotx
in cvx exec reopen
/MO [[0 0 1]
[0 1 0]
[1 0 0]] def % 90 roty
in cvx exec

showpage

[tlvp]

On Thu, 13 Sep 2012 19:07:37 -0700 (PDT), luser- -droog wrote:

> ... I've got to go
> newsserver shopping again ...

http://aioe.org/ ... news.eternal-september.org ... freenews.netfront.net

There are others, too. HTH. Cheers, -- tlvp

[luser.droog]

tlvp wrote:

> On Thu, 13 Sep 2012 19:07:37 -0700 (PDT), luser- -droog wrote:
>
>> ... I've got to go
>> newsserver shopping again ...
>
> http://aioe.org/ ... news.eternal-september.org ... freenews.netfront.net
>
> There are others, too. HTH. Cheers, -- tlvp

Awesome! thnx. And eternal-september has its own calendrical system!

So I did some further reorganizing and added normal-visibility and
simplified the model->object transform. Even without lighting and color,
I think it looks pretty cool now.

625(1)04:49 PM:ps 1> cat 3d5.ps

%!
%Steinmertz, take 5: NEEDFORSPEED
(mat.ps)run

%/forall { pstack()= forall } bind def

%The horizontal chopping cylinder

/cylxz <<
/inside { % x y z . bool
exch pop
dup mul exch dup mul add R^2 le } bind
>> def

%The vertical chopping cylinder

/cylyz <<
/inside { % x y z . bool
dup mul exch dup mul add R^2 le
exch pop } bind
>> def

%Check that the vertices of a face
%are inside a specified chopping cylinder

/checkface { % [ v1 v2 v3 v4 ] dict
begin
{ aload pop inside } forall

4 copy eq 3 1 roll eq eq { % all-eq

pop pop pop
}{ % uneq in>=1 out>=1
pop pop pop %TODO: clipface

%pop false
} ifelse
end } bind def

%Write the vertex array [v1 v2 v3 v4]
%to outfile with embedded drawing command

/writeface {
out ([) writestring
[ v1 v2 v3 v4 ]{ % forall vertices
out ([) writestring
{ % forall coords of vertex
out exch outbuf cvs writestring
out ( ) writestring
} forall
out ( ) writestring
out (]) writestring
} forall
out (] drawface\n) writestring
out flushfile
} def

%Generate the faces of the Forward cylinder
%and eliminate faces outside the chopping cylinders

/steinmertz-gen { 20 dict begin
/N exch def
/R exch def
/R^2 R dup mul def
/h exch def
/dz 1 N div def
/dt 360 dz mul def
/hdz h dz mul def

0 dz 1 dz sub {
h mul h 2 div sub /z exch def

0 dt dup N 3 add mul { /t exch def

/v1 [ t cos R mul
t sin R mul
z ] def
/v4 [ t cos R mul
t sin R mul
z hdz add ] def
/t t dt add def
/v2 [ t cos R mul
t sin R mul
z ] def
/v3 [ t cos R mul
t sin R mul
z hdz add ] def

/face [ v1 v2 v3 v4 ] def

face cylxz checkface {
face cylyz checkface {
doface

} if
} if
} for
} for
end } bind def

%Action performed by steinmertz-gen
%on each face that survives the chopping
/doface {
usecache? {
writeface
}{
face drawface
} ifelse
} def


%This controls the parameters of the
%wireframe approximation of the cylinder
/genfaces {

% Generate the faces
%h R N

4 2 89 steinmertz-gen
} def



%Default modeling transform
/I3 3 ident def
/MO I3 def % model->object
/model {} def
/OW I3 def % object->world

/fill? true def
/wire? true def
/wirecolor
%{} def
{1 setgray} def
%Perform modeling transform
%Perform object->world transform
%Perform perspective projection
%Check visibility
%Fill if fill?
%Draw Outline if wire?

/drawface { DICT begin % [ v1 v2 v3 v4 ]

/face exch def

/action { moveto /action { lineto } def } def

face { % [ x y z ]

aload 4 1 roll model 4 3 roll astore

%MO matmul
dup OW matmul 0 get exch copy pop
} forall

{ %exitloop
face visible not { exit } if
face {
aload pop

proj %X Y
action
} forall
closepath

%set color
fill? { gsave fill grestore } if

wire? { gsave wirecolor stroke grestore } if

newpath
exit } loop

end } dup 0 2 dict put bind def

/E [ 0 0 10 ] def %eye point
/crackE { % set pointers into E
/ex E 0 1 getinterval cvx def
/ey E 1 1 getinterval cvx def
/ez E 2 1 getinterval cvx def
} def crackE
/E^ E dup [ exch mag neg dup dup ] { div } vop def %eye center vec

%Check visibility from the eye
/visible { % [ v1 v2 v3 v4 ] . bool
dup 0 get 1 index 1 get 2 index 3 get % [] v1 v2 v4
2 index { sub } vop % [] v1 v2 v14
3 1 roll exch { sub } vop % [] v14 v12
cross /normal exch def %normal ==
%/normal [ normal mag dup dup ] { div } vop def
dup 0 get E { sub } vop % [] ve1
/ev exch def %ev ==
ev normal dot %dup =
0 gt % [] bool
exch pop
} def

/proj { DICT begin /z exch def /y exch def /x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
end } dup 0 10 dict put bind def

/setuppage {

300 400 translate
1 90 dup dup scale div setlinewidth

} def

/usecache? true def
usecache? { % Generate Data File to cache the faces

/data (stein.fac) def
{
/in data (r) file def
} stopped {
clear
/out data (w) file def
/outbuf 128 string def

genfaces

flush out closefile
/out null def
/in data (r) file def
} if
/reopen { in closefile /in data (r) file def } def

/drawshape { in cvx exec reopen } def
}{ % Generate fresh data for each shape
/drawshape { genfaces } def
} ifelse

/up 90 rotx def
/left 90 roty def

/ang 0 def
{
/OW ang roty def
/ang ang 10 add def

setuppage

/model {} def % x y z -> x y z
%/MO I3 def
drawshape
/model { exch neg } def % x y z -> x z -y
%/MO up def
drawshape
/model { 3 1 roll exch neg } def % x y z -> z y -x
%/MO left def
drawshape

showpage
} loop

[luser.droog]

luser.droog wrote:

> tlvp wrote:
>
>> On Thu, 13 Sep 2012 19:07:37 -0700 (PDT), luser- -droog wrote:
>>
>>> ... I've got to go
>>> newsserver shopping again ...
>>
>> http://aioe.org/ ... news.eternal-september.org ... freenews.netfront.net
>>
>> There are others, too. HTH. Cheers, -- tlvp
>
> Awesome! thnx. And eternal-september has its own calendrical system!
>
> So I did some further reorganizing and added normal-visibility and
> simplified the model->object transform. Even without lighting and color,
> I think it looks pretty cool now.
>
> 625(1)04:49 PM:ps 1> cat 3d5.ps

Back to the drawing board, I realized that nothing really changes along
the "z" axis of a cylinder. So there's really no need to iterate over it.
So I tried one big face per theta but then there's no known "inside"
point to seed the clip. So I tried two faces per theta with vertices like
6 on a domino. Then I realized that I could avoid clipping against one
of the cylinders by simply setting the outer z coords on the face by
solving for z using x or y of one of the cylinders, the x or y already
calculated for the face. After further playing, I realized I could take
both and use the max for the upper z and min for the lower z. Thus the
faces of the cylinder being drawn are bounded by the other two.
No need to clip. But I leave it here anyway so you all can see how
heinous doing vector clipping is (assuming a demonstration of this fact
is at all necessary to anyone).

But while most of the shape is simple and beautiful, there are distracting
artifacts at the 4-face joins that look like little swastikas a higher
resolutions. And that just sucks. Maybe I should ask for help with that
on SO... (???)

Anyway, I snagged a simpler light calculation from Casselman than the
one from Rogers, and set the hue to the x/y angle of v1 of the face,
const .5 saturation. And wow! If you can overlook (or remove! ?)
the swastikas, this is quite pretty to look at. Maybe a dark background
would help...

522(1)11:30 PM:ps 0> cat 3d5.ps

%!
%Steinmertz, take 5: NEEDFORSPEED

% O(N) generator makes caching largely irrelevant!
% Bounded faces makes clipping unnecessary
% (which wasn't quite working somehow)
% Improved visibility now shows the Front of the object!
% TODO: eliminate those ugly swastika artifacts

(mat.ps)run

%/forall { pstack()= forall } bind def

%/if { pstack()= if } bind def

/min { 2 copy gt { exch } if pop } def
/max { 2 copy lt { exch } if pop } def

%The vertical chopping cylinder
/cylyz <<
/inside { % x y z . bool
dup mul exch dup mul add R^2 le
exch pop } bind

/doparams {
/D edge 1 get edge 0 get { sub } vop
[ 1 index mag dup dup ] { div } vop def
/a D 1 get dup mul D 2 get dup mul add def
/b 2 edge 0 get 1 get D 1 get mul mul
2 edge 0 get 2 get D 2 get mul mul add def
/c edge 0 get 1 get dup mul
edge 0 get 2 get dup mul add 1 sub def
/t [ b dup mul 4 a c mul mul sub dup 0 lt { neg } if
sqrt 2 a mul div dup neg
b neg 2 a mul div add exch
b neg 2 a mul div add ] def }

>> def

%The horizontal chopping cylinder
/cylxz <<
/inside { % x y z . bool
exch pop
dup mul exch dup mul add R^2 le } bind

/doparams {
/D edge 1 get edge 0 get { sub } vop
%[ 1 index mag dup dup ] { div } vop
def
/a D 0 get dup mul D 2 get dup mul add def
/b 2 edge 0 get 0 get D 0 get mul mul
2 edge 0 get 2 get D 2 get mul mul add def
/c edge 0 get 0 get dup mul
edge 0 get 2 get dup mul add 1 sub def
/t [ b dup mul 4 a c mul mul sub dup 0 lt { neg } if
sqrt 2 a mul div dup neg
b neg 2 a mul div add exch
b neg 2 a mul div add ] def }
>> def

%cd:chopping-cylinder
/clipface { /i exch def
/v1 face 0 get def
/v2 face 1 get def
/v3 face 2 get def
/v4 face 3 get def
[ [v1 v2] [v4 v3] [v1 v2] [v4 v3] [v1 v2] [v2 v1] [v3 v4] ]
i 2 getinterval { /edge exch def
edge 1 get aload pop inside not {
doparams
t aload pop
2 copy 0 lt exch 0 lt xor { dup 0 ge { exch } if pop }{

2 copy lt { exch } if pop

} ifelse
dup 1 gt { pop 1.01 } if
[ exch dup dup ] D { mul } vop
edge 0 get { add } vop
edge 1 get copy pop
} if
} forall
true

} def

%Check that the vertices of a face
%are inside a specified chopping cylinder
/checkface { % [ v1 v2 v3 v4 ] dict
begin
{ aload pop inside } forall

%pop true false false true
%4 copy eq 3 1 roll eq eq { % all-eq
% pop pop pop
%}{ % uneq in>=1 out>=1
true { %don't clip
pop pop pop pop true
}{
{ %if v4true
pop pop pop
3 clipface
}{ %else
{ %if v3true
pop pop
2 clipface
}{ %else
{ %if v2true
pop
1 clipface
}{ %else
{ %if v1true
0 clipface
}{ %else none true
false
} ifelse
} ifelse
} ifelse
} ifelse % i(v_i inside)
} ifelse
%} ifelse

end } bind def

%Write the vertex array [v1 v2 v3 v4]
%to outfile with embedded drawing command
/writeface {
out ([) writestring
[ v1 v2 v3 v4 ]{ % forall vertices
out ([) writestring
{ % forall coords of vertex
out exch outbuf cvs writestring
out ( ) writestring
} forall
out ( ) writestring
out (]) writestring
} forall
out (] drawface\n) writestring
out flushfile
} def

/fudge 1.03 def

%Generate the faces of the Forward cylinder
%and eliminate faces outside the chopping cylinders
/steinmertz-gen {

/N exch def
/R exch def
/R^2 R dup mul def
/h exch def
/dz 1 N div def
/dt 360 dz mul def
/hdz h dz mul def

0 dt 360 { /t exch def

%0 dt 180 { /u exch def

/v1 [
R t cos mul %u cos mul
R t sin mul %u cos mul
0
%R^2 2 index dup mul sub sqrt neg
%R^2 2 index dup mul sub sqrt neg max
] def
/v2 [ %v1 aload pop neg
R t cos mul
R t sin mul
R^2 2 index dup mul sub sqrt
R^2 2 index dup mul sub sqrt min
] def
/t t dt fudge mul add def
/v3 [
R t cos mul
R t sin mul
R^2 2 index dup mul sub sqrt
R^2 2 index dup mul sub sqrt min
] def
/v4 [ %v3 aload pop neg
R t cos mul
R t sin mul
0

] def
/face [ v1 v2 v3 v4 ] def

face cylyz checkface { doface } if

/v1 [
R t cos mul %u cos mul
R t sin mul %u cos mul
0
%R^2 2 index dup mul sub sqrt neg
%R^2 2 index dup mul sub sqrt neg max
] def
/v2 [ %v1 aload pop neg
R t cos mul
R t sin mul
R^2 2 index dup mul sub sqrt neg
R^2 2 index dup mul sub sqrt neg max
] def
/t t dt fudge mul sub def
/v3 [
R t cos mul
R t sin mul
R^2 2 index dup mul sub sqrt neg
R^2 2 index dup mul sub sqrt neg max
%R^2 2 index dup mul sub sqrt neg max
] def
/v4 [ %v3 aload pop neg
R t cos mul
R t sin mul
0

] def
/face [ v1 v2 v3 v4 ] def

face cylyz checkface { doface } if

%} for
} for

} def

%Action performed by steinmertz-gen
%on each face that survives the chopping
/doface {
usecache? {
writeface
}{
face drawface
} ifelse
} def


%This controls the parameters of the
%wireframe approximation of the cylinder
/genfaces {
% Generate the faces
%h R N

4 2 300 steinmertz-gen

} def



%Default modeling transform
/I3 3 ident def
/MO I3 def % model->object
/model {} def
/OW I3 def % object->world

/fill? true def

/wire? false def

colorface

fill? { gsave fill grestore } if
wire? { gsave wirecolor stroke grestore } if

/flushpage where { pop flushpage } if

newpath
exit } loop
end } dup 0 2 dict put bind def

/E [ 3 3 10 ] def %eye point

/crackE { % set pointers into E
/ex E 0 1 getinterval cvx def
/ey E 1 1 getinterval cvx def
/ez E 2 1 getinterval cvx def
} def crackE
/E^ E dup [ exch mag neg dup dup ] { div } vop def %eye center vec

%Check visibility from the eye
/visible { % [ v1 v2 v3 v4 ] . bool
dup 0 get 1 index 1 get 2 index 3 get % [] v1 v2 v4
2 index { sub } vop % [] v1 v2 v14
3 1 roll exch { sub } vop % [] v14 v12
cross /normal exch def %normal ==
%/normal [ normal mag dup dup ] { div } vop def
dup 0 get E { sub } vop % [] ve1
/ev exch def %ev ==
ev normal dot %dup =

0 lt % [] bool
exch pop
} def

/light [ -5 10 20 ] def

/colorface {
%normal light dot 1 add 3 div setgray
face 0 get aload pop pop exch atan 360 div
.5
normal light dot 1 add 3 div sethsbcolor

} def

/proj { DICT begin /z exch def /y exch def /x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
end } dup 0 10 dict put bind def


/setuppage {
300 400 translate

1 90 dup dup scale div 2 div setlinewidth
} def

/usecache? false def

usecache? { % Generate Data File to cache the faces
/data (stein.fac) def
{
/in data (r) file def
} stopped {

pop pop pop%clear

/out data (w) file def
/outbuf 128 string def
genfaces
flush out closefile
/out null def
/in data (r) file def
} if
/reopen { in closefile /in data (r) file def } def
/drawshape { in cvx exec reopen } def
}{ % Generate fresh data for each shape
/drawshape { genfaces } def
} ifelse

/up 90 rotx def
/left 90 roty def

/ang 0 def
{
/OW ang roty def
/ang ang 10 add def

setuppage

/model {} def % x y z -> x y z
%/MO I3 def
drawshape

/model { neg exch } def % x y z -> x -z y
%/MO up def
drawshape
/model { neg 3 1 roll exch } def % x y z -> -z y x
%/MO left def
drawshape

pstack flush
showpage
} loop

[luser.droog]

luser.droog wrote:

> 522(1)11:30 PM:ps 0> cat 3d5.ps
> %!
> %Steinmertz, take 5: NEEDFORSPEED
> % O(N) generator makes caching largely irrelevant!
> % Bounded faces makes clipping unnecessary
> % (which wasn't quite working somehow)
> % Improved visibility now shows the Front of the object!
> % TODO: eliminate those ugly swastika artifacts

TODO: Start spelling Steinmetz correctly.

[luser.droog]

Done and Done.

The little swastikas were happening at the corners
where the quadrilateral faces degenerated to triangles.
That meant that one or the other of the edges being
used to calculate the normal vector was of length 0.
Thus the normal vector was 0, which produces a zero
dot product against any vector (like the vector from
the eye).

I fixed it by checking for degenerate edges and
selecting a different vertex if it is. There's still
a missing triangle at each of the 3-point joins.
But that doesn't bother me so much.


533(1)11:30 AM:ps 0> cat 3d5.ps
%!
%Steinmetz, take 5: NEEDFORSPEED

% O(N) generator makes caching largely irrelevant!
% Bounded faces makes clipping unnecessary
% (which wasn't quite working somehow)
% Improved visibility now shows the Front of the object!

% Eliminated those ugly swastika artifacts
(mat.ps)run

/min { 2 copy gt { exch } if pop } def
/max { 2 copy lt { exch } if pop } def

%Write the vertex array [v1 v2 v3 v4]
%to outfile with embedded drawing command
/writeface {
out ([) writestring
[ v1 v2 v3 v4 ]{ % forall vertices
out ([) writestring
{ % forall coords of vertex
out exch outbuf cvs writestring
out ( ) writestring
} forall
out ( ) writestring
out (]) writestring
} forall
out (] drawface\n) writestring
out flushfile
} def

/fudge 1.03 def
%Generate the faces of the Forward cylinder
%and eliminate faces outside the chopping cylinders

/steinmetz-gen {

%face cylyz checkface {
doface
%} if
%face ==

%face cylyz checkface {
doface
%} if

%} for
} for
} def

%Action performed by steinmetz-gen

%on each face that survives the chopping
/doface {
usecache? {
writeface
}{
face drawface
} ifelse
} def


%This controls the parameters of the
%wireframe approximation of the cylinder
/genfaces {
% Generate the faces
%h R N

4 2 500 steinmetz-gen

} def



%Default modeling transform
/I3 3 ident def
/MO I3 def % model->object
/model {} def
/OW I3 def % object->world

/fill? true def

/wire? true def
/wirecolor
%{} def

%{1 setgray} def
{normal light dot 1 add 4 div setgray} def

aload pop /v4 exch def /v3 exch def /v2 exch def /v1 exch def
%dup 0 get 1 index 2 get 2 index 3 get % [] v1 v2 v4
%2 index { sub } vop % [] v1 v2 v14
%3 1 roll exch { sub } vop % [] v14 v12
%cross
v1 v4{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if
v1 v2{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if

cross
/normal exch def %normal ==
%/normal [ normal mag dup dup ] { div } vop def
dup 0 get E { sub } vop % [] ve1
/ev exch def %ev ==
ev normal dot %dup =
0 lt % [] bool
exch pop
} def

/light [ -5 10 20 ] def

/colorface {
%normal light dot 1 add 3 div setgray

face 0 get 45 rotz matmul 0 get
aload pop pop %exch
dup 0 eq { pop .001 } if
atan 360 div 3 div fn .33 mul add
%.5
face 0 get 60 rotz matmul 0 get
aload pop pop exch
dup 0 eq { pop .001 } if
atan 360 div 2 div .25 add

normal light dot 1 add 3 div

3 { 3 1 roll dup 1 gt { pop 1 } if dup 0 lt { pop 0 } if } repeat

/fn 0 def

/model {} def % x y z -> x y z
%/MO I3 def
drawshape

/fn 1 def

/model { neg exch } def % x y z -> x -z y
%/MO up def
drawshape

/fn 2 def

[luser.droog]

Added identity and rotation(3D) matrices.

%!
%mat.ps
%Matrix and Vector math routines

/ident { 1 dict begin /n exch def
[
1 1 n { % [ i
[ exch % [ [ i
1 1 n { % [ [ i j
1 index eq { 1 }{ 0 } ifelse % [ [ i b
exch % [ [ b i
} for % [ [ b+ i
pop ] % [ [ b+ ]
} for % [ [b+]+ ]
]
end } def

/rotx { 1 dict begin /t exch def
[ [ 1 0 0 ]
[ 0 t cos t sin neg ]
[ 0 t sin t cos ] ]
end } def

/roty { 1 dict begin /t exch def
[ [ t cos 0 t sin ]
[ 0 1 0 ]
[ t sin neg 0 t cos ] ]
end } def

/rotz { 1 dict begin /t exch def
[ [ t cos t sin neg 0 ]
[ t sin t cos 0 ]
[ 0 0 1 ] ]
end } def

0 1 P 1 sub { /p exch def % rows of A
[
0 1 R 1 sub { /r exch def % cols of B
0
0 1 Q 1 sub { /q exch def % terms of sum

[luser.droog]

Pretty much done with this, I think. Still, I kept playing with it.
Since the colors were set by the position of the face and the thing has
4x radial symmetry, it didn't take a full rotation to "see the whole thing".

So I jazzed it up!

I've got it set up now to change the settings to a random configuration
for each cylinder. And it draws the back first and then the front.

So run it with 'gs -dNOGC -dNOPAUSE 3d5a.ps' and you'll visualize the
solid and the circles and the squares in no time flat.

So, hold onto yer socks.

%!
%Steinmetz, take 5: NEEDFORSPEED
% O(N) generator makes caching largely irrelevant!
% Bounded faces makes clipping unnecessary
% (which wasn't quite working somehow)
% Improved visibility now shows the Front of the object!

% Eliminated those ugly artifacts
% Random Configurations (with optional "less-random" mode)

(mat.ps)run

/min { 2 copy gt { exch } if pop } def
/max { 2 copy lt { exch } if pop } def

%Write the vertex array [v1 v2 v3 v4]
%to outfile with embedded drawing command
/writeface {
out ([) writestring
[ v1 v2 v3 v4 ]{ % forall vertices
out ([) writestring
{ % forall coords of vertex
out exch outbuf cvs writestring
out ( ) writestring
} forall
out ( ) writestring
out (]) writestring
} forall
out (] drawface\n) writestring
out flushfile
} def

/fudge

1 def
%1.03 def

4 2 300 steinmetz-gen

} def



%Default modeling transform
/I3 3 ident def
/MO I3 def % model->object
/model {} def
/OW I3 def % object->world

/fill? true def

/fillcolor {} def
/wire? false def
/wirecolor
%{} def
%{0 setgray} def

%{1 setgray} def
{normal light dot 1 add 4 div

%1 exch sub

setgray} def
%Perform modeling transform
%Perform object->world transform
%Perform perspective projection
%Check visibility
%Fill if fill?
%Draw Outline if wire?
/drawface { DICT begin % [ v1 v2 v3 v4 ]
/face exch def
/action { moveto /action { lineto } def } def

face { % [ x y z ]
aload 4 1 roll model 4 3 roll astore
%MO matmul
dup OW matmul 0 get exch copy pop
} forall

{ %exitloop
face visible not { exit } if
face {
aload pop
proj %X Y
action
} forall
closepath

colorface fillcolor

fill? { gsave fill grestore } if
wire? { gsave wirecolor stroke grestore } if
/flushpage where { pop flushpage } if
newpath
exit } loop
end } dup 0 2 dict put bind def

/E [ 3 3 10 ] def %eye point
/crackE { % set pointers into E
/ex E 0 1 getinterval cvx def
/ey E 1 1 getinterval cvx def
/ez E 2 1 getinterval cvx def
} def crackE
/E^ E dup [ exch mag neg dup dup ] { div } vop def %eye center vec

/vistest { lt } def

%Check visibility from the eye
/visible { % [ v1 v2 v3 v4 ] . bool
dup
aload pop /v4 exch def /v3 exch def /v2 exch def /v1 exch def
%dup 0 get 1 index 2 get 2 index 3 get % [] v1 v2 v4
%2 index { sub } vop % [] v1 v2 v14
%3 1 roll exch { sub } vop % [] v14 v12
%cross
v1 v4{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if
v1 v2{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if
cross
/normal exch def %normal ==
%/normal [ normal mag dup dup ] { div } vop def
dup 0 get E { sub } vop % [] ve1
/ev exch def %ev ==
ev normal dot %dup =

0 vistest % [] bool
exch pop
} def

/light [ -7 12 30 ] def

/colorface {
%normal light dot 1 add 3 div setgray

%face 0 get 90 rotz matmul 0 get
%aload pop pop %exch
%dup 0 eq { pop .001 } if
%atan 360 div 3 div
fn .33 mul %add
%.5
face 0 get -90 rotz matmul 0 get

aload pop pop exch
dup 0 eq { pop .001 } if
atan 360 div 2 div .25 add
normal light dot 1 add 3 div
3 { 3 1 roll dup 1 gt { pop 1 } if dup 0 lt { pop 0 } if } repeat
sethsbcolor
} def

/proj { DICT begin /z exch def /y exch def /x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul

end } dup 0 10 dict put bind def


/setuppage {
310 440 translate
1 120 dup dup scale div 2 div setlinewidth

} def

/usecache? false def
usecache? { % Generate Data File to cache the faces
/data (stein.fac) def
{
/in data (r) file def
} stopped {
pop pop pop%clear
/out data (w) file def
/outbuf 128 string def
genfaces
flush out closefile
/out null def
/in data (r) file def
} if
/reopen { in closefile /in data (r) file def } def
/drawshape { in cvx exec reopen } def
}{ % Generate fresh data for each shape
/drawshape { genfaces } def
} ifelse

/up 90 rotx def
/left 90 roty def

/front { /vistest { lt } def } def
/back { /vistest { gt } def } def

/drawsolid {

/fn 0 def
/model {} def % x y z -> x y z

drawshape

/fn 1 def
/model { neg exch } def % x y z -> x -z y

morerandom? { randconfig } if

drawshape

/fn 2 def
/model { neg 3 1 roll exch } def % x y z -> -z y x

morerandom? { randconfig } if
drawshape
} def

/rbool { rand 2 mod 0 eq } def

/randconfig {
true rbool
rbool { exch } if /fill? exch def /wire? exch def
[{} {} {} {} {} {0 setgray} {0 setgray}
{1 setgray} {1 setgray} {1 setgray}

{normal light dot 1 add 4 div setgray}

{normal light dot 1 add 4 div 1 exch sub setgray}
{ 1 currentgray sub setgray }
{ currenthsbcolor exch pop .5 exch sethsbcolor }
{ currenthsbcolor 3 2 roll .33 add dup 1 gt {1 sub} if
3 1 roll exch pop .5 exch sethsbcolor }
{ currenthsbcolor 1 exch sub sethsbcolor }
{ currentrgbcolor exch 1 exch sub exch setrgbcolor }
{ currenthsbcolor setrgbcolor }
{ currentrgbcolor sethsbcolor }
{ currentrgbcolor 3 1 roll exch setrgbcolor }
//colorface //colorface //colorface //colorface //colorface
] dup
dup length rand exch mod get /wirecolor exch def
dup length rand exch mod get /fillcolor exch def
currentlinewidth 1 rand 2 mod 2 div add mul setlinewidth
/fudge 1 rand 10 mod 100 div add def
} def

/morerandom? true def

/ang 0 def
{
/OW ang roty def

/ang ang 3 add def

setuppage

%/fill? true def
%/wire? true def
%/wirecolor { 1 currentgray sub setgray } def
back
randconfig
drawsolid

%/fill? false def
%/wire? true def
%/wirecolor {} def
%2 setlinewidth
%currentlinewidth 2 mul setlinewidth
front
randconfig
drawsolid

pstack flush %cya

%pause for applause
500000 { 100 dup cos exch sin exch atan pop } repeat
showpage
} loop

[tlvp]

(LOL)! . Suppress that "r", though, and the Smithsonian's Uta C. Merzbach
(she without the "t" before the "z") might take offense at its loss :-) .

Cheers, -- tlvp

PS: I've let the takes through and beyond Take 5 go without comment so far
because you seem to be guiding yourself pretty well on your own, w/o
benefit of my misguided and digressionary kibitzing :-) . -- tlvp

[luser.droog]

No problem there. But if there is anything else you want to see it do,
let me know. If I can figure out more about the projection math,
enough to rotate the projection plane, then I could move that little
tear in the coloring right at the front where the z-axis pokes through.
It doesn't look terrible to me, but it is an infelicity not being able
to control it. As it is, no matter how I rotate the points to select the
hue, there's always that tear at the front.

[tlvp]

That tear at the front is probably another indication of the validity of
what's called the "index theorem (aka Hairy Ball Theorem) for non-vanishing
vector fields on a 2-sphere", sometimes paraphrased as: "Ya can't comb the
hair on a coconut without leaving a swirl" :-) . Your "tear" = the index
theorem's "swirl". (Yeah -- go ahead, Google it, you'll find it :-) !)

Best you can try to do is hide it in back somewhere, where it can't show.

Cheers, -- tlvp

[luser.droog]

I decided to go the other way and make it appear bold and "justified".
This one has a new "cage" mode where the back always has a fill and the
front is only wires. And it draws coordinate axes in the middle.
The linewidth varies so sometimes the interior is more revealed, sometimes
more obscured.

557(1)12:41 AM:ps 0> cat 3d5b.ps
%!
%Steinmetz, take 5: NEEDFORSPEED

% O(N) generator makes caching largely irrelevant!
% Bounded faces makes clipping unnecessary
% (which wasn't quite working somehow)
% Improved visibility now shows the Front of the object!

% Eliminated those ugly artifacts
% Random Configurations (with optional "less-random" mode)
(mat.ps)run

/min { 2 copy gt { exch } if pop } def
/max { 2 copy lt { exch } if pop } def

%Write the vertex array [v1 v2 v3 v4]
%to outfile with embedded drawing command
/writeface {
out ([) writestring
[ v1 v2 v3 v4 ]{ % forall vertices
out ([) writestring
{ % forall coords of vertex
out exch outbuf cvs writestring
out ( ) writestring
} forall
out ( ) writestring
out (]) writestring
} forall
out (] drawface\n) writestring
out flushfile
} def

/fudge

1 def
%1.03 def

%Generate the faces of the Forward cylinder
%and eliminate faces outside the chopping cylinders

/steinmetz-gen {

%face cylyz checkface {
doface
%} if
%face ==

%face cylyz checkface {
doface
%} if

%} for
} for
} def

%Action performed by steinmetz-gen

%on each face that survives the chopping
/doface {
usecache? {
writeface
}{
face drawface
} ifelse
} def


%This controls the parameters of the
%wireframe approximation of the cylinder
/genfaces {
% Generate the faces
%h R N

4 2 300 steinmetz-gen

} def



%Default modeling transform
/I3 3 ident def
/MO I3 def % model->object
/model {} def
/OW I3 def % object->world

/fill? true def

/fillcolor {} def

/wire? false def
/wirecolor
%{} def

%{0 setgray} def
%{1 setgray} def
{normal light dot 1 add 4 div
%1 exch sub

setgray} def
%Perform modeling transform
%Perform object->world transform
%Perform perspective projection
%Check visibility
%Fill if fill?
%Draw Outline if wire?
/drawface { DICT begin % [ v1 v2 v3 v4 ]
/face exch def
/action { moveto /action { lineto } def } def

face { % [ x y z ]
aload 4 1 roll model 4 3 roll astore
%MO matmul
dup OW matmul 0 get exch copy pop
} forall

{ %exitloop
face visible not { exit } if
face {
aload pop
proj %X Y
action
} forall
closepath

colorface fillcolor

fill? { gsave fill grestore } if
wire? { gsave wirecolor stroke grestore } if
/flushpage where { pop flushpage } if
newpath
exit } loop
end } dup 0 2 dict put bind def

/E [ 3 3 10 ] def %eye point
/crackE { % set pointers into E
/ex E 0 1 getinterval cvx def
/ey E 1 1 getinterval cvx def
/ez E 2 1 getinterval cvx def
} def crackE
/E^ E dup [ exch mag neg dup dup ] { div } vop def %eye center vec

/vistest { lt } def

%Check visibility from the eye
/visible { % [ v1 v2 v3 v4 ] . bool
dup

aload pop /v4 exch def /v3 exch def /v2 exch def /v1 exch def
%dup 0 get 1 index 2 get 2 index 3 get % [] v1 v2 v4
%2 index { sub } vop % [] v1 v2 v14
%3 1 roll exch { sub } vop % [] v14 v12
%cross
v1 v4{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if
v1 v2{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if

cross
/normal exch def %normal ==
%/normal [ normal mag dup dup ] { div } vop def
dup 0 get E { sub } vop % [] ve1
/ev exch def %ev ==
ev normal dot %dup =

0 vistest % [] bool
exch pop
} def

/light [ -7 12 30 ] def

/colorface {
%normal light dot 1 add 3 div setgray

%face 0 get 90 rotz matmul 0 get
%aload pop pop %exch
%dup 0 eq { pop .001 } if
%atan 360 div 3 div
fn .33 mul %add
%.5

face 0 get -90 rotz matmul 0 get
aload pop pop exch

dup 0 eq { pop .001 } if
atan 360 div 2 div .25 add

normal light dot 1 add 3 div

3 { 3 1 roll dup 1 gt { pop 1 } if dup 0 lt { pop 0 } if } repeat

sethsbcolor
} def

/proj { DICT begin /z exch def /y exch def /x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
end } dup 0 10 dict put bind def


/setuppage {

310 440 translate
1 120 dup dup scale div 2 div setlinewidth
2 setlinejoin

} def

/usecache? false def
usecache? { % Generate Data File to cache the faces
/data (stein.fac) def
{
/in data (r) file def
} stopped {
pop pop pop%clear
/out data (w) file def
/outbuf 128 string def
genfaces
flush out closefile
/out null def
/in data (r) file def
} if
/reopen { in closefile /in data (r) file def } def
/drawshape { in cvx exec reopen } def
}{ % Generate fresh data for each shape
/drawshape { genfaces } def
} ifelse

/up 90 rotx def
/left 90 roty def

/front { /vistest { lt } def } def
/back { /vistest { gt } def } def

/drawsolid {
/fn 0 def

/model {} def % x y z -> x y z

drawshape

/fn 1 def

/model { neg exch } def % x y z -> x -z y

morerandom? { randconfig } if
drawshape

/fn 2 def

/model { neg 3 1 roll exch } def % x y z -> -z y x

morerandom? { randconfig } if
drawshape
} def

/rbool { rand 2 mod 0 eq } def

/randconfig {

cage? {
/vistest load 0 get /gt eq { %back
/fill? true def
/wire? rbool def
}{ %front
/fill? false def
/wire? true def
} ifelse
}{ %at least one of wire?,fill?

true rbool
rbool { exch } if /fill? exch def /wire? exch def

} ifelse

[{} {} {} {} {} {0 setgray} {0 setgray}
{1 setgray} {1 setgray} {1 setgray}
{normal light dot 1 add 4 div setgray}
{normal light dot 1 add 4 div 1 exch sub setgray}
{ 1 currentgray sub setgray }
{ currenthsbcolor exch pop .5 exch sethsbcolor }
{ currenthsbcolor 3 2 roll .33 add dup 1 gt {1 sub} if
3 1 roll exch pop .5 exch sethsbcolor }
{ currenthsbcolor 1 exch sub sethsbcolor }
{ currentrgbcolor exch 1 exch sub exch setrgbcolor }
{ currenthsbcolor setrgbcolor }
{ currentrgbcolor sethsbcolor }
{ currentrgbcolor 3 1 roll exch setrgbcolor }
//colorface //colorface //colorface //colorface //colorface
] dup
dup length rand exch mod get /wirecolor exch def
dup length rand exch mod get /fillcolor exch def

currentlinewidth .2 rand 6 mod 3 div add mul setlinewidth

/fudge 1 rand 10 mod 100 div add def
} def

/morerandom? true def

/cage? true def

/axes { %neg:black/white pos:white/black
gsave currentlinewidth 2 mul setlinewidth
-2 0 0 proj moveto
0 0 0 proj lineto
gsave currentlinewidth 2 mul setlinewidth stroke grestore
0 setgray currentpoint stroke moveto
2 0 0 proj lineto
gsave currentlinewidth 2 mul setlinewidth stroke grestore
1 setgray stroke

0 -2 0 proj moveto
0 0 0 proj lineto
gsave currentlinewidth 2 mul setlinewidth stroke grestore
0 setgray currentpoint stroke moveto
0 2 0 proj lineto
gsave currentlinewidth 2 mul setlinewidth stroke grestore
1 setgray stroke

0 0 -2 proj moveto
0 0 0 proj lineto
gsave currentlinewidth 2 mul setlinewidth stroke grestore
0 setgray currentpoint stroke moveto
0 0 2 proj lineto
gsave currentlinewidth 2 mul setlinewidth stroke grestore
1 setgray stroke
grestore

} def

/ang 0 def
{
/OW ang roty def

/ang ang 3 add def

setuppage

%/fill? true def
%/wire? true def
%/wirecolor { 1 currentgray sub setgray } def
back
randconfig
drawsolid

axes

[luser.droog]

Here's a silly little scene with a cube, a tetrahedron, and axis lines.
But now, the object stays still and the camera moves around it.

I suppose next I should try to learn how to do Z-buffering.

%!
/olddiv /div load def
/div { dup 0 eq { pop pop 100000 }{ olddiv } ifelse } def
(mat.ps)run

/I3 3 ident def

/disp <<
/cam [ 0 0 10 ] % Camera position
/theta [ 0 0 0 ] % Rotation sequence
/eye [ 0 0 20 ] % Eye relative to image surface
/Rot I3
>> def

/makerot {
theta 0 get roty
theta 1 get rotx matmul
theta 2 get rotz matmul
} def

% Ax Ay Az
/proj { DICT begin
%3 array astore
%dup == flush
cam {sub}vop %Camera translation
%pstack()=
Rot matmul %Camera rotation
0 get aload pop % Dx Dy Dz
eye aload pop % Dx Dy Dz Ex Ey Ez
%pstack()=
4 3 roll div % Dx Dy Ex Ey Ez/Dz
exch neg % Dx Dy Ex Ez/Dz -Ey
4 3 roll add % Dx Ex Ez/Dz Dy-Ey
1 index mul % Dx Ex Ez/Dz Ez(Dy-Ey)/Dz
4 1 roll 3 1 roll % Ez(Dy-Ey)/Dz Ez/Dz Dx Ex
sub mul exch % Ez(Dx-Ex)/Dz Ez(Dy-Ey)/Dz
%pstack ()=
end } dup 0 disp put def

/L 10 def
/axes {
[
[ [ 0 0 L neg ] [ 0 0 L ] ]
[ [ 0 L neg 0 ] [ 0 L 0 ] ]
[ [ L neg 0 0 ] [ L 0 0 ] ]
] { aload pop proj moveto proj lineto } forall
} def

/v [[ 1 1 -1 ] %cube vertices
[ -1 1 -1 ]
[ -1 -1 -1 ]
[ 1 -1 -1 ]
[ 1 1 1 ]
[ -1 1 1 ]
[ -1 -1 1 ]
[ 1 -1 1 ]] def
/fv [[ 0 1 2 3 ] %cube faces out of vertices
[ 0 4 5 1 ]
[ 1 5 6 2 ]
[ 2 6 7 3 ]
[ 3 7 4 0 ]
[ 4 7 6 5 ]
] def
/ev [[0 1][1 2][2 3][3 0] %cube edges
[0 4][1 5][2 6][3 7]
[4 5][5 6][6 7][7 4]] def
/tv [[ 0 5 7 ] %tetrahedron faces
[ 0 5 2 ]
[ 2 5 7 ]
[ 0 2 7 ]] def

/R 20 def
/H -3 def
/ang 0 def
%stepon
{
%300 400 translate
300 700 translate
1 70 dup dup scale div setlinewidth
%(1)=

{
disp begin
/cam [
ang sin R mul
H
ang cos R mul %neg
] def
/theta [
ang %neg
%180 add
H R atan %neg
0
] def
%2 copy get ang add put
/Rot makerot def
end
} exec%pop

tv {
{ v exch get proj } forall
moveto lineto lineto closepath
1 0 0 setrgbcolor fill
} forall

tv {
{ v exch get proj } forall
moveto lineto lineto closepath
0 setgray stroke
} forall

{
ev {
{ v exch get proj } forall
moveto lineto %lineto lineto closepath
} forall
0 1 0 setrgbcolor stroke
} exec%pop

axes 0 0 1 setrgbcolor stroke
/flushpage where { pop flushpage
80000 { .2 sin .2 cos atan pop } repeat
} if
showpage

/ang ang 3 add def

} loop

currentfile flushfile
{
[ [1 1 1] [2 2 2] [3 3 3] ]
dup { {add}vop } vop % yay! vop composes!
pstack
%:r!gsnd -q %
[[2 2 2] [4 4 4] [6 6 6]]
0 0 0 proj moveto
1 0 0 proj lineto
1 1 0 proj lineto
0 1 0 proj lineto
closepath
0 0 1 proj lineto
1 0 1 proj lineto
1 0 0 proj lineto
1 0 1 proj moveto
1 1 1 proj lineto
1 1 0 proj lineto
1 1 1 proj moveto
0 1 1 proj lineto
0 1 0 proj lineto
0 1 1 proj moveto
0 0 1 proj lineto
} pop

[luser.droog]

Finally, a rotating view from a little above without skew distortions.
It's less "imposing", but more accurate.

For some reason I had to reverse the visibility test again (gt <=> lt).


549(0)09:27 AM:ps 0> cat 3d5c.ps

%!
%Steinmetz, take 5: NEEDFORSPEED
% O(N) generator makes caching largely irrelevant!
% Bounded faces makes clipping unnecessary
% (which wasn't quite working somehow)
% Improved visibility now shows the Front of the object!
% Eliminated those ugly artifacts
% Random Configurations (with optional "less-random" mode)

% Revolving Camera removes skew distortions

4 2 90 steinmetz-gen

/E [ 3 3 10 ] def %eye point %Replaced by disp/cam get

/crackE { % set pointers into E
/ex E 0 1 getinterval cvx def
/ey E 1 1 getinterval cvx def
/ez E 2 1 getinterval cvx def
} def crackE
/E^ E dup [ exch mag neg dup dup ] { div } vop def %eye center vec

/vistest { lt } def
%Check visibility from the eye
/visible { % [ v1 v2 v3 v4 ] . bool
dup
aload pop /v4 exch def /v3 exch def /v2 exch def /v1 exch def
%dup 0 get 1 index 2 get 2 index 3 get % [] v1 v2 v4
%2 index { sub } vop % [] v1 v2 v14
%3 1 roll exch { sub } vop % [] v14 v12
%cross
v1 v4{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if
v1 v2{sub}vop
dup mag 0 eq { pop v1 v3{sub}vop } if
cross
/normal exch def %normal ==
%/normal [ normal mag dup dup ] { div } vop def
dup 0 get

%E
disp /cam get

{ sub } vop % [] ve1
/ev exch def %ev ==
ev normal dot %dup =
0 vistest % [] bool
exch pop
} def

/light

%[ -7 12 30 ] def
[ 3 2 7 ] def

/colorface {
%normal light dot 1 add 3 div setgray
%face 0 get 90 rotz matmul 0 get
%aload pop pop %exch
%dup 0 eq { pop .001 } if
%atan 360 div 3 div
fn .33 mul %add
%.5
face 0 get -90 rotz matmul 0 get
aload pop pop exch
dup 0 eq { pop .001 } if
atan 360 div 2 div .25 add
normal light dot 1 add 3 div
3 { 3 1 roll dup 1 gt { pop 1 } if dup 0 lt { pop 0 } if } repeat
sethsbcolor
} def

{
/proj { DICT begin /z exch def /y exch def /x exch def
1 ez z sub div
x ez mul z ex mul sub
1 index mul
y ez mul z ey mul sub
3 2 roll mul
end } dup 0 10 dict put bind def

} pop

/disp <<
/cam [ 0 0 10 ] % Camera position
/theta [ 0 0 0 ] % Rotation sequence
/eye [ 0 0 20 ] % Eye relative to image surface
/Rot I3

/Rad 50 def
/Ht 5 def

>> def

/makerot {
theta 0 get roty
theta 1 get rotx matmul
theta 2 get rotz matmul
} def

% Ax Ay Az
/proj { DICT begin

3 array astore
%dup == flush
cam {sub}vop %Camera translation
%pstack()=
Rot matmul %Camera rotation
0 get aload pop % Dx Dy Dz
eye aload pop % Dx Dy Dz Ex Ey Ez
%pstack()=
4 3 roll div % Dx Dy Ex Ey Ez/Dz
exch neg % Dx Dy Ex Ez/Dz -Ey
4 3 roll add % Dx Ex Ez/Dz Dy-Ey
1 index mul % Dx Ex Ez/Dz Ez(Dy-Ey)/Dz
4 1 roll 3 1 roll % Ez(Dy-Ey)/Dz Ez/Dz Dx Ex
sub mul exch % Ez(Dx-Ex)/Dz Ez(Dy-Ey)/Dz
%pstack ()=

end } dup 0 disp put bind def


/revcam {
disp begin
/ang exch def
/cam [
ang neg sin Rad mul
Ht
ang neg cos Rad mul %neg
] def
/theta [
ang neg
%180 add
Ht Rad atan %neg

0
] def
%2 copy get ang add put
/Rot makerot def
end

} def


/setuppage {
310 -600 translate
1 250 dup dup scale div 2 div setlinewidth

2 setlinejoin
} def

/usecache? false def
usecache? { % Generate Data File to cache the faces
/data (stein.fac) def
{
/in data (r) file def
} stopped {
pop pop pop%clear
/out data (w) file def
/outbuf 128 string def
genfaces
flush out closefile
/out null def
/in data (r) file def
} if
/reopen { in closefile /in data (r) file def } def
/drawshape { in cvx exec reopen } def
}{ % Generate fresh data for each shape
/drawshape { genfaces } def
} ifelse

/up 90 rotx def
/left 90 roty def

/front { /vistest { gt } def } def
/back { /vistest { lt } def } def

/drawsolid {
/fn 0 def
/model {} def % x y z -> x y z
drawshape

/fn 1 def
/model { neg exch } def % x y z -> x -z y
morerandom? { randconfig } if
drawshape

/fn 2 def
/model { neg 3 1 roll exch } def % x y z -> -z y x
morerandom? { randconfig } if
drawshape
} def

/rbool { rand 2 mod 0 eq } def

/randconfig {
cage? {

/vistest load 0 get /lt eq { %back

/fill? true def
/wire? rbool def
}{ %front
/fill? false def
/wire? true def
} ifelse
}{ %at least one of wire?,fill?
true rbool
rbool { exch } if /fill? exch def /wire? exch def
} ifelse

[{} {} {} {} {} %{0 setgray} {0 setgray}
%{1 setgray} {1 setgray} {1 setgray}

{normal light dot 1 add 4 div setgray}
{normal light dot 1 add 4 div 1 exch sub setgray}
{ 1 currentgray sub setgray }
{ currenthsbcolor exch pop .5 exch sethsbcolor }
{ currenthsbcolor 3 2 roll .33 add dup 1 gt {1 sub} if
3 1 roll exch pop .5 exch sethsbcolor }
{ currenthsbcolor 1 exch sub sethsbcolor }
{ currentrgbcolor exch 1 exch sub exch setrgbcolor }
{ currenthsbcolor setrgbcolor }
{ currentrgbcolor sethsbcolor }
{ currentrgbcolor 3 1 roll exch setrgbcolor }
//colorface //colorface //colorface //colorface //colorface
] dup
dup length rand exch mod get /wirecolor exch def
dup length rand exch mod get /fillcolor exch def
currentlinewidth .2 rand 6 mod 3 div add mul setlinewidth
/fudge 1 rand 10 mod 100 div add def
} def

/morerandom? false def

%/randconfig {} /rndcfg /randconfig load def def
%rndcfg
%1536830211 srand
1887638 srand
/ang 0 def
{
%/OW ang roty def
ang revcam

/ang ang 3 add def

rrand =

setuppage

%/fill? true def
%/wire? true def
%/wirecolor { 1 currentgray sub setgray } def
back
randconfig
drawsolid

axes

%/fill? false def
%/wire? true def
%/wirecolor {} def
%2 setlinewidth
%currentlinewidth 2 mul setlinewidth
front
randconfig
drawsolid

pstack flush %cya

%disp === [ v1 dup == aload pop proj ] ==

%pause for applause
500000 { 100 dup cos exch sin exch atan pop } repeat
showpage

%exit
} loop