Hans-Bernhard Broeker wrote:
>As threatened, I've begun regression-testing prospective new
>patches against the existing code.
>
>And already found the first one: the current edition of the "with pixels &
>friends" patch (after fixing a conflict in misc.c...) causes off-by-one
>differences to the current head of CVS. Affected plots include a
>considerable subset of 'surface.dem' (starting at the first that has x**3
>in it, then every plot up to the second plot of the sinc function), the
>see-through sphere plot in world.dem, and all the pm3d "flush" option's
>demos.
>
I found and verified the change. I recall now I needed to retain
precision for the map3d_xy() routine beyond unsigned int for a
hyperplane model of corner pixel coordinates in the image routine. So I
modified the map3d_xy() routine in util3d.c in the following way:
==========
/* Function to map from user 3D space to normalized 'camera' view
* space, and from there directly to terminal coordinates */
void
map3d_xy(
double x, double y, double z,
unsigned int *xt, unsigned int *yt)
{
#ifdef WITH_IMAGE
double xtd, ytd;
map3d_xy_double(x, y, z, &xtd, &ytd);
*xt = (unsigned int) xtd;
*yt = (unsigned int) ytd;
}
void
map3d_xy_double(
double x, double y, double z,
double *xt, double *yt)
{
#endif
int i, j;
double v[4], res[4], /* Homogeneous coords. vectors. */
w = trans_mat[3][3];
v[0] = map_x3d(x); /* Normalize object space to -1..1 */
v[1] = map_y3d(y);
v[2] = map_z3d(z);
v[3] = 1.0;
for (i = 0; i < 2; i++) { /* Dont use the third axes (z). */
res[i] = trans_mat[3][i]; /* Initiate it with the weight factor */
for (j = 0; j < 3; j++)
res[i] += v[j] * trans_mat[j][i];
}
for (i = 0; i < 3; i++)
w += v[i] * trans_mat[i][3];
if (w == 0)
w = 1e-5;
#ifdef WITH_IMAGE
*xt = ((res[0] * xscaler / w) + xmiddle);
*yt = ((res[1] * yscaler / w) + ymiddle);
#else
*xt = (unsigned int) ((res[0] * xscaler / w) + xmiddle);
*yt = (unsigned int) ((res[1] * yscaler / w) + ymiddle);
#endif
}
=========
That is, it is a routine, map3d_xy_double() that retains the precision
and then map3d_xy() uses that routine and makes the cast to an unsigned
int. I certainly didn't expect that to make a change. The difference
here is likely that the result of
((res[0] * xscaler / w) + xmiddle)
is not cast directly from the register contents, but instead is saved in
memory as a variable, retrieved and then cast. Something must happen in
that transition to result in a difference of 1.
What do you think? Certainly, duplicating this code is a no-no. The
core translation could be defined as a macro and included in two
separate functions, one where xt and yt are doubles, one where xt and yt
are unsigned ints... or as I've argued in the discussion initiated by
Ethan, just plain ints.
Dan
--
Dan Sebald
email: daniel DOT sebald AT ieee DOT org
URL: http://acer-access DOT com/~dsebald AT acer-access DOT com/
|