How to “trap” my surface patches to prevent the ba

In response to a challenge in comp.lang.postscript, I'm working-up my 3D chops trying to render a cylinder as projected rectangular patches. But I'm still seeing the wire-frame even after I comment-out the line-drawing, because the patches don't butt-up flush.

The cylinder is modeled along the z-axis by double-looping over z (-2 .. 2, step 4/N) and theta (0 .. 360, step 360/N). The four points of the rectangle are:

v1 = (Rcos T, Rsin T, z)
v4 = (Rcos T, Rsin T, z+dz)
v2 = (Rcos (T+dT), Rsin (T+dt), z)
v3 = (Rcos (T+dT), Rsin (T+dt), z+dz)

Then we apply a model->world rotation to all four points. Then we take the vectors v1->v4 and v1->v2 and do a cross product to get the normal vector for the patch. Take a dot product with the eye vector to check if the patch is on "this side" of the shape; if not, skip the drawing and procede to the next patch (fall off the bottom of the loop). Then we apply a perspective projection to each point and draw the quadrilateral with regular postscript 2D moveto and lineto. One last calculation on the normal vector to set the graylevel and then fill.

So the question is: Is there a usual way to deal with this? Is it a 3D problem or just a numerical problem (floating-point round-off kind of stuff)? Do I just add a little fudge-factor to my dz and dT when calculating the points? Or stroke the edges explicitly? These last 2 options both produce the desired result but I can't say that I'm satisfied with them. While each make be used on an individual illustration, it doesn't really solve the problem, you know?

I took a dump of the points being used. Here's the first few from N=12. It appears to me that, as predicted, v2 and v3 coincide precisely with v1 and v4 of the next piece on the band. These are the 2D "user coordinates" passed to moveto and lineto to produce the individual quadrilaterals. Since the CTM doesn't change, these points should map to the same pixels, right? So it does appear to be a very similar issue to the linked question. But I'm using Postscript precisely because I don't want to muck-about with writing my own rasterization routine :). I really think that the solution from the linked question, mapped to Postscript, would be to reverse the orientation of alternating checkerboard squares, at least for even N. That way, all corresponding edges are drawn in the same direction (as each other).

[-2.64550757 2.08465409]
[-3.00470281 1.69015563]
[-2.7090168 1.69015563]
[-2.38403082 2.08465409]

[-3.00470281 1.69015563]
[-3.28940701 0.936108589]
[-2.96660638 0.936108589]
[-2.7090168 1.69015563]

[-3.28940701 0.936108589]
[-3.4 -0.0666666701]
[-3.0666666 -0.0666666701]
[-2.96660638 0.936108589]

[-3.4 -0.0666666701]
[-3.28940701 -1.05890918]
[-2.96660638 -1.05890918]
[-3.0666666 -0.0666666701]

[-3.28940701 -1.05890918]
[-3.00470281 -1.78584146]
[-2.7090168 -1.78584146]
[-2.96660638 -1.05890918]

I've added a simple light model and tweaked it to bring out more mids. Jpeg output doesn't exhibit the problem, presumably due to the lossy compression. So here's a PNG snapshot.

The effect is much more apparent if I use the eye-vector as the light source. Here's xpost on the left showing the problem and gs on the right showing a modification where dz and dt are multiplied by a fudge factor of 1.06.

And the code: [Do not use this code. There is an error in the matmul routine. Corrected routines available here. Completed challenge available here.]

%!
%A shaded cylinder! Woohoo!

%(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 ]

end } def


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

% 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

% x y z -> x' y' z'
/model {
%ang roty
%ang .25 mul rotx
%alpha rotz
beta roty
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

showpage

Alright, I've come up with something that sits a little easier in the gut.

6% fudging just feels to horrible to bear.

But Ken suggested that rounding could be involved. That means taking control of the rounding should gain one some measure of control over the problem. And it looks like it's true.

So I tried prepending all moveto and lineto calls with a call to prep:

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

The comments show the various pieces I tried. Rounding on both device coordinates eliminated all horizontal bleed-lines and leaves very thin vertical bleeds. This seems to make sense assuming Ghostscript rasterizes by horizontal scanlines: it has an easier time with the horizontal ones with just a little help, but near-verticals are tougher.

But then I combined this with fudging. And I found that rounding just the device-y 'ordinate and fudging the patch dimensions by 2% eliminates all bleeds. It really lit up this batcave.

2% is an acceptable level of fudging, I think. (?)

Unfortunately, all the above requires tweaking when you adjust the value of N (the number of slices). The simplest fix to cover the whole surface is to stroke the edges in the same color as the fill. The only difficult point here is making sure the linewidth is appropriate for the scale. And the easy way to do that is to set them both together. For very high resolutions, this should probably be adjusted in some way to account for N.

1 70 dup dup scale div setlinewidth

Here's one of the images generated by the final program, a Steinmetz solid with coordinate axes and random colors, in a slightly skewed perspective (its right foot sticks out a little).