It is currently Thu, 27-07-17, 18:53 GMT

All times are UTC




Post new topic Reply to topic  [ 8 posts ] 
Author Message
PostPosted: Mon, 09-11-09, 21:55 GMT 
Offline
Site Admin
User avatar

Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4467
Location: Hamburg, Germany
Hi friends,

=========EDIT ========== EDIT =========EDIT================

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,
t00fri wrote:
Finally, let me note that within my algorithm, the multiplicative spherical correction factor (spherecorr in nms) precisely cancels 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.

=========EDIT ========== EDIT =========EDIT================


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.

http://forum.celestialmatters.org/viewtopic.php?t=339

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!

Image

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.

Image

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.

Image

+++++++++++++++++++++++++++++
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

Code:
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:

Code:
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:

Code:
// work out the vector of int divisors of texture width: width

int* divVec = new int[height];    
divVec[0] = 1;
int k = 1;
for (int m = height; m > 1; m--)
{
    if (width%m == 0)
    {         
        divVec[k] = width/m;
        k++;         
     }
}   


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
Code:
 = 1.0 / sin(PI * (y + 0.5) / (double)height);

in terms of the above integer divisors (divVec):

Code:
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;
            }   
       }   
       return bs;
}


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:

Code:
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];   
    int n;
   
    if (width%binSize == 0)
        n = width/binSize ;
    else
    {      
        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);
    }    
    else
    {   
        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);
             }   
        }          
    }
    return row;
}


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)
Code:
 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,

Code:
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!!

Code:
if (x < 2 * bs)
    dx = (double)(h[y][x] - h[y][x + bs]) / (double)bs;
else   
    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...

Fridger


Last edited by t00fri on Tue, 17-11-09, 17:01 GMT, edited 21 times in total.

Top
 Profile  
 
 Post subject:
PostPosted: Mon, 09-11-09, 22:55 GMT 
Offline
Site Admin
User avatar

Joined: Thu, 30-08-07, 22:52 GMT
Posts: 2726
Location: France, South, not far from Montpellier
Too much maths for me (as usual! :roll: ) but I strongly approve the process here!!

Soon your tools will be really outstanding, king of flawless gems! :wink: :o


Top
 Profile  
 
 Post subject:
PostPosted: Mon, 09-11-09, 23:38 GMT 
Offline
Site Admin
User avatar

Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4467
Location: Hamburg, Germany
;-)

Cheers,
Fridger

_________________
Image


Top
 Profile  
 
 Post subject:
PostPosted: Tue, 10-11-09, 1:34 GMT 
Offline
User avatar

Joined: Wed, 05-09-07, 0:09 GMT
Posts: 57
Location: Seattle, WA, USA
Fridger,

First of all, there's a huge improvement at the poles--very nice work there.

But, I think that the binning has some bad side effects in areas a few degrees away from the poles. The gradients seem unrealistically diminished outside of the regions with artifacts. Ideally, the pinch-corrected map should match the uncorrected one very closely outside the latitudes where artifacts appear. So maybe the filtering should only be applied in problem areas?

--Chris


Top
 Profile  
 
 Post subject:
PostPosted: Tue, 10-11-09, 19:58 GMT 
Offline
Site Admin
User avatar

Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4467
Location: Hamburg, Germany
chris wrote:
Fridger,

First of all, there's a huge improvement at the poles--very nice work there.

But, I think that the binning has some bad side effects in areas a few degrees away from the poles. The gradients seem unrealistically diminished outside of the regions with artifacts. Ideally, the pinch-corrected map should match the uncorrected one very closely outside the latitudes where artifacts appear. So maybe the filtering should only be applied in problem areas?

--Chris


Chris,

thanks for your constructive critique.

First off, let me emphasize that I am not doing some empirical fiddling here that can be switched on/off at will. Once we are ready to accept that we want to keep the effective resolution in x constant for all values of y, then there is comparatively little freedom left. Clearly, with this working hypothesis, we should not be loosing any physical information towards the poles since in-between we are working with optimal resolution.
By construction, the reduction of resolution towards the poles follows pretty well the theoretical expression for a spherical geometry.

Here is a first plot, illustrating in dark gray, the regime where in my algorithm I am recovering automatically the binSize = 1 as before.
Image

Moreover, here is a comparison of the exact formula for the binSize (following from requiring constant resolution for all y) (blue curve):

Code:
binSize_exact = 1/sin(Pi*(y+0.5)/height) >= 1


with the best approximation of binSize in terms of the (13) integer divisors of the cylindrical texture width (red histogram).

Image

As you can see, the approximation is pretty good throughout. It's perhaps seen even better in a log-log plot of the left-hand section with the North pole on the left (y=1):

Image

+++++++++++++++++++++++++++++++++
So, where do we have some freedom??
+++++++++++++++++++++++++++++++++

A normalmap displays the RGB-coded normal vectors, as we all know. Ideally, this vector involves the TRUE partial derivative and NOT the naive Newton's difference quotient
Code:
 dx = (f(y, x) - f(y, x+h)) / h

that only tends towards the derivative for h-> 0!

Once h is generalized towards playing the role of the binSize h <=> bs >=1
clearly that naive formula should be improved. That's what I did by using the very good 5-point expression as explained above.

So, OF COURSE, this is an improvement that in principle is independent from my rebinning algorithm. The 5-point gradient formula will also modify the normalmap (towards the better) for ALL values of the binSize, including binSize=1. Everyone with experience in numerical math is aware that numerical differentiation is delicate and simple first order differences are unreliable...

I thought I explained this in detail above and thus do not quite understand why you were wondering about some global modification of my normalmap compared to the old formula.

