trueharmoniccolours.co.uk

The code I’ve already given for doing wavelet transforms is very simple but I treat horizontal and vertical passes as seperate functions. It occurred to me as I looked at those functions that they are, essentially, the same. The only real difference is how many pixels you step over to get to the next or previous one. This was a massive realisation to me. It meant that I could do all my wavelet transforms through a single function.

Here is the Foward transform function:

const int32_t postDwtNum	= dwtMax / 2;

const int32_t nextReadOffset2		= nextReadOffset * 2;
//const int32_t nextWriteOffset2		= nextWriteOffset * 2;

int32_t readOffset	= startOffset;

int32_t s			= readOffset;
int32_t d			= readOffset	+ (postDwtNum * nextWriteOffset);

const DataType d1	= m_SignalData[readOffset + nextReadOffset]	- (((m_SignalData[readOffset]	+ m_SignalData[readOffset + nextReadOffset2])	/ (DataType)2));
const DataType s1	= m_SignalData[readOffset]					+ (((d1							+ d1)											/ (DataType)4));

out.m_SignalData[d]	= d1;
out.m_SignalData[s]	= s1;

s	+= nextWriteOffset;
d	+= nextWriteOffset;

readOffset	+= nextReadOffset2;

int dwt = 2;
while( dwt < dwtMax - 2 )
{
	const DataType d1	= m_SignalData[readOffset + nextReadOffset]	- (((m_SignalData[readOffset]				+ m_SignalData[readOffset + nextReadOffset2])	/ (DataType)2));
	const DataType s1	= m_SignalData[readOffset]					+ (((out.m_SignalData[d - nextWriteOffset]	+ d1)											/ (DataType)4));

	out.m_SignalData[d]	= d1;
	out.m_SignalData[s]	= s1;

	s	+= nextWriteOffset;
	d	+= nextWriteOffset;

	readOffset	+= nextReadOffset2;

	dwt += 2;
}
{
	const DataType d1	= m_SignalData[readOffset + nextReadOffset]	- (((m_SignalData[readOffset]				+ m_SignalData[readOffset])	/ (DataType)2));
	const DataType s1	= m_SignalData[readOffset]					+ (((out.m_SignalData[d - nextWriteOffset]	+ d1)						/ (DataType)4));

	out.m_SignalData[d]	= d1;
	out.m_SignalData[s]	= s1;
}
return true;

What a relisation this was! This meant that I could now generalise the algorith to handle multi dimension signals. A 1D signal is pretty simple and is a single call to this function per wavelet "level"

int32_t level	= 0;
while( level < numLevels )
{
	Signal< DataType >::FDWT( 0, 1, 1, GetWidth() >> level, out );
	level++;
}

A 2D signal is a little more complicated as you need to do it for each horizontal line of an image and then each vertical line of an image. The code however still goes through that one function and looks something like this:

// Now run DWT.
int32_t level	= 0;
while( level < numLevels )
{	
	int32_t y		= 0;
	int32_t yMax	= GetHeight() >> level;
	while( y < yMax )
	{
		out.Signal< DataType >::FDWT( y * GetWidth(), 1, 1, GetWidth() >> level, out2 );			// Horizontals
		y++;
	}
	
	int32_t x		= 0;
	int32_t xMax	= GetWidth() >> level;
	while( x < xMax )
	{
		out2.Signal< DataType >::FDWT( x, GetWidth(), GetWidth(), (GetHeight() >> level), out );	// Verticals
		x++;
	}
	level++;
}

Now this then led me on to thinking about a 3D signal. To do the forward transform a 3D signal would now mean for each depth slice performing the 2D transform as above. After that you would then perform the transform a long a depth line or constant row and column. Thats a hell of a lot of calculations but I "think" (Though I have not tested this at all) that the code would look something like this:

int32_t level	= 0;
while( level < numLevels )
{
	int32_t z		= 0;
	int32_t zMax	= GetDepth();
	while( z < zMax )
	{
		int32_t y		= 0;
		int32_t yMax	= GetHeight();
		while( y < xMax )
		{
			out2.Signal< DataType >::FDWT( (y * GetWidth()) + (z * widthHeight), 1, 1, GetWidth() >> level, out );	// Horizontals
			x++;
		}

		int32_t x		= 0;
		int32_t xMax	= GetWidth();
		while( x < xMax )
		{
			out.Signal< DataType >::FDWT( x, GetWidth(), GetWidth(), GetHeight() >> level, out2 ); // Verticals
			x++;
		}
		z++;
	}

	int32_t y		= 0;
	int32_t yMax	= GetHeight();
	while( y < yMax )
	{
		int32_t x		= 0;
		int32_t xMax	= GetWidth();
		while( x < xMax )
		{
			out2.Signal< DataType >::FDWT( x + (y * GetWidth()), GetWidth() * GetHeight(),  GetWidth() * GetHeight(), GetDepth() >> level, out ); // "Depthicals"
			x++;
		}
		y++;
	}
	

	level++;
}

This 3D method gives a potential compression schemes for video. Though I'd hate to think the kind of horsepower you'd need to have available to process a video real-time.

Allin, I cannot believe I didn't spot this simplification before hand ...

So it turned out that the refinement pass was actually pretty easy to do.

Lets say I’m compressing a pixel with a value of -129 and I’m using a threshold of 128. I push into a vector the value |129| – 128. ie the absoloute value minus the threshold. I subsequently do this for every pixel above my threshold.

pixelsForRefinement.push_back( absVal - threshold );

I then halve the threshold and do the same sort of processing again. Now I run through the values that were pushed into my vector on the pass before hand and I do a simple check. Is the value in the vector (in the example case 1) greater than 128? No so I write a 0 to the refinement stream. If it IS greater than I write a 1 to the refinement stream and subtract the threshold value from the value stored in the vector.

// We now need to refine the pixels that were added on the previous pass
uint32_t refinement		= 0;
uint32_t refinementMax	= lastNumPixelsForRefinement;
while( refinement < refinementMax )
{
	const int16_t val				= pixelsForRefinement[refinement];
	const BitEncoding::Bit bit		= (val >= threshold) ? BitEncoding::Bit_One : BitEncoding::Bit_Zero;
	pixelsForRefinement[refinement]	-= (threshold * (uint32_t)bit);
	
	bitStreamOut.refineBitStream	+= bit;

	refinement++;
}
lastNumPixelsForRefinement	= pixelsForRefinement.size();

Effectively this means that for each pass through I refine the already WDR’d pixels by the current threshold.

The result is a perfect re-construction from the wavelet information (I shan’t bother posting up a before and after picture as they look the same ;)).

Now to get some form of compression I need to be able to early out from the algorithm and, basically, only store high frequency information (as well as the low frequency image in the first place). This posed a slight problem, however as I need to be able to perform the compression on each wavelet step as I go along. The system I had built thus far performed the compression on all the LH etc segments for each wavelet decomposition level. I needed to do some code re-organisation to get it do them in the correct order (ie LH, HL and HH). It wasn’t too difficult.

At this point though I realised I could MASSIVELY simplify the wavelet transforms. This I will explain in the next article.