Sign in to follow this  
RonHiler

Problems with Mollweide Projection

Recommended Posts

Hello all, I'm having an issue with my Mollweide map projection, and I can't figure it out. Maybe a fresh set of eyes will help. I have drawn the projection to the screen, seeded it with a few points, then told it to draw points that were to the north and south of the seed points for a few iterations, in different colors. Here is what I got: As you can see, it's not what I would have expected. It seems okay along the edges, but as you get nearer the center, the lines get "snakey" rather than nicely curved as they ought to be. (there is also an issue at the south pole, but that's a different problem, and one I can deal with, so ignore that). It looks to me like some kind of rounding error, but I can't figure out where it's coming from in the code. Here is a part of the relevant routine:
//---------------------------------------------------------------------------
int WorldBuilderClass::FindMapPointIndexByDirection(int CurrentIndex, DIRECTION Direction)
	{
	int X, Y;           //current Mollweide coord
	int PointOnRow;

	//get the current X, Y coords
	X = WorldBuilderPoint[CurrentIndex].GetData(WBP_POINT_LOCATION_X);
	Y = WorldBuilderPoint[CurrentIndex].GetData(WBP_POINT_LOCATION_Y);

	if (Direction == NORTH)
		{
		PointOnRow = CurrentIndex - FirstPointInRow[Y];
		return FirstPointInRow[Y-1] + RoundFloatToInt((double)PointOnRow/(double)MapExtents[Y]*(double)MapExtents[Y-1]);
		}

	if (Direction == SOUTH)
		{
		PointOnRow = CurrentIndex - FirstPointInRow[Y];
		return FirstPointInRow[Y+1] + RoundFloatToInt((double)PointOnRow/(double)MapExtents[Y]*(double)MapExtents[Y+1]);
		}
	}





CurrentIndex is the point we are calculating from. It is simply an index into a vector of a sequential list of all the points that runs from 0 (at the north pole) to about 411,000 at the south pole. X and Y are the location of the point broken down into column and row, without regard to the row centering (thus, the first point of a row is always at X=0, regardless of how many points are in that row). FirstPointInRow[] is a vector that lists, for each row, the index of the first point. The first row (the north pole) has one point, and FirstPointInRow[0] = 0. The second row has 90 points, and FirstPointInRow[1] = 1, while FirstPointInRow[2] = 91, and so on. I have confirmed the values in this vector are correct. Finally, MapExtents[] is another vector which is a list of the number of points in each row (1, 90, etc). I have also confirmed the values in this row to my satisfaction. The return value is the index of the calculated point in the desired direction. Any ideas? I just don't see where the error is coming from. I think the equation is right. Oh, here is the routine for RoundFloatToInt(), which does exactly what it says:
//---------------------------------------------------------------------------
inline int WorldBuilderClass::RoundFloatToInt(double Value)
	{
	int result;

	//warning: truncating a negative float is undefined in C++!!!
	//         e.g. -1.7 might truncate to -1 or -2 depending on computer/compiler
	//         Therefore the somewhat odd equation for negative values to
	//         get around this little problem

    if (Value > 0)
	{
	if (Value < (INT_MAX - 1))
            {
            result = (int)(Value + 0.5);
            }
        else
            {
            result = INT_MAX;
            }
        }
    else if (Value < 0)
        {
        if (Value > (INT_MIN + 1))
            {
            result = -(int)(-Value + 0.5);	//convert to positive, truncate, then convert back to negative
            }
        else
            {
            result = INT_MIN;
            }
        }
    else
        {
        result = 0;
        }

    return result;
    }





Share this post


Link to post
Share on other sites
No one has any thoughts, huh?

Was my explanation too confusing? If so I can try to explain it beter. Or did I post in the wrong forum? I figured it was a math problem, but perhaps the Graphics forum would be a better place to ask?

Share this post


Link to post
Share on other sites
I don't think this is something many people will be familiar with, so you'll have to wait and see if someone has helpful information for you. Give it a little more time, and maybe someone will have some advice. But, perhaps my unproven, untested analysis and suggestion below will be helpful?