This is the reason why the red component in my above normalmap comparison
http://forum.celestialmatters.org/userp ... acts_1.jpg seems to diminish for all latitudes. But so what. What's the reference for your above claim?
I printed out and compared the various expressions (2-point, 5-point) approximating the x-gradient, and surely there are substantial differences. Why did you think that the simple 2-point formula is the right reference?


Finally, let me note that within my algorithm, the multiplicative spherical correction factor (spherecorr in nms) precisely cancels against the denominator bs in the x-derivative formula (dx in nms). Without my generalization away from bs = 1, this does not happen. ;-)


Fridger


Top
 Profile  
 
 Post subject:
PostPosted: Tue, 10-11-09, 22:48 GMT 
Offline
User avatar

Joined: Wed, 05-09-07, 0:09 GMT
Posts: 57
Location: Seattle, WA, USA
t00fri wrote:
This is the reason why the red component in my above normalmap comparison
http://forum.celestialmatters.org/userp ... acts_1.jpg seems to diminish for all latitudes. But so what. What's the reference for your above claim?
I printed out and compared the various expressions (2-point, 5-point) approximating the x-gradient, and surely there are substantial differences. Why did you think that the simple 2-point formula is the right reference?

I don't think that simple differencing is an ideal way to estimate the gradient, and your wider filter may be superior. But something--either the new derivative estimator or the binning (more likely)--is changing the output in a matter that systematically suppresses the gradient in x. It's not just that it's different, but that it's always biased toward zero. The goal is to remove artifacts, not to produce a low-pass filter, unless you can provide some reason that the original data requires this.

--Chris


Top
 Profile  
 
 Post subject:
PostPosted: Tue, 10-11-09, 23:59 GMT 
Offline
Site Admin
User avatar

Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4467
Location: Hamburg, Germany
chris wrote:
t00fri wrote:
This is the reason why the red component in my above normalmap comparison
http://forum.celestialmatters.org/userp ... acts_1.jpg seems to diminish for all latitudes. But so what. What's the reference for your above claim?
I printed out and compared the various expressions (2-point, 5-point) approximating the x-gradient, and surely there are substantial differences. Why did you think that the simple 2-point formula is the right reference?

I don't think that simple differencing is an ideal way to estimate the gradient, and your wider filter may be superior. But something--either the new derivative estimator or the binning (more likely)--is changing the output in a matter that systematically suppresses the gradient in x. It's not just that it's different, but that it's always biased toward zero. The goal is to remove artifacts, not to produce a low-pass filter, unless you can provide some reason that the original data requires this.

--Chris


Can you quantify? What exactly do you find disturbing? Have you got an image to show?
I should need some concrete image aspect that I can hang on to...in order to find out.
In any case, I know by direct comparison, that the x-gradient is lowered by replacing the 2-point difference with the 5-point formula. I can make a respective plot. But since 5-point differences should be way better, a flatter x-gradient does not seem to be bad by itself. The corresponding Celestia displays are also nice.

Fridger

Addendum:

As you probably know, the error of the 5-point formula wrto the true derivative f '(x) of a function f(x) is

Code:
O(f^(5)(x) * h^4 / 30)

with f^(5)(x) being the 5-th derivative of f(x) (!!). This means in particular
that my improved approximation calculates the derivative exactly for ANY polynomial surface up to and including degree 4, irrespective of the used value of h! Notably, identifying h=binSize, it implies that one is not falsifying texture information up to degree 4 polynomial surfaces by using a large binSize near the poles.

++++++++++++++++++
When you plot instead the corresponding naive 2-point approximant for larger values of h you will get complete nonsense, certainly not anything reminding you of the desired derivative of a 4-degree polynomial surface....
++++++++++++++++++

If you like to see some analytical example (MAPLE), let me know.

Incidentally, there is also a corresponding 5-point formula for 2D surfaces (x,y). That seems the ideal tool for effecting really good approximations for normal vectors...


Top
 Profile  
 
 Post subject:
PostPosted: Wed, 11-11-09, 20:14 GMT 
Offline
Site Admin
User avatar

Joined: Fri, 31-08-07, 7:01 GMT
Posts: 4467
Location: Hamburg, Germany
Chris,

here is a simple example, illustrating the VAST superiority of my 5-point approximation for the x-derivative over the naive 2-point form that you also used in your original nm16 code...

Take this function f(x) that might be imagined to describe a section in x-direction through a crater on the surface of Mars, say:

Image

As I stated in my mail above, the 5-point approximation produces the exact result for the x-derivative, since that "crater function f(x)" is a polynomial of degree 4. The result is also independent of the stepsize h.

Here is what one gets:

Image

++++++++++++++++++
The red curve is the exact derivative of my above "crater function" f(x), while the blue dashed-dotted curve is the 5-point approximation. Indeed, both curves are identical (as anticipated) ;-)

On the other hand the simple 2-point approximation gives a nonsensical result for h=10. This is the green curve.

This example should make clear that the larger binSize values required for the no-pinch correction algorithm do call for a replacement of the naive 2-point difference by 5-point, 7-point, or better approximations. All these formulae are simple enough to be fast improvements. Of course then the normalmaps will get modified globally, since before they were calculated with inappropriate difference approximations for the gradients.
++++++++++++++++++

Fridger

PS: The next step of improvements comes when one considers elevation data that are noisy. Forming numerical derivatives from noisy data is a classical subject of numerical math that I remember well from my respective University courses decades ago ;-) . But I suppose, nowadays one just clicks some Web pages and copies the formulae.


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 8 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 1 guest


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
Powered by phpBB® Forum Software © phpBB Group