waveletprocessing.cs

Image Fusion Techniues

// Waveblend - complex dualtree based image fusion
// (C) Copyright 2004 -- Sebastian Nowozin <[e-mail]>
//
// This file is part of Waveblend.
//
// Waveblend is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation; version 2 of the License.
//
// Waveblend is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// The license is included with the distribution in the file 'LICENSE'.
//

/* WaveletProcessing.cs
 *
 * A lot of this code is based on the excellent Matlab code by
 * Shihua Cai, Keyong Li, Farras Abdelnour and Professor Ivan Selesnick, which
 * is available at
 * http://taco.poly.edu/WaveletSoftware/index.html
 *
 * Thank you very much for this work, without which this program would not
 * have been possible.
 */

using System;

public class
WaveletProcessing : ICloneable
{
	public object Clone ()
	{
		return (new WaveletProcessing ());
	}

	public void Convolve1d (IWavelet wave, ref double[] inVal,
		ref double[] outTrend, ref double[] outFluctuation)
	{
		SimpleArray saIn = new SimpleArray (inVal);
		SimpleArray saTrend = new SimpleArray (outTrend);
		SimpleArray saFluctuation = new SimpleArray (outFluctuation);

		Convolve (wave, saIn, saTrend, saFluctuation);
	}

	private void ConvolveAFBWorking (IWavelet wave, IWaveletArray inVal,
		IWaveletArray outTrend, IWaveletArray outFluctuation)
	{
		double[] loFilt = wave.ScalingNumbersAnalysis;
		double[] hiFilt = wave.WaveletNumbersAnalysis;
		int valClen = loFilt.Length + inVal.Length - 1;

		double[] trend = new double[valClen];
		double[] fluct = new double[valClen];

		int N = inVal.Length;
		int L = loFilt.Length / 2;

		for (int k = 0 ; k < valClen ; ++k) {
			for (int j = 0 ; j < loFilt.Length ; ++j) {
				int vIdx = k - j;

				if (vIdx < 0 || vIdx >= inVal.Length)
					continue;

				trend[k] += loFilt[j] * inVal[vIdx];
				fluct[k] += hiFilt[j] * inVal[vIdx];
			}
		}

		// Debug output
		foreach (double val in trend)
			Console.WriteLine ("trend: {0}", val);
		Console.WriteLine ("");
		foreach (double val in fluct)
			Console.WriteLine ("fluct: {0}", val);
	}

	// Convolution
	//
	// Preconditions:
	// inVal.Length % 2 == 0
	// wave.SupportSize >= inVal.Length
	// outTrend.Length == outFluctuation.Length == inVal.Length / 2
	//
	private void Convolve (IWavelet wave, IWaveletArray inVal,
		IWaveletArray outTrend, IWaveletArray outFluctuation)
	{
		/*Console.WriteLine ("Convolve (\"{0}\", inVal[{1}], ..)",
			wave.Name, inVal.Length);*/

		double[] loFilt = wave.ScalingNumbersAnalysis;
		double[] hiFilt = wave.WaveletNumbersAnalysis;

		ConvolveArr (loFilt, hiFilt, inVal, outTrend, outFluctuation);
	}

	private void ConvolveArr (double[] loFilt, double[] hiFilt,
		IWaveletArray inVal,
		IWaveletArray outTrend, IWaveletArray outFluctuation)
	{
		int valClen = loFilt.Length + inVal.Length - 1;

		int N = inVal.Length;
		int L = loFilt.Length / 2;
		int ON = N / 2;

		for (int k = 0 ; k < ((valClen + 1) / 2) ; ++k) {
			for (int j = 0 ; j < loFilt.Length ; ++j) {
				int vIdx = 2*k - j;

				if (vIdx < 0 || vIdx >= inVal.Length)
					continue;

				vIdx += L;
				vIdx %= inVal.Length;

				if (k >= N)
					continue;

				outTrend[k % ON] += loFilt[j] * inVal[vIdx];
				outFluctuation[k % ON] += hiFilt[j] * inVal[vIdx];
			}
		}
	}


	public void Convolve2d (IWavelet wave, double[,] inVal,
		double[,] outTrend, double[,] outH,
		double[,] outD, double[,] outV)
	{
		Convolve2dArr (wave.ScalingNumbersAnalysis, wave.WaveletNumbersAnalysis,
			wave.ScalingNumbersAnalysis, wave.WaveletNumbersAnalysis,
			inVal, outTrend, outH, outD, outV);
	}

	public void Convolve2dArr (double[] loFiltCol, double[] hiFiltCol,
		double[] loFiltRow, double[] hiFiltRow,
		double[,] inVal,
		double[,] outTrend, double[,] outH,
		double[,] outD, double[,] outV)
	{
		if ((inVal.GetLength (0) % 2) != 0 ||
			(inVal.GetLength (1) % 2) != 0)
		{
			throw (new ArgumentException (String.Format
				("Input array has to be of even dimensions, got ({0}x{1}).",
				inVal.GetLength (0), inVal.GetLength (1))));
		}

		double[,] rowTrends =
			new double[inVal.GetLength (0), inVal.GetLength (1) / 2];
		double[,] rowFlucts =
			new double[inVal.GetLength (0), inVal.GetLength (1) / 2];

		// From F create T_{Row} and F_{Row}
		for (int row = 0 ; row < inVal.GetLength (0) ; ++row) {
			ConvolveArr (loFiltCol, hiFiltCol,
				new DirectionalArray (DirectionalArray.Direction.Row, inVal, row),
				new DirectionalArray (DirectionalArray.Direction.Row, rowTrends, row),
				new DirectionalArray (DirectionalArray.Direction.Row, rowFlucts, row));
		}

		// From T_{Row} create H and A (Trend)
		for (int col = 0 ; col < rowTrends.GetLength (1) ; ++col) {
			ConvolveArr (loFiltRow, hiFiltRow,
				new DirectionalArray (DirectionalArray.Direction.Column, rowTrends, col),
				new DirectionalArray (DirectionalArray.Direction.Column, outTrend, col),
				new DirectionalArray (DirectionalArray.Direction.Column, outH, col));
		}

		// From F_{Row} create V and D
		for (int col = 0 ; col < rowFlucts.GetLength (1) ; ++col) {
			ConvolveArr (loFiltRow, hiFiltRow,
				new DirectionalArray (DirectionalArray.Direction.Column, rowFlucts, col),
				new DirectionalArray (DirectionalArray.Direction.Column, outV, col),
				new DirectionalArray (DirectionalArray.Direction.Column, outD, col));
		}
	}

	//
	// 2D DT Complex Discrete Wavelet Transform
	//

	// Returns two trees, each which hold 3 complex subbands per level plus
	// the trend signal. That are six real signals per level.
	//
	// The returned array 'dt' has the following structure:
	//
	// dt[0] holds the first 3 subband directions
	// dt[1] holds the second 3 subband directions
	//
	// dt[0][0-2] are the three the ComplexArray H, D, V bands of plane 1
	// dt[1][0-2] are the other three the ComplexArray H, D, V bands of
	//   plane 2
	// dt[0].GetTrend (0/1) gets the first or second trend (lowpass subband)
	//   signals of plane 1
	// dt[1].GetTrend (0/1) gets the first or second trend (lowpass subband)
	//   signals of plane 2

	public Dualtree2dComplex[] DualtreeTransform2dComplex (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, double[,] inVal, int scale)
	{
		Dualtree2dComplex[] trees = new Dualtree2dComplex[2];
		Dualtree2dReal[] treesReal = new Dualtree2dReal[2];

		// Orientation index t2
		for (int t2 = 0 ; t2 < 2 ; ++t2) {
			// Real/imaginary index t1
			for (int t1 = 0 ; t1 < 2 ; ++t1) {
				/*Console.WriteLine ("DualtreeTansform2dComplex: creating real tree ({0}, {1})",
					t1 == 0 ? "Real" : "Imaginary", t2 == 0 ? "Orientation 1" : "Orientation 2");*/

				treesReal[t1] = DualtreeTransform2dRealSingleTree (wave1rst,
					wave, inVal, scale, t1, t2);
			}

			trees[t2] = new Dualtree2dComplex (treesReal[0], treesReal[1]);
		}

		// Now trees[0] is 1 orientation, trees[1] is 2 orientation, each
		// holding 3 subbands.

		Dualtree2dComplex tree1 = trees[0];
		Dualtree2dComplex tree2 = trees[1];

		while (tree1 != null) {
			for (int m = 0 ; m <= 2 ; ++m) {
				DualtreeTransform2dComplexSubtract (tree1[m], 0, tree2[m], 1);
				DualtreeTransform2dComplexSubtract (tree2[m], 0, tree1[m], 1);
			}

			tree1 = tree1.Next;
			tree2 = tree2.Next;
		}

		return (trees);
	}

	private void DualtreeTransform2dComplexSubtract
		(Dualtree2dComplex.ComplexArray c1, int ri1,
		Dualtree2dComplex.ComplexArray c2, int ri2)
	{
		/*Console.WriteLine ("DualtreeTransform2dComplexSubtract (c1, {0}, c2, {1})",
			ri1, ri2);*/
		double sqrt2 = Math.Sqrt (2.0);

		for (int y = 0 ; y < c1.GetLength (0) ; ++y) {
			for (int x = 0 ; x < c1.GetLength (1) ; ++x) {
				double u = (c1[y, x, ri1] + c2[y, x, ri2]) / sqrt2;

				c2[y, x, ri2] = (c1[y, x, ri1] - c2[y, x, ri2]) / sqrt2;
				c1[y, x, ri1] = u;
			}
		}
	}

	// Complex 2d DT INVERSION
	// Warning: this changes the original data stored in the tree.
	//
	public double[,] DualtreeTransform2dComplexInverse (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, Dualtree2dComplex[] trees)
	{
		Dualtree2dComplex tree1 = trees[0];
		Dualtree2dComplex tree2 = trees[1];

		while (tree1 != null) {
			for (int m = 0 ; m <= 2 ; ++m) {
				DualtreeTransform2dComplexSubtract (tree1[m], 0, tree2[m], 1);
				DualtreeTransform2dComplexSubtract (tree2[m], 0, tree1[m], 1);
			}

			tree1 = tree1.Next;
			tree2 = tree2.Next;
		}

		// The magic output array.
		double[,] result = new double[trees[0].GetTrend (0).GetLength (0) * 2,
			trees[0].GetTrend (0).GetLength (1) * 2];

		// Orientation index t2 (tree index)
		for (int t2 = 0 ; t2 < 2 ; ++t2) {
			// Real/imaginary index t1
			for (int t1 = 0 ; t1 < 2 ; ++t1) {
				double[,] trend = DualtreeTransform2dRealInverseSingleTree
					(wave1rst, wave, trees[t2].ProduceRealDualtree (t1),
					t1, t2, true);

				for (int y = 0 ; y < trend.GetLength (0) ; ++y)
					for (int x = 0 ; x < trend.GetLength (1) ; ++x)
						result[y, x] += trend[y, x];
			}
		}

		for (int y = 0 ; y < result.GetLength (0) ; ++y)
			for (int x = 0 ; x < result.GetLength (1) ; ++x)
				result[y, x] /= 4.0;

		return (result);
	}


	//
	// 2D DT Real Discrete Wavelet Transform
	//

	public Dualtree2dReal[] DualtreeTransform2dReal (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, double[,] inVal, int scale)
	{
		Dualtree2dReal[] trees = new Dualtree2dReal[2];

		for (int t = 0 ; t < 2 ; ++t)
			trees[t] = DualtreeTransform2dRealSingleTree (wave1rst,
				wave, inVal, scale, t);

		// Do the final add/subtraction and normalization step for each
		// fluctuation subsignal
		Dualtree2dReal tree1 = trees[0];
		Dualtree2dReal tree2 = trees[1];

		while (tree1 != null) {
			for (int m = 1 ; m <= 3 ; ++m)
				DualtreeTransform2dSubtract ((double[,]) tree1[m],
					(double[,]) tree2[m]);

			tree1 = tree1.Next;
			tree2 = tree2.Next;
		}

		return (trees);
	}

	private void DualtreeTransform2dSubtract (double[,] s1,
		double[,] s2)
	{
		double sqrt2 = Math.Sqrt (2.0);

		for (int y = 0 ; y < s1.GetLength (0) ; ++y) {
			for (int x = 0 ; x < s1.GetLength (1) ; ++x) {
				double u = (s1[y,x] + s2[y,x]) / sqrt2;

				s2[y,x] = (s1[y,x] - s2[y,x]) / sqrt2;
				s1[y,x] = u;
			}
		}
	}

	private Dualtree2dReal DualtreeTransform2dRealSingleTree (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, double[,] inVal, int scale, int tree)
	{
		return (DualtreeTransform2dRealSingleTree (wave1rst, wave, inVal,
			scale, tree, tree));
	}

	private Dualtree2dReal DualtreeTransform2dRealSingleTree (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, double[,] inVal, int scale, int t1, int t2)
	{
		Dualtree2dReal root = DualtreeTransform2dRealSub (wave1rst, t1, t2, inVal);
		Dualtree2dReal walk = root;

		for (int k = 1 ; k < scale ; ++k) {
			double[,] trend = (double[,]) walk[0];

			Dualtree2dReal subtree = DualtreeTransform2dRealSub (wave, t1, t2,
				trend);

			walk.Next = subtree;
			walk = subtree;
		}
		walk.LastLevel = true;

		return (root);
	}

	private Dualtree2dReal DualtreeTransform2dRealSub (IWaveletDTCW wave,
		int t1, int t2, double[,] inVal)
	{
		Dualtree2dReal dt = new Dualtree2dReal ();

		int yDim = inVal.GetLength (0) / 2;
		int xDim = inVal.GetLength (1) / 2;

		double[,] trend = new double[yDim, xDim];
		double[,] H = new double[yDim, xDim];
		double[,] D = new double[yDim, xDim];
		double[,] V = new double[yDim, xDim];

		Convolve2dArr (wave[FilterType.Analysis, t1, PassType.Trend],
			wave[FilterType.Analysis, t1, PassType.Fluctuation],
			wave[FilterType.Analysis, t2, PassType.Trend],
			wave[FilterType.Analysis, t2, PassType.Fluctuation],
			inVal, trend, H, D, V);

		dt[0] = trend;
		dt[2] = H;
		dt[3] = D;
		dt[1] = V;

		return (dt);
	}

	// Warning: this changes the original data stored in the tree.
	public double[,] DualtreeTransform2dRealInverse (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, Dualtree2dReal[] trees)
	{
		// Sum and difference (why?)
		Dualtree2dReal tree1 = trees[0];
		Dualtree2dReal tree2 = trees[1];

		while (tree1 != null) {
			for (int m = 1 ; m <= 3 ; ++m)
				DualtreeTransform2dSubtract ((double[,]) tree1[m],
					(double[,]) tree2[m]);

			tree1 = tree1.Next;
			tree2 = tree2.Next;
		}

		double[][,] outTrends = new double[2][,];

		for (int n = 0 ; n < 2 ; ++n)
			outTrends[n] = DualtreeTransform2dRealInverseSingleTree
				(wave1rst, wave, trees[n], n);

		double[,] outData = new double[outTrends[0].GetLength (0),
			outTrends[0].GetLength (1)];

		for (int y = 0 ; y < outTrends[0].GetLength (0) ; ++y) {
			for (int x = 0 ; x < outTrends[0].GetLength (1) ; ++x) {
				outData[y, x] = (outTrends[0][y, x] +
					outTrends[1][y, x]) / 2.0; // FIXME: Farras uses sqrt(2.0) here, why?
			}
		}

		return (outData);
	}

	private double[,] DualtreeTransform2dRealInverseSingleTree
		(IWaveletDTCW wave1rst, IWaveletDTCW wave, Dualtree2dReal root,
		int tree)
	{
		return (DualtreeTransform2dRealInverseSingleTree (wave1rst, wave,
			root, tree, true));
	}

	private double[,] DualtreeTransform2dRealInverseSingleTree
		(IWaveletDTCW wave1rst, IWaveletDTCW wave, Dualtree2dReal root,
		int tree, bool first)
	{
		return (DualtreeTransform2dRealInverseSingleTree (wave1rst, wave,
			root, tree, tree, first));
	}

	private double[,] DualtreeTransform2dRealInverseSingleTree
		(IWaveletDTCW wave1rst, IWaveletDTCW wave, Dualtree2dReal root,
		int t1, int t2, bool first)
	{
		/*Console.WriteLine ("DualtreeTransform2dRealInverseSingleTree (tree {0}|{1}, first {2})",
			t1, t2, first);*/

		double[,] subTrend;