I'm not sure you have really explained your code very well. I don't see anything in the code you posted that resembles the Mollweide projection, so the projection must be contained in your calculation of FirstPointInRow/MapExtents, which isn't in your post. The outer shape of the projection area does look correct, however, so I'll believe your comment about them being correct. You don't say anything about WorldBuilderPoint. You don't really say how you are using the result from FindMapPointIndexByDirection to create your image. Are you telling us that CurrentIndex is simply an index into a linear framebuffer/pixelbuffer/bitmap image? I assume so. My assumptions are:

1) The "index" is an index into a linear framebuffer/pixelbuffer/bitmap and basically tells you which memory address to color, and therefore is relative to the start of the framebuffer rather than the start of the row.

2) CurrentIndex is the index for your seed point. And from assumption #1 is relative to the start of the frame buffer not the row.

3) If you're marching north or south from a seed point, you appropriately update CurrentIndex after each pixel is colored...outside of the code you listed. And...you're basically doing something like CurrentIndex += FindMapPointIndexByDirection(CurrentIndex, NORTH | SOUTH);

4) WorldBuilderPoint.GetData(...) returns the (x,y) *pixel* coordinate (integer values) within the frame buffer for the current index. I am blindly assuming that this gives you a correct pixel coordinate since I see no code.

So, if I got most of that right, and given that FirstPointInRow and MapExtents are correct, I think the mapping in FindMapPointIndexByDirection is quite wrong. I think PointOnRow is computed incorrectly. And I think you are not doing what you think you are doing on the "return *" lines. Try this instead for north:


if (Direction == NORTH)
{
// Find index into frame buffer that is at the beginning of the frame
// buffer row
IndexToStartOfFrameBufferRow = (Y-1)*frame_buffer_width;

// compute the actual point relative to the start of the row
PointOnRow = CurrentIndex - FirstPointInRow[Y] - IndexToStartOfFrameBufferRow;

// update IndexToStartOfFrameBufferRow to reflect moving up one row
IndexToStartOfFrameBufferRow -= frame_buffer_width;

// finally, compute and return updated CurrentIndex. Note that I
// moved MapExtents[Y-1] into the numerator.
return FirstPointInRow[Y-1] + RoundFloatToInt((double)MapExtents[Y-1]*(double)
PointOnRow/(double)MapExtents[Y]) + IndexToStartOfFrameBuffer);
}


And for south:


if (Direction == SOUTH)
{
// Find index into frame buffer that is at the beginning of the frame
// buffer row. Yep, still Y-1 here, and not Y+1!
IndexToStartOfFrameBufferRow = (Y-1)*frame_buffer_width;

// compute the actual point relative to the start of the row
PointOnRow = CurrentIndex - FirstPointInRow[Y] - IndexToStartOfFrameBufferRow;

// update IndexToStartOfFrameBufferRow to reflect moving down one row. Gotta
// add frame_buffer_width here instead of subtract, to move south a row.
IndexToStartOfFrameBufferRow += frame_buffer_width;

// finally, compute and return updated CurrentIndex. Note that I
// moved MapExtents[Y+1] into the numerator.
return FirstPointInRow[Y+1] + RoundFloatToInt((double)MapExtents[Y+1]*(double)
PointOnRow/(double)MapExtents[Y]) + IndexToStartOfFrameBuffer);
}


Here, frame_buffer_width is sizeof(pixel_data) times the number of pixels in each row of the framebuffer. This is of course a fixed # that does not change per row, e.g., it is always 1024*sizeof(pixel_data) if you have a framebuffer that is 1024 pixels wide. sizeof(pixel_data) depends on what your frame buffer stores. If it is RGB color only, and 1 byte (8 bits) per color component, then sizeof(pixel_data) would be 3*1 = 3 bytes or 24 bits. Depending on API, it might be padded out to whole words, e.g., a pixel would have to be 32 bits/4 bytes on a 32-bit operating system. Er, you're already writing out colors and displaying them so you'll realize all this I am absolutely certain!

Note that the return value in my code samples *is* the new CurrentIndex, so you would use it as:

// marching north:
CurrentIndex = FindMapPointIndexByDirection(CurrentIndex, North);

// marching south:
CurrentIndex = FindMapPointIndexByDirection(CurrentIndex, South);

My code here is untested and subject to my interpretation of what you're doing, but I think you should find this works out better. Let me know!

Share this post


Link to post
Share on other sites
Thanks Graham,

Sorry, I guess I didn't explain my code very well. It's hard sometimes to figure out what is relevant and what is not.

