Please ignore this post
, since I found that my no-pinch code for nms.cpp underlying the presented results contained a most stupid typo-based
bug. It's too bad, since apparently it has drastic negative consequences for my pinch correction algorithm presented in this post...
In the course of the following discussion, I even did particularly emphasize
that the correction factor ~ 1/sin() due to the spherical body geometry cancels against the no-pinch step-width in the x-gradient,
Finally, let me note that within my algorithm, the multiplicative spherical correction factor (spherecorr in nms)
against the denominator bs
in the x-derivative formula (dx in nms). Without my generalization away from bs = 1, this does not happen.
but unfortunately ... after the numerous copy and paste operations that
I did in the code, while trying many things out,
the stepsize bs has accidentally
remained in the denominator of the x-gradient formula, while indeed I
had cancelled by hand spherecorr against bs in the later line
towards the end of the nms code:
dx *= bumpheight * spherecorr; ===> dx *= bumpheight;
So this stupid bug almost means a fresh start for the whole algorithm.
Just stay tuned, please, there is already new light at the end of the tunnel
I'll replace this post VERY soon with new results!
The encouraging results below are of course true results, yet --due to this bug-- they did not emerge from a mathematically rigorous analysis. In the sense of an empirical
pinch correction, they may, however, retain some value.
The encouraging recent topography data from the Japanese lunar explorer mission KAGUYA have refocussed my interest to normalmaps and related issues. My texture tools for making highest quality normalmaps from scientific elevation data should be sufficiently known by people around here
A nasty and persistent problem in these ambitious endeavours has been the longstanding polar pinch effects
. In this thread, I want to report on an algorithm that I developed to suppress these effects strongly without sacrificing detail... This thread reports on work in progress. Let's see how far I got.
A major input were the lessons that I learned from tackling the analogous problem in case of optimized tesselations of (CMOD) shape models for non-spherical bodies, like Phobos.
Please note that this discussion is NOT meant to be a pedagogical presentation, but just some "raw" development blog.
For concreteness, let me focus here on the 8k normalmap for Mars
, as generated with my nms tool from the the MOLA-064 elevation data. Of course, I tested my algorithm for all available high-quality normalmaps (with similar results), but this would lead too far here.
Let me first flash a few "before and after" images
, to make you a bit curious, before explaining what I actually did...
The first image, shows the encouraging effect of my no-pinch algorithm around the North pole of Mars. To expose things as clearcut as possible, I displayed ONLY the normalmap with Celestia, assigning an "empty" base texture in the .ssc file (Texture ""). Hence, whatever you see in the following images is due to the normalmap!
This first "before and after" comparison is around the North pole of Mars
. On the left is the normalmap display from the current version of my nms tool. No pinch corrections implemented! The nasty radial pinch effects are clearly seen. On the right you see what my new algorithm did to the pinch! It's not perfect yet, but pretty much encouraging, I would say.
The next comparison shows the same for the South pole of Mars.
Note that the comparisons were done under exactly equal conditions!
Another way of displaying the improvement is via comparing the normalmaps (here close to the south-pole rim): The top image displays clearly a lot of vertical bar-type artefacts very close to the South pole rim, while my no-pinch algorithm produces a clean and smooth appearance instead, without sacrifice in resolution.
So what did I actually do??
( For now, I just added my modifications to the nms code, but this may change eventually.)
While textures for Celestia are typically represented as simple cylindrical projections
(aka equirectangular maps), these projections are not really natural in case of mostly spherical bodies. While the width of the texture is identified with the equatorial
circumference of the body, this leads to a massive discrepancy towards the polar regions. The simple cylindrical map retains always the same width for all latitudes, the actual circumference, however, shrinks proportional to sin(90 - latitude )
towards the poles (latitude = +-90 degrees). Since the cylindrical map retains always the same number of pixels for all latitudes, the resolution increases absurdly towards the poles and consequently produces the familiar radial polar pinch fluctuations...
As a first step, I therefore required that the polar and equatorial resolutions should remain the same. Since resolution is #pixels/length, this means
n0 / width = n / circumference (latitude) = n / (sin( 90 - latitude) * width)
Here n0 is the number of pixels in longitude of the cylindrical map, while n is the reduced number of pixels
that is imposed from retaining the same resolution for all latitudes
n = sin (90 - latitude) * n0 <= n0
So we just need to do a consistent pixel rebinning
to achieve this aim.
In the corresponding code, I started with a trivial method that finds all integer divisors of the texture width and stores it in an array, like so:
// work out the vector of int divisors of texture width: width
int* divVec = new int[height];
divVec = 1;
int k = 1;
for (int m = height; m > 1; m--)
if (width%m == 0)
divVec[k] = width/m;
Note this also applies for textures where width is not of power of two size
These integer divVec values are the possible candidate binSizes in the rebinning process.
Next I optimally approximated the new, exact binSizes in longitude for latitude y
= 1.0 / sin(PI * (y + 0.5) / (double)height);
in terms of the above integer divisors (divVec):
int binSize(int* divVec, int length, int height, int y)
// best approximation of spherecorr factor by
// integer divisors of texture size 'width'
double spherecorr = 1.0 / sin(PI * (y + 0.5)/(double)height);
int bs = 1;
double old = (double) (4 * height * height);
for (int l = 0; l < length; l++)
double delta = spherecorr - (double)divVec[l];
if (delta * delta < old)
bs = divVec[l];
old = delta * delta;
The next step consisted in a modified 2 byte read method of one row in latitude, such that an appropriate rebinning with binSize is implemented towards the poles:
short* readRowS16(FILE* in, int width, int binSize)
// Modified 16bit read method, with a rebinning of the simple
// cylindrical map with bin size 'binSize'.
// This is useful for avoiding pinch effects close to the poles.
short* row = new short[width];
if (width%binSize == 0)
n = width/binSize ;
cerr << "\nPinch correction error: width % binSize != 0\n";
cerr << "No Pinch correction applied\n";
binSize = 1;
if (binSize == 1)
for (int i = 0; i < width; i++)
row[i] = readS16(in);
for (int m = 0; m < n; m++)
double s = 0;
for (int iav = 0; iav < binSize; iav++)
// Average over 'binSize' pixels in each row
s += (double)readS16(in);
for (int iav = 0; iav < binSize; iav++)
// Substitute the average into the binSize pixels, which
// amounts to a rebinning
int i = m * binSize + iav;
row [i] = (short)(s/(double)binSize);
It should be self-explanatory.
This method identifies binSize pixels in every row with the pixel average over the binSize pixels , thus effectively enlarging the pixel width, thereby reducing the resolution in x as desired.
Finally, we must also modify the formulae for the gradients in x (longitude), taking into account the larger binSize towards the poles.
Previously, I have used a simple two-point approximation for the (partial) x-derivative
(Newton's difference quotient that tends towards the derivative for h-> 0)
dx = (f(y, x) - f(y, x+h)) / h , with the step width h=1. (i.e. pixel width =1)
Since now the step-width bs needs to be much bigger towards the poles,
dx = (f(y,x) - f(y,x+bs) ) / bs; bs >(>) 1
I actually used a far better approximation (5 point
) for the x-derivative that should work well for much larger values of the step-width = binSize bs!!
if (x < 2 * bs)
dx = (double)(h[y][x] - h[y][x + bs]) / (double)bs;
dx = -(double)(-h[y][(x + 2 * bs) % width] + 8 * h[y][(x + bs)% width] - 8 * h[y][x - bs] + h[y][x - 2 * bs]) / (double)(12 * bs);
Still the new nms code is significantly faster than the original one, since much fewer points are worked upon in the polar regime...
So far so good...
Let me know about anything I have overlooked or some further improvements, of course...