# Plasma Fractals Coming Out Well, Just With Visible Border Lines Along Subdivisions

This topic is 3725 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi, I've been working on generating plasma fractals for about a day now and the images are coming out SO close to being perfect, only there are visible lines in the image where you can see each rectangular subdivision. So on one side of the border is a darker color where the other side is rather bright. I would say except for the visible borders, the images are perfect plasmas. BTW, I am using the Diamond-Square method. Any ideas off the top of anyone's head? Here's an example image:

##### Share on other sites
looks like an incorrect implementation of diamond square. If done correctly there shouldn't be any borders.

-me

##### Share on other sites
Some of the code is a bit sloppy (converting from pixel position vectors to vertex position vectors), sorry. [grin]

void CPlasmaFractal::Process(int numIterations /* = -1 */){	m_MaxIterations = (numIterations == -1 ? m_MaxIterations : numIterations);	m_VertexCount = Width() * Height();	m_CurrentIteration = 1;		m_Values.Resize(Width(), Height(), false); 	for (DWORD y = 0; y < Height(); y++) 		for (DWORD x = 0; x < Width(); x++) 			m_Values(x, y) = -0.5;	//Seed initial corners 	m_Values(0,0)						= 0.5; 	m_Values(Width() - 1, 0)			= 0.5; 	m_Values(0, Height() - 1)			= 0.5; 	m_Values(Width() - 1, Height() - 1)	= 0.5;	Generate(m_P1, m_P2, m_P3, m_P4, eCorner_Center);}//////////////////////////////////////////////////////////////////////////void CPlasmaFractal::Generate(Vec2 p1, Vec2 p2, Vec2 p3, Vec2 p4, ECorner corner){// 	if (Distance<float>(p1, p2) < 1.0 || Distance<float>(p2, p3) < 1.0 || // 		Distance<float>(p1, p4) < 1.0 || Distance<float>(p3, p4) < 1.0)// 		return;	//Do the center first	Vec2 mid = Vec2(p1.Midpoint(p2).X, p1.Midpoint(p4).Y);	m_Roughness					*= m_RoughnessMultiplier;	m_RandomRangeLow			*= m_Roughness;	m_RandomRangeHigh			*= m_Roughness;	double range = Random(m_RandomRangeHigh, m_RandomRangeLow);	//ASSERT(range != 0.0);	GetValueByVector(mid) = GetMidPointValue(p1, p2, p3, p4) + range;	//	//Now the other points	//	//store temp variables for the parameters that will be altered within subdivisions	float	roughness = m_Roughness;	float	tempupper = m_RandomRangeHigh,			templower = m_RandomRangeLow;	//find the top value from the top left and top right	Vec2 midT = p1.Midpoint(p2);	double val = GetMidPointValue(p1, p2, mid/*, Get4thPoint(midT, eCorner_TopSide)*/);	val += Random(m_RandomRangeHigh, m_RandomRangeLow);	//ASSERT(val != 0.0);	GetValueByVector(midT)	= val;	m_Roughness = roughness;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;	//find the right value from the top right and bottom right	Vec2 midR = p2.Midpoint(p3);	val = GetMidPointValue(p2, p3, mid/*, Get4thPoint(midR, eCorner_RightSide)*/);	val += Random(m_RandomRangeHigh, m_RandomRangeLow);	GetValueByVector(midR)	= val;	m_Roughness = roughness;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;	//find the bottom value from the bottom left and bottom right	Vec2 midB = p4.Midpoint(p3);	val = GetMidPointValue(p4, p3, mid/*, Get4thPoint(midB, eCorner_BottomSide)*/);	val += Random(m_RandomRangeHigh, m_RandomRangeLow);	GetValueByVector(midB)	= val;	m_Roughness = roughness;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;	//find the left value from the top left and bottom left	Vec2 midL = p1.Midpoint(p4);	val = GetMidPointValue(p1, p4, mid/*, Get4thPoint(midL, eCorner_LeftSide)*/);	val += Random(m_RandomRangeHigh, m_RandomRangeLow);	GetValueByVector(midL)	= val;	m_Roughness = roughness;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;	if (m_CurrentIteration >= m_MaxIterations)		return;	int tempiteration = ++m_CurrentIteration;	//	//Now subdivide	//	//(Comment one of the corners out to get the Sierpinski triangle)	//TopLeft	m_Roughness = roughness;	m_CurrentIteration = tempiteration;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;	Generate(p1, midT, mid, midL, eCorner_TopLeft);	//TopRight	m_Roughness = roughness;	m_CurrentIteration = tempiteration;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;	Generate(midT, p2, midR, mid, eCorner_TopRight);	//BottomRight	m_Roughness = roughness;	m_CurrentIteration = tempiteration;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;		Generate(mid, midR, p3, midB, eCorner_BottomRight);	//BottomLeft	m_Roughness = roughness;	m_CurrentIteration = tempiteration;	m_RandomRangeHigh = tempupper;	m_RandomRangeLow = templower;	Generate(midL, mid, midB, p4, eCorner_BottomLeft);}//////////////////////////////////////////////////////////////////////////double CPlasmaFractal::AverageSquare(double c1, double c2, double c3, double c4){	return (c1 + c2 + c3 + c4) / 4.0;}//////////////////////////////////////////////////////////////////////////double CPlasmaFractal::AverageTriangle(double c1, double c2, double c3){	return (c1 + c2 + c3) / 3.0;}//////////////////////////////////////////////////////////////////////////double CPlasmaFractal::GetMidPointValue(Vec2 p1, Vec2 p2, Vec2 p3, Vec2 p4){	return AverageSquare(	GetValueByVector(p1),							GetValueByVector(p2),							GetValueByVector(p3),							GetValueByVector(p4));}//////////////////////////////////////////////////////////////////////////double CPlasmaFractal::GetMidPointValue(Vec2 p1, Vec2 p2, Vec2 p3){	return AverageTriangle(	GetValueByVector(p1),							GetValueByVector(p2),							GetValueByVector(p3));}//////////////////////////////////////////////////////////////////////////double CPlasmaFractal::GetMidPointValue(Vec2 p1, Vec2 p2){	double val = (GetValueByVector(p1) + GetValueByVector(p2)) / 2;	//Decrease roughness// 	m_Roughness *= m_RoughnessMultiplier;// // 	m_RandomRangeLow	*= m_Roughness;// 	m_RandomRangeHigh	*= m_Roughness;	double range = Random(m_RandomRangeHigh, m_RandomRangeLow);	////////////////////////////////////////////////////////////////	////	////	TODO!!!!:	Find a way so that different iterations	////				have the same vertices all across (i.e.	////				the 10th iteration looks like the 3rd,	////				just with the inbetweens different	////////////////////////////////////////////////////////////////	return range + val;}//////////////////////////////////////////////////////////////////////////Vec2 CPlasmaFractal::Get4thPoint(Vec2 mid, ECorner corner){	Vec2 coord = GetVertexCoordByVector(mid);	switch (corner)	{		case eCorner_TopSide:			coord.Y -= 1.0f;			break;		case eCorner_RightSide:			coord.X += 1.0f;			break;		case eCorner_BottomSide:			coord.Y += 1.0f;			break;		case eCorner_LeftSide:			coord.X -= 1.0f;			break;	}		//Subtracting a negative is the same as adding a positive	if (coord.X < 0)		coord.X = (Width() - 1) + coord.X;	else if (coord.X >= Width())		coord.X = Width() - coord.X;	if (coord.Y < 0)		coord.Y = (Height() - 1) + coord.Y;	else if (coord.Y >= Height())		coord.Y = Height() - coord.Y;	coord = GetVectorFromVertexCoord(coord);	return coord;}//In the header fileVec2	GetVertexCoordByVector(Vec2 p){	float	xp = p.X / m_P3.X * double(Width() - 1),		yp = p.Y / m_P3.Y * double(Height() - 1);	return Vec2(xp, yp);}Vec2	GetVectorFromVertexCoord(Vec2 p){	float	xp = p.X * m_P3.X / double(Width() - 1),		yp = p.Y * m_P3.Y / double(Height() - 1);	return Vec2(xp, yp);}double& GetValueByVector(Vec2 p){	Vec2 v = GetVertexCoordByVector(p);	return m_Values(v.X, v.Y);}

Thank you so much for the info!

##### Share on other sites
I'm pretty sure I had this exact problem. I solved it by keep track of which pixels had been written to already and not writing to them again. I'm pretty sure that's not how the algorithm's supposed to work but it worked great for me.

*Actually I looked back at my code and it seems I didn't do that at all. I simply rewrote it to do the diamond square algorithm iterively and it all worked out.

##### Share on other sites
Hmm, I've been trying to avoid the iterative process. If all else fails, I guess I need to fix it somehow, so iterative might be the way. [cool]

But in the meantime, my stubborn self is trying to do it recursivly. [grin]

I've been thinking that if I figure out all four quadrants AND THEN subdivide, that will save me from the uninitialized values. But I'm not sure.

Any ideas on that?

(BTW, thanks for ALL the help so far. Insight is always good.)

##### Share on other sites
Quote:
 Original post by MountainCodeMonkeyBut in the meantime, my stubborn self is trying to do it recursivly.

Recursive is the "right" way. I'll have to dig around and find my DS algorithm implementation... not enough time to dig through your code atm, sorry

-me

##### Share on other sites
for what it's worth. here's my working version:

I haven't looked at it for a while, but I found that I had to make the bitmap +1 in both width and height to get an odd number of pixels along each axis. Otherwise the algorithm didn't work out correctly. That's what the patchupData method is all about.

generate is the thing that kicks off the generation.

doDiamond and doSquare should be self explanatory.

my class sits on top of a bitmap class, IIRC so there may be calls to base class methods not shown here.

Now you're very close with yours and will be better off just debugging yours. perhaps you just need that +1 to width and height? I don't know. maybe something in here can help shed light on your problem

Again sorry I don't have to time to help you fix this more directly =)

void MidpointDisplacementGenerator::patchupData(){	if ( m_dataHolder )	{		//lose the +1 width and +1 height		for ( int h = 0; h < m_settingsHolder.m_height; ++h )		{			memcpy( &m_dataHolder[h*m_settingsHolder.m_width*4], &m_data[h*m_settings.m_width*4], sizeof(unsigned char) * m_settingsHolder.m_width * 4 );		}		//move pointers and data around so we're all happy		unsigned char *tmp = m_dataHolder;		m_dataHolder = m_data;		m_data = tmp;		m_settings = m_settingsHolder;	}}void MidpointDisplacementGenerator::generate( std::list<MapGenerator::InterestPoint> &interestPoints ){	if ( m_settings.m_height*m_settings.m_width > 0 && m_settings.m_height == m_settings.m_width )	{		if ( m_data && m_settings.m_width%2 == 0 )		{			m_dataHolder = m_data;			m_settings.m_height += 1;			m_settings.m_width += 1;			//it's ok because m_dataHolder is holding the old pointer			m_data = new unsigned char[m_settings.m_height*m_settings.m_width*4];			memset( m_data, 0, sizeof(unsigned char)*m_settings.m_height*m_settings.m_width*4 );		}		assert( m_settings.m_width%2 == 1);		assert(m_data);		setDataAt(0, 0, rand()%256);		setDataAt(0, m_settings.m_height-1, rand()%256);		setDataAt(m_settings.m_width-1, 0, rand()%256);		setDataAt(m_settings.m_width-1, m_settings.m_height-1, rand()%256);		doPass(0);		patchupData();	}}void MidpointDisplacementGenerator::doPass( int pass ){	if ( pass < m_numPasses )	{		int stride = (int)( ((double)m_settings.m_width-1) / pow(2,pass));		for ( int h = 0; h < m_settings.m_height; h += stride )		{			for ( int w = 0; w < m_settings.m_width; w += stride )			{				doDiamond( MidpointDisplacementGenerator::Point(w,h), stride, pass );			}		}		for ( int h = 0, i = 0; h < m_settings.m_height; h += stride/2, ++i )		{			for ( int w = (i%2 == 0)? stride/2: 0; w < m_settings.m_width; w += stride )			{				doSquare( MidpointDisplacementGenerator::Point(w,h), stride/2, pass );			}		}		doPass(pass+1);	}}void MidpointDisplacementGenerator::doDiamond( const MidpointDisplacementGenerator::Point &point, int width, int pass ){	int avg = 0;		if ( unsigned char *dataAt = getDataAt( point.x, point.y ) )	{		avg += dataAt[0];		if ( dataAt = getDataAt( point.x, point.y + width ) )		{			avg += dataAt[0];			if ( dataAt = getDataAt( point.x + width, point.y ) )			{				avg += dataAt[0];				if ( dataAt = getDataAt( point.x + width, point.y + width ) )				{					avg += dataAt[0];					avg = (int)( (float)avg * 0.25f );					int lo = (int)(m_loRand/pow(2,pass));					int hi = (int)(m_hiRand/pow(2,pass));					int height = avg + getGaussianRand(lo, hi);					height = std::min( height, m_settings.m_hi );					height = std::max( height, m_settings.m_lo );					assert( height >= 0 && height < 256 );					setDataAt( point.x + width/2, point.y + width/2, (unsigned char)height );				}			}		}	}}void MidpointDisplacementGenerator::doSquare( const MidpointDisplacementGenerator::Point &point, int width, int pass ){	int avg = 0;	int count = 0;	unsigned char *dataAt = getDataAt( point.x - width, point.y );	if ( dataAt )	{		avg += dataAt[0];		count++;	}	dataAt = getDataAt( point.x + width, point.y );	if ( dataAt )	{		avg += dataAt[0];		count++;	}	dataAt = getDataAt( point.x, point.y - width );	if ( dataAt )	{		avg += dataAt[0];		count++;	}	dataAt = getDataAt( point.x, point.y + width );	if ( dataAt )	{		avg += dataAt[0];		count++;	}	assert(count > 0);	avg /= count;	int lo = (int)(m_loRand/pow(2,pass));	int hi = (int)(m_hiRand/pow(2,pass));	int height = avg + getGaussianRand(lo, hi);	height = std::min( height, m_settings.m_hi );	height = std::max( height, m_settings.m_lo );	assert( height >= 0 && height < 256 );	setDataAt( point.x, point.y, (unsigned char)height );}int MidpointDisplacementGenerator::getGaussianRand(int lo, int hi){	if ( lo < hi )	{		int range = hi - lo;		int val = (rand() % range + lo);		val += (rand() % range + lo);		val += (rand() % range + lo);		val += (rand() % range + lo);		float fVal = (float)val * 0.25f;		return (int) fVal;	}	else	{		return 0;	}}

##### Share on other sites
All I know is that effect you have their would make a really cool photo filter.

##### Share on other sites
It's too much code for me to check, especially since I don't have the framework to test it. If I had to make a bet, it would be that you are being inconsistent about your vertex convention. You are probably interpreting the 4 vertices as

1 2
4 3

in some places and as

1 2
3 4

in others.

If you go through the code carefully with the debugger (or even by hand), you'll probably find the problem.

##### Share on other sites
Quote:
Original post by Palidine
Quote:
 Original post by MountainCodeMonkeyBut in the meantime, my stubborn self is trying to do it recursivly.

Recursive is the "right" way. I'll have to dig around and find my DS algorithm implementation... not enough time to dig through your code atm, sorry

-me

maybe...the diamond square algorithm actually says to subdivide faces...but unless you want arbitrary dynamic viewing depth, this would be a pretty inefficient way to implement it...better to make a grid of your target resolution and then just recursively consider points on that grid that are closer together as being neighbors

##### Share on other sites
Well, technically, the bitmap isn't something I render to. I have an array of doubles (sized m_VertexCount) that store values in the range [0.0, 1.0] during generation and then when the ChildView of my app calls its OnPaint() function, it runs through the array by the vertices and colors the corresponding pixel accordingly. So technically, I have one pixel for every vertex and currently the image/vertex resolution is 513x513 (powers of 2 plus 1 in general).

I'll look through the given code samples here and see if I can't figure out what's going on with it.

Thanks a lot for all the help people!!

(Plus, I might be able to use this somewhere down the line in making Lego Universe, so you might be someone who contributed to a large game!! [grin])

##### Share on other sites
OK, so just to test my diamond-square algorithm with NO randomization yielded interesting results. The blend is pretty smooth, yet you can see the borders ever so slightly (or blatently if these things stand out to you like they do me).

I can see that my edges are pretty erroneous. I've tried both the vertex on the opposite side (wrapping around basically) and I've tried just cutting out that 4th vertex and averaging the triangle instead of a quad; same results really.

So, here are the first seven iterations of the regular blending, no randomization and the top-left corner seeded to 0, the top-right and bottom-left corners seeded to 0.5 and the bottom-right corner seeded to 1.0.

(These are sized down to 25%, click the images for the full size)

Iteration 1:

Iteration 2:

Iteration 3:

Iteration 4:

Iteration 5:

Iteration 6:

Iteration 7:

Does this point out any obvious errors in my code that veterans of the algorithm recognize or know about?

##### Share on other sites
I've had a look at my code properly this time it turned out to be recursive after all.

rather than call the function again for every new square generated when you subdivide. I processed all the squares at once then the function called itself with a higher level of detail so the function processes four times as many squares as the last one.

It's kinda hard to explain but it works. I've posted the code below. Should be fairly self explanatory and wholely unoptimised. I was written a quite a while ago.

void diamond_square(int mapsize,int detail,float range,float rough){	// mapsize and detail must be a power of 2. eg 256	float ave;	float r;	int midx;	int midy;	int x;	int y;	int step;			step=mapsize/detail;		for(x=0;x<mapsize;x=x+step)	{		for(y=0;y<mapsize;y=y+step)		{			// Compute midpoint of square			midx=x+step/2;			midy=y+step/2;			ave=(map[x][y]+map[x+step][y]+map[x][y+step]+map[x+step][y+step])/4;			r=randnum(-1*range,range);			map[midx][midy]=ave+r;									// Compute the four diamond mid-points			// left			midx=x;			midy=y+step/2;			ave=(map[x][y]+map[x][y+step]+map[x+step/2][y+step/2])/3;			if(x-step/2>=0)			{				ave=(ave+map[x-step/2][y+step/2])/2;			}			r=randnum(-1*range,range);			map[midx][midy]=ave+r;						// top			midx=x+step/2;			midy=y;			ave=(map[x][y]+map[x+step][y]+map[x+step/2][y+step/2])/3;			if(y-step/2>0)			{				ave=(ave+map[x+step/2][y-step/2])/2;			}			r=randnum(-1*range,range);			map[midx][midy]=ave+r;									// right			midx=x+step;			midy=y+step/2;			ave=(map[x+step][y]+map[x+step][y+step]+map[x+step/2][y+step/2])/3;			if(midx+step/2<mapsize)			{				ave=(ave+map[midx+step/2][midy])/2;			}			r=randnum(-1*range,range);			map[midx][midy]=ave+r;						// bottom			midx=x+step/2;			midy=y+step;			ave=(map[x][y+step]+map[x+step][y+step]+map[x+step/2][y+step/2])/3;			if(midy+step/2<mapsize)			{				ave=(ave+map[midx][midy+step/2])/2;			}			r=randnum(-1*range,range);			map[midx][midy]=ave+r;								}	}	if(step>3)	{		diamond_square(mapsize,detail*2,range*expf(-1*rough),rough);	}	}

If you recursively call each one of the four new squares created the it will process the same corner each time you end up calculating the new points based on the averages of points that have not been calculated yet. That's why you have to process everything at the same level of detail at once and also why this algorithm can quite naturally be implemented itertively. In fact, in my code, you could replace the call to itself with a for loop quite simply.

[Edited by - DaveG on November 1, 2007 10:16:59 AM]

##### Share on other sites
Hey man, thanks a lot! I haven't had much time the last couple days to work on it, but this sounds like a winner. I thought I had tried to do something along the lines of what you said, but maybe not quite.

I'll give it a try for sure though!

Thanks!

##### Share on other sites
DaveG's method is what I also used when I had this problem (though the implementation might differ). The problem is indeed that certain points are not initialized prior being used in the calculations with a depth-first approach. Doing level-by-level you'll be able to make sure that all points used to calculate a level are initialized.

##### Share on other sites

This topic is 3725 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.