Unfortunatley, some of your assumptions are wrong, so your code fragments aren't right. I'll try to clarify a bit.

Quote:

1) The "index" is an index into a linear framebuffer/pixelbuffer/bitmap and basically tells you which memory address to color, and therefore is relative to the start of the framebuffer rather than the start of the row.

No. None of the code I have shown has anything to do with the drawing routines, it's a level above that. I have omitted the code controlling the drawing of the map, as I don't think it is relevant. So don't think in terms of pixel buffers or frame buffers.

Index is simply a reference into a set of points (actually a vector of a class that holds data about each individual point). The points go from 0 (which will be found at the north pole) to roughly 411,000 (which is found at the south pole). It is strictly sequential, there is no projection of any sort associated with these points.

The Mollweide projection, as you surmised, comes from the MapExtents vector, which tells the program how many of those points are in each row (so, since there is exactly 1 point in the first row, MapExtents[0] is 1. MapExtents[1] is 90, there being 90 points in the second row, and so on for all the rest of the rows.

The calculation for MapExtents is probably irrelevent, but just in case, it looks like this:

WorldBuilderRow = 0.0;
while (WorldBuilderRow < 2.0 * WORLDBUILDER_WORLD_RADIUS)
{
//Get the circumference of the horizontal circle at this lattitude
Circumference = 2.0 * PI * sqrt((2.0 * WORLDBUILDER_WORLD_RADIUS * WorldBuilderRow) - (WorldBuilderRow*WorldBuilderRow));
NumPointsInRow = RoundFloatToInt(Circumference);
if (NumPointsInRow == 0)
NumPointsInRow = 1;
MapExtents.push_back(NumPointsInRow);
WorldBuilderRow += (2.0/PI);
}


The radius (a constant) is set such that there ends up being exactly 512 rows, with the circumference at the equator as 1024 points. Anyway, that is where the projection comes from.

FirstPointInRow is simply a shortcut. It is the equivelent of:


Total = 0;
for (unsigned int i=0; i<MapExtents[CurrentRow]-1; i++)
{
Total += MapExtents[i];
}
FirstPointInRow.push_back(Total);


Where CurrentRow would be the row you are looking for the index of the first point of. It's just a shortcut so I don't have to keep running through that loop every time I want to know the index of the first point of a given row. Does that make more sense?

Quote:

2) CurrentIndex is the index for your seed point. And from assumption #1 is relative to the start of the frame buffer not the row.

The first sentence is right, the second is not. It is neither relative to the frame buffer nor the row. CurrentIndex is relative to the sequential set of points (in other words, it is simply a number that runs anywhere from 0 to 411,000).

Quote:

3) If you're marching north or south from a seed point, you appropriately update CurrentIndex after each pixel is colored...outside of the code you listed. And...you're basically doing something like CurrentIndex += FindMapPointIndexByDirection(CurrentIndex, NORTH | SOUTH);

More or less the right idea, although it's not a += function. What we want to return is the absolute index, they are not additive. What it comes out to is (and I'm paraphrasing my code here, as it's actually much more complex than this):

NeighborIndex = FindMapPointIndexByDirection(CurrentIndex, NORTH);

Once I have the NeighborIndex value, I can then set the color of the point.

Quote:

4) WorldBuilderPoint.GetData(...) returns the (x,y) *pixel* coordinate (integer values) within the frame buffer for the current index. I am blindly assuming that this gives you a correct pixel coordinate since I see no code.

No. Definitely NOT pixel coordinates, although there is a bit of a relationship. The centering of the points on the rows is done in the drawing routines. Here, the X and Y values are the column and row, just like pixel coordinates, but the X value is not centered. Therefore the first point (which would be WorldBuilderPoint[0]) would return 0,0 for X and Y. The second point (WorldBuilderPoint[1]) would return 0, 1 (this would be the first point on the second row), while WorldBuilderPoint[2] would return 1, 1.

I know, it's a horrendously complicated system, and I'm asking a lot for someone to come in and look at a fragment of code and try to see a problem that I myself can't see.

Maybe the problem DOES lie in the equation. The equation I came up with (in pseudo code] was:

NewIndex = FirstPointInRow[DesiredRow] + X/MapExtents[StartingRow]*MapExtents[DesiredRow]

Anyway, I hope my clarifications help. Like I said, it's kind of hard to know how deep to go into it when describing complex code like this.


