Problems with Mollweide Projection

214

Author

July 07, 2008 05:35 AM

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;
    }

Creation is an act of sheer will

RonHiler

214

Author

July 07, 2008 10:59 AM

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?

Creation is an act of sheer will

grhodes_at_work

1,385

July 07, 2008 05:13 PM

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!

Graham Rhodes Moderator, Math & Physics forum @ gamedev.net

RonHiler

214

Author

July 07, 2008 06:36 PM

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;		}	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.

Creation is an act of sheer will

RonHiler

214

Author

July 08, 2008 08:11 AM

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

Creation is an act of sheer will

grhodes_at_work

1,385

July 08, 2008 09:26 AM

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!

Graham Rhodes Moderator, Math & Physics forum @ gamedev.net

grhodes_at_work

1,385

July 08, 2008 09:46 AM

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

Graham Rhodes Moderator, Math & Physics forum @ gamedev.net

RonHiler

214

Author

July 08, 2008 11:06 AM

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!

Creation is an act of sheer will

grhodes_at_work

1,385

July 08, 2008 11:25 AM

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.

Graham Rhodes Moderator, Math & Physics forum @ gamedev.net

RonHiler

214

Author

July 09, 2008 03:14 PM

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.

Creation is an act of sheer will

Problems with Mollweide Projection

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Problems with Mollweide Projection

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines