[PATCH] Better terrain normal calculation

[sanderd17]

The following patch contains a diff for better normal calculation. It's better in the following ways:

• Mathematically equal (go take a piece of paper and work it out ;)
• 4 square roots less, a number of native operations and comparisons
• one Query to the terrain less
• less code in general
• better precision due to the elimination of the square roots

[sanderd17]

Some aditional notes:

• border cases (within the if-statements in the original code) are handled by the CalcPosition(Fixed) method. There is no reason to check for them twice
• 46 lines of code less to be exact ;)

[sanderd17]

Oops, forgot that the vectors n0 - n3 aren't normalised when the mean is calculated. This skews the result.

So it's not mathematically equivalent. But I'd even say the old method is mathematically skewed (for not normalizing the vectors after taking the cross product).

[scroogie]

Imho you can take out even more if you respect the regularity of the Grid. (up - down) will lead to 1 component being 0, because they have equal X coordinate, and another component being 2*d where d is the Grid distance. The same goes for left-right. The 2*d factor can be left out, as you normalize anyway, so the normal should be something like (left.Z - right.Z; up.Z - down.Z; 2*d).normalize().

[sanderd17]

You're right (except for the actual formula, as Y is the up coordinate).

I've updated the patch. Although now it's a bit more code than the previous patch, there are less operations done.

[scroogie]

Yeah, I know, I always forget the convention 0ad uses. You could also calculate the height directly, eliminating the 8 multiplications for the X and Z coordinates, but I guess readability is more important than that. :)

Ah and btw, nice find, care to sum up the difference in op count? :)

[sanderd17]

Which multiplications?

And counting the ops, that's a bit much. But OK.

A cross product on itself is at least 12 operations (and 5 overflow checks that probably contain some additional operations). A normalize is a length calculation (5 operations + a sqrt) and 3 divisions. So in total, I have 4 cross products (48 ops) 4 normalizes (32 ops + 4 sqrts) + 6 additional ops saved.

And I used 3 ops extra.

So that gives a net profit of at least 83 ops and 4 sqrts per normal calculated.

[scroogie]

The op count was just out of curiosity, as it seemed a good gain. Sorry if that bugged you. :)

About the multiplications: As the normal is only dependent on the Y coordinate, you wouldn't need the X and Z coordinates (in calcposition). So you could replace the 4 calls to calcposition and 4 vectors with 4 floats and direct heightmap lookups. Something like this:

+       float lz,rz,uz,dz;
+       i = clamp(i, (ssize_t)0, m_MapSize-1);
+       j = clamp(j, (ssize_t)0, m_MapSize-1);
+       lz = (i>0) ? float(HEIGHT_SCALE * m_Heightmap[j*m_MapSize + (i-1)]) : float(HEIGHT_SCALE * m_Heightmap[j*m_MapSize]);
+       rz = (i<(m_MapSize-1)) ? float(HEIGHT_SCALE * m_Heightmap[j*m_MapSize + i + 1]) : float(HEIGHT_SCALE * m_Heightmap[j*m_MapSize + i]);
+       uz = (j>0) ? float(HEIGHT_SCALE * m_Heightmap[(j-1)*m_MapSize + i]) : float(HEIGHT_SCALE * m_Heightmap[j*m_MapSize+i]);
+       dz = (j<(m_MapSize-1)) ? float(HEIGHT_SCALE * m_Heightmap[(j+1)*m_MapSize + i]) : float(HEIGHT_SCALE * m_Heightmap[j*m_MapSize+i]);
+       normal.X = lz-rz;
+       normal.Y = 2.*float(TERRAIN_TILE_SIZE);
+       normal.Z = uz-dz;
+       normal.Normalize();  

But to be honest I don't know if such a patch would be accepted. I'm not a contributor, I just found this ticket coincidentally and thought my first suggestion would improve your patch.

[historic_bruno]

Replying to sanderd17:

Oops, forgot that the vectors n0 - n3 aren't normalised when the mean is calculated. This skews the result.

So it's not mathematically equivalent. But I'd even say the old method is mathematically skewed (for not normalizing the vectors after taking the cross product).

So the correct thing to do would be to normalize them to avoid the skewing?

The test_CalcNormal* tests pass with or without the normalization of n0-n3 (maybe it doesn't test cases where that matters, or the error is within the allowed 0.01 delta), but they fail with the attached patch.

[sanderd17]

Damn, that's right. Sorry, I made a mistake. My method seems to prefer terrain with a bigger slope, so it's not mathematically equivalent. I guess we'll leave it this way.

But now that I'm inspecting it a bit better, shouldn't CCmpTerrain->CalcNormal use a fixed version of the method Terrain->CalcExactNormal? As that method is less CPU intensive (also when ported to fixeds), and it seems to do a better job at calculating the normal (as it takes account of the triangle you're on). That method does almost what I did, a single sqrt, and a single cross product.

[historic_bruno]

I don't know if one is better or worse, they are different methods of calculating the normal, the same distinction exists for finding terrain height (see GetVertexGroundLevel vs. GetExactGroundLevel)

[sanderd17]

Well, the CalcNormal from CCmpTerrain gets as arguments coordinates in the world, from those coordinates, it gets the nearest vertex, from which it asks the normal. That normal takes 4 sqrts to calculate with the current method.

While it could also ask the normal on the triangle it's on (CalcExactNormal), which would only take some minor calculations and a single sqrt. And it looks more exact to me.

E.g. when a unit is near a cliff, the nearest vertex can be on the edge of the cliff, if you say the cliff is almost vertical, the normal will point at almost 45°, while the center of the unit might even be on flat terrain, and using the exactNormal calculation, the normal will indeed point up.

This could fix the "Horse legs sticking out" problem you mentioned in #1988. As the horse center point will normally not be on such steep terrain.

[wraitii]

Hadn't noticed that ticket…

Seems good to me but for using CalcPosition, which is about as inefficient as you can get here. Given that only height info is needed, I recommend forgetting all about the rest and using

ssize_t hi = clamp(i, (ssize_t)0, m_MapSize-1);
ssize_t hj = clamp(j, (ssize_t)0, m_MapSize-1);
float leftHeight = float(m_Heightmap[hj*m_MapSize + hi-1]);
float rightHeight = float(m_Heightmap[hj*m_MapSize + hi+1]);
float upHeight = float(m_Heightmap[(hj-1)*m_MapSize + hi]);
float downHeight = float(m_Heightmap[(hj+1)*m_MapSize + hi]);

and then changing whatever.Y to whateverHeight, but the rest is good.

This indeed calculates vertex normal and not surface normal, which is calculated by CalcExactNormal, which could also use some work I believe (mainly because there is a ton of redundancy by calling CalcExactNormal()).

[scroogie]

Hi wraitii, I did that in comment #9.

[wraitii]

Oh, right, sorry, had only read the patch. Your solution is better than mine because I forgot that this might try to reach wrong/Out of range values.

[wraitii]

Closing this as it is apparently invalid per Sanderd17's comment.