Share this post


Link to post
Share on other sites
Well, for starters, I'm an idiot, because in the function FindMapPointIndexByDirection(), the values for X and PointOnRow are exactly the same. Basically I recalculated the value that I had already calculated earlier and stored in the point class. So I've taken out the calculation for PointOnRow and replaced it in the formula with X. Not that it made any difference other than to slightly speed up that routine.

Anyway, I WILL find the bug, and when I do I'll post a followup here on the off chance anyone ever comes along with the same problem. Not too likely, I suppose, but nevertheless :)

Share this post


Link to post
Share on other sites
Quote:
Original post by RonHiler
Thanks Graham,

Sorry, I guess I didn't explain my code very well. It's hard sometimes to figure out what is relevant and what is not.

Unfortunatley, some of your assumptions are wrong, so your code fragments aren't right. I'll try to clarify a bit.


This is awesome! I'm glad I made assumptions, since it did enable you to better describe your system, and keep the discussion rolling!

Share this post


Link to post
Share on other sites
I wonder if your whole problem is simply that (for north):


return FirstPointInRow[Y-1] + RoundFloatToInt((double)PointOnRow/(double)MapExtents[Y]*(double)MapExtents[Y-1]);


should be


return FirstPointInRow[Y-1] + RoundFloatToInt((double)MapExtents[Y-1]*(double)PointOnRow/(double)MapExtents[Y]);


(With a similar change made to the south code)

I'd try that simple mod to the code you have, and see what happens.

Your original image actually reminds me somewhat of so-called "interrupted" Mollweide projections. Not quite the same though...

Another thing you might try is a brute force approach using the formal definition of the Mollweide projection. It'll be slower, but could give you a theoretical point of comparison. The Wikipedia and Mathworld pages are decent:

Mollweide Projection @ Wikipedia
Mollweide Projection @ Mathworld
Some illustrations of interrupted Mollweide projections

Share this post


Link to post
Share on other sites
Quote:
Original post by grhodes_at_work
I wonder if your whole problem is simply that (for north):


return FirstPointInRow[Y-1] + RoundFloatToInt((double)PointOnRow/(double)MapExtents[Y]*(double)MapExtents[Y-1]);


should be


return FirstPointInRow[Y-1] + RoundFloatToInt((double)MapExtents[Y-1]*(double)PointOnRow/(double)MapExtents[Y]);


(With a similar change made to the south code)

I'd try that simple mod to the code you have, and see what happens.

Errr, those are exactly equivalent equations, unless I'm losing it :) I tried it for S&G, and got the same results (as far as I could tell). Unless there is some weird rounding thing going on related to the order of the operations?

Quote:
Another thing you might try is a brute force approach using the formal definition of the Mollweide projection. It'll be slower, but could give you a theoretical point of comparison. The Wikipedia and Mathworld pages are decent:

Good idea, thanks! I will give those equations a try if I can't garner any more information from the debugging I am doing. In fact, I could calculate them both ways and take a difference to see exactly where the divergences are happening. That should be pretty helpful. Awesome!

Share this post


Link to post
Share on other sites
Quote:
Original post by RonHiler
Errr, those are exactly equivalent equations, unless I'm losing it :) I tried it for S&G, and got the same results (as far as I could tell). Unless there is some weird rounding thing going on related to the order of the operations?


Oops! Guess you're right. I never memorized operator precedence, and try to avoid dealing with it by writing equations in strict left-to-right sequence and explicitly using parentheses. I find this easier to read, in any case.

How confident are you, really, about your final drawing routine? Oh, well, I'd love to hear how this progresses.

Share this post


Link to post
Share on other sites
Hey Graham,

I figured out the problem. The solution, however, is another matter :)

It's a bit hard to explain, so hopefully I can get the gist across. Basically, the "snakeyness" is a result of an accumulation of rounding errors.

As a test, I picked a starting point on the map and found 1) the number of points on that row and 2) the point along the row we were on, and then did a division. It turned out to be 0.4017 or something along those lines.

I then ONLY took the north points from that point (using the same routine we discussed before) and did the same calculation for each row. In theory that division number should *always* be the same (since it's just a ratio of the point on the row over the number of points in the row). As it turned out the division ranged anywhere from 0.3985 to 0.4100 or thereabouts (I don't have it open in front of me, but it was around there). That discounts the two or three rows just before the north pole, where the rounding got way out of whack, heh.

The reason it is doing that is there is no persistant state information. The algorithm is taking the ratio based on the current row to calcuate the point on the next row (which of course will not be exactly the same ratio, just as close as it can get). When it goes to the NEXT row up, it starts all over, using the slightly off ratio of the second row to calculate the third, which will be slightly off from the second row (and thus possibly even FURTHER off from the first row), and so on. Thus building up all the various rounding errors. Therefore the snaking rather than the nice curve.

Does that make sense?

The question is how to solve that problem?

I COULD calculate the ratio on the first row (whatever it happens to be) and then use that for each subsequent calcualation (same ratio for row 2, 3, 4, and so on) rather than using the current row we were on. Basically pass it as a parameter to the function (with a default value of -1 which would mean to calculate the ratio rather than use the parameter).

I dunno, I'm still bouncing ideas around in my head. I'm open to suggestions if you should have any.

Share this post


Link to post
Share on other sites
Hey Graham,

Sorry to necro this thread, but I did say that I'd follow up with the solution. While I haven't coded it yet, I have redesigned the algorithm. I don't know for sure if it will work, but I think it will.

This is actually the third attempt. The first one you saw. The second was even more complex and failed miserably. This one is simpler, and hopefully gets around the accumulated rounding error issue (I hope). See what you think and tell me if you agree.

What I plan to do is, for every individual point, calculate it's best match point on the equator (essentially turn the molleweide projection into a mercator projection one point at a time). I will do so by calculating a range that the point covers if it were to lie on the equator, like so:

Range1 = (PointInRow/NumPointsInRow) * NumPointsInEquator
Range2 = ((PointInRow+1)/NumPointsInRow) * NumPointsInEquator

That should give us the range of points on the equator that the test point would cover if the row it was on were stretched to the size of the equator row.

Then we just take an average of Range1 and Range2, convert it to an int, and subtract one (to go from 1 based to 0 based indexing).

EquatorMatchPoint = RoundFloatToInt((Range1 + Range2)/2.0) - 1

We store those values (one for each point on the map).

To find a neighbor point, we take the point we start on and look at its EquatorMatchPoint (EMP). Then we look at the points on the row above us (if we are looking for the North point) and find the closest match to that EMP, and that will be the neighbor. In some cases there will be multiple matches, and in such cases (I think) we simply take the centrist point of all those that match (if there are three matches, we'd take the second). [Alternately, I could instead store EMP as a double rather than an int, which would always give a single "best" match point, then convert that to an int to get the index point, that may be a better idea].

Since we are relating every point back to the equator row at the beginning, we should eliminate the issue of accumulated rounding errors we were getting before (when we were relating each point to the point on the prior row).

What do you think? Any suggestions? This is now going to be the second time I've ripped the code out and restarted it, and if I get through coding it again and it doesn't work, I may start to despair :)

I really sincerly appreciate your help with this issue, and your letting me bounce ideas off of you. You've been a great help to me!

Share this post


Link to post
Share on other sites
I didn't try to deeply understand your idea, and am kind of slammed at the moment, but in theory if you can compute everything relative to the same fixed reference, you should stand a better chance of minimizing the error, at least if you keep everything as int's. Floats, or doubles, could cause trouble, actually.

One thing I didn't think about before is...there is an open source library called GDAL for doing all sorts of geospatial transformations, and it supports (it seems) the Mollweide projection. You might take a look at that, if you are open to using someone else's code as a reference. I cannot promise GDAL will be totally intuitive at the start (do you feel the foreshadowing?):

GDAL

(The GDAL Warp API might prove useful, I guess...)

Share this post


Link to post
Share on other sites
I got it!

The algorithm above works, with a bit of modification (as written, it still gives a bit of drift [albiet not as much as before]. You have to take your ratio of the point you start from, find the closest ratio to that on the equator row, then apply THAT ratio to find the closest ratio on the row above or below the row you started from. A bit convoluted, but it works fine, and keeps everything nice and relative to the equator row).

These are random seed points, which are then drawn south to north only:



I still have a little bit of error in there, but it doesn't seem major. The basic algorithm appears to be sound. Probably a fencepost or possibly a rounding error. Nothing I can't track down and fix pretty readily.

Thanks again for the assist, Graham!

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this