waveletprocessing.cs

Image Fusion Techniues

		if (root.LastLevel) {
			subTrend = (double[,]) root[0];
		} else {
			subTrend = DualtreeTransform2dRealInverseSingleTree (wave1rst,
				wave, root.Next, t1, t2, false);
		}

		double[,] H = (double[,]) root[2];
		double[,] D = (double[,]) root[3];
		double[,] V = (double[,]) root[1];

		double[,] upTrend = DualtreeTransform2dRealInverseSub (subTrend, H, D, V,
			first ? wave1rst : wave, t1, t2);

		return (upTrend);
	}

	private double[,] DualtreeTransform2dRealInverseSub (double[,] trend,
		double[,] H, double[,] D, double[,] V,
		IWaveletDTCW wave, int t1, int t2)
	{
		double[,] outVal = new double[trend.GetLength (0) * 2,
			trend.GetLength (1) * 2];

		/*Console.WriteLine ("DualtreeTransform2dRealInverseSub: to {0}x{1}",
			outVal.GetLength (0), outVal.GetLength (1));*/

		Inverse2dArr (wave[FilterType.Synthesis, t1, PassType.Trend],
			wave[FilterType.Synthesis, t1, PassType.Fluctuation],
			wave[FilterType.Synthesis, t2, PassType.Trend],
			wave[FilterType.Synthesis, t2, PassType.Fluctuation],
			outVal, trend, H, D, V);

		return (outVal);
	}

	//
	// 1D DT Real Discrete Wavelet Transform
	//

	public Dualtree[] DualtreeTransform1d (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, double[] inVal, int scale)
	{
		SimpleArray saIn = new SimpleArray (inVal);

		Dualtree[] trees = new Dualtree[2];

		for (int n = 0 ; n < 2 ; ++n)
			trees[n] = DualtreeTransform1dSingleTree (wave1rst,
				wave, saIn, scale, n);

		return (trees);
	}

	private Dualtree DualtreeTransform1dSingleTree (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, SimpleArray saIn, int scale, int tree)
	{
		// Build tree, starting from root node
		Dualtree root = DualtreeTransform (wave1rst, tree, saIn);
		Dualtree walk = root;

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

			Dualtree subtree = DualtreeTransform (wave, tree,
				new SimpleArray (trend));

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

		return (root);
	}

	private Dualtree DualtreeTransform (IWaveletDTCW wave, int tree,
		IWaveletArray inVal)
	{
		Dualtree dt = new Dualtree ();

		double[] trend = new double[inVal.Length / 2];
		double[] fluctuation = new double[inVal.Length / 2];

		ConvolveArr (wave[FilterType.Analysis, tree, PassType.Trend],
			wave[FilterType.Analysis, tree, PassType.Fluctuation],
			inVal, new SimpleArray (trend), new SimpleArray (fluctuation));

		dt[0] = trend;
		dt[1] = fluctuation;

		return (dt);
	}


	public double[] DualtreeTransform1dInverse (IWaveletDTCW wave1rst,
		IWaveletDTCW wave, Dualtree[] trees)
	{
		double[][] outTrends = new double[2][];

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

		double[] outData = new double[outTrends[0].Length];
		for (int n = 0 ; n < outTrends[0].Length ; ++n)
			outData[n] = (outTrends[0][n] + outTrends[1][n]) / 2.0; //Math.Sqrt (2.0);

		return (outData);
	}

	private double[] DualtreeTransform1dInverseSingleTree
		(IWaveletDTCW wave1rst, IWaveletDTCW wave, Dualtree root,
		int tree, bool first)
	{
		double[] subTrend;
		double[] subFluctuation = (double[]) root[1];

		if (root.LastLevel) {
			subTrend = (double[]) root[0];
		} else {
			subTrend = DualtreeTransform1dInverseSingleTree (wave1rst, wave,
				root.Next, tree, false);
		}

		double[] upTrend = DualtreeTransformInverse (subTrend,
			subFluctuation, first ? wave1rst : wave, tree);

		return (upTrend);
	}

	// Inverses the trend and fluctuation using the wavelet 'wave' to yield
	// the upsampled trend signal as output. 'tree' denotes if its tree 0 or
	// tree 1.
	private double[] DualtreeTransformInverse (double[] trend,
		double[] fluctuation, IWaveletDTCW wave, int tree)
	{
		double[] outVal = new double[trend.Length * 2];

		InverseArr (wave[FilterType.Synthesis, tree, PassType.Trend],
			wave[FilterType.Synthesis, tree, PassType.Fluctuation],
			new SimpleArray (outVal),
			new SimpleArray (trend), new SimpleArray (fluctuation));

		return (outVal);
	}

	public double[,] tileCombine2d (double[,] trend, double[,] H,
		double[,] D, double[,] V, bool normalize)
	{
		int x = trend.GetLength (0);
		int y = trend.GetLength (1);
		double[,] res = new double[2 * x, 2 * y];

		if (normalize) {
			ImageMap tMap = new ImageMap (trend);
			tMap.Normalize ();
			trend = tMap.ValueArray;

			tMap = new ImageMap (H);
			tMap.Normalize ();
			H = tMap.ValueArray;

			tMap = new ImageMap (D);
			tMap.Normalize ();
			D = tMap.ValueArray;

			tMap = new ImageMap (V);
			tMap.Normalize ();
			V = tMap.ValueArray;
		}

		placeCopy (ref res, 0, y, trend);
		placeCopy (ref res, 0, 0, H);
		placeCopy (ref res, x, 0, D);
		placeCopy (ref res, x, y, V);

		return (res);
	}

	private void placeCopy (ref double[,] res, int wx, int wy, double[,] source)
	{
		for (int y = 0 ; y < source.GetLength (0) ; ++y)
			for (int x = 0 ; x < source.GetLength (1) ; ++x)
				res[y + wy, x + wx] = source[y, x];
	}

	public void Inverse1d (IWavelet wave, double[] function,
		double[] trend, double[] fluctuation)
	{
		SimpleArray saOut = new SimpleArray (function);
		SimpleArray saTrend = new SimpleArray (trend);
		SimpleArray saFluctuation = new SimpleArray (fluctuation);

		Inverse (wave, saOut, saTrend, saFluctuation);
	}

	private void InverseSFBWorking (IWavelet wave, IWaveletArray function,
		IWaveletArray trend, IWaveletArray fluctuation)
	{
		/*Console.WriteLine ("InverseSFB (\"{0}\", function[{1}], ..)",
			wave.Name, function.Length);*/

		double[] loFilt = wave.ScalingNumbersSynthesis;
		double[] hiFilt = wave.WaveletNumbersSynthesis;

		// Upsample
		double[] trend0 = new double[2 * trend.Length];
		double[] fluct0 = new double[2 * fluctuation.Length];
		for (int n = 0 ; n < trend.Length ; ++n) {
			trend0[2 * n] = trend[n];
			fluct0[2 * n] = fluctuation[n];
		}


		int N = 2 * trend.Length;
		int L = loFilt.Length;

		int valClen = loFilt.Length + trend0.Length - 1;
		double[] outValHi = new double[valClen];
		double[] outValLo = new double[valClen];

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

				//Console.WriteLine ("vIdx = {0}, k = {1}", vIdx, k);
				if (vIdx < 0 || vIdx >= trend0.Length)
					continue;

				outValLo[k] += loFilt[j] * trend0[vIdx];
				outValHi[k] += hiFilt[j] * fluct0[vIdx];
				//outFluctuation[k % ON] += hiFilt[j] * inVal[vIdx];
			}
		}

		double[] outVal = new double[valClen];
		for (int n = 0 ; n < valClen ; ++n)
			outVal[n] = outValHi[n] + outValLo[n];

		for (int n = 0 ; n < (L - 2) ; ++n)
			outVal[n] += outVal[n + N];

		int shift = L / 2 - 1;
		for (int n = 0 ; n < N ; ++n)
			function[n] = outVal[(n + shift) % N];

		/*
		foreach (double val in outVal)
			Console.WriteLine ("out: {0}", val);
			*/
	}

	private void Inverse (IWavelet wave, IWaveletArray function,
		IWaveletArray trend, IWaveletArray fluctuation)
	{
		/*Console.WriteLine ("InverseSFB (\"{0}\", function[{1}], ..)",
			wave.Name, function.Length);*/

		double[] loFilt = wave.ScalingNumbersSynthesis;
		double[] hiFilt = wave.WaveletNumbersSynthesis;

		InverseArr (loFilt, hiFilt, function, trend, fluctuation);
	}

	private void InverseArr (double[] loFilt, double[] hiFilt,
		IWaveletArray function,
		IWaveletArray trend, IWaveletArray fluctuation)
	{
		int N = 2 * trend.Length;
		int L = loFilt.Length;

		int valClen = loFilt.Length + N - 1;
		double[] outVal = new double[N];

		// Reset
		for (int n = 0 ; n < N ; ++n)
			function[n] = 0.0;

		// Circular shift in final array assignment
		int shift = -(L / 2 - 1);
		if (shift < 0)
			shift += N;

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

				if (vIdx < 0 || vIdx >= N)
					continue;

				if (vIdx % 2 == 1)
					continue;

				double tv = trend[vIdx / 2];
				double fv = fluctuation[vIdx / 2];

				int dk = k;
				if (k >= N && k < (N + L - 2))
					dk -= N;
				if (dk >= N)
					continue;

				function[(dk + shift) % N] += loFilt[j] * tv + hiFilt[j] * fv;
			}
		}
	}

	// Reconstructs the upsampled signal in function from the four subsignals
	// trend, H, D, V.
	public void Inverse2d (IWavelet wave, double[,] function,
		double[,] trend, double[,] H,
		double[,] D, double[,] V)
	{
		Inverse2dArr (wave.ScalingNumbersAnalysis, wave.WaveletNumbersAnalysis,
			wave.ScalingNumbersAnalysis, wave.WaveletNumbersAnalysis,
			function, trend, H, D, V);
	}

	public void Inverse2dArr (double[] loFiltCol, double[] hiFiltCol,
		double[] loFiltRow, double[] hiFiltRow,
		double[,] function,
		double[,] trend, double[,] H,
		double[,] D, double[,] V)
	{
		if ((function.GetLength (0) % 2) != 0 ||
			(function.GetLength (1) % 2) != 0)
		{
			throw (new ArgumentException ("Output array has to be of even dimensions."));
		}

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


		// Create F_{Row} from V and D
		for (int col = 0 ; col < rowFlucts.GetLength (1) ; ++col) {
			InverseArr (loFiltRow, hiFiltRow,
				new DirectionalArray (DirectionalArray.Direction.Column, rowFlucts, col),
				new DirectionalArray (DirectionalArray.Direction.Column, V, col),
				new DirectionalArray (DirectionalArray.Direction.Column, D, col));
		}

		// Create T_{Row} from H and A (Trend)
		for (int col = 0 ; col < rowTrends.GetLength (1) ; ++col) {
			InverseArr (loFiltRow, hiFiltRow,
				new DirectionalArray (DirectionalArray.Direction.Column, rowTrends, col),
				new DirectionalArray (DirectionalArray.Direction.Column, trend, col),
				new DirectionalArray (DirectionalArray.Direction.Column, H, col));
		}

		// Create F from T_{Row} and F_{Row}
		for (int row = 0 ; row < function.GetLength (0) ; ++row) {
			InverseArr (loFiltCol, hiFiltCol,
				new DirectionalArray (DirectionalArray.Direction.Row, function, row),
				new DirectionalArray (DirectionalArray.Direction.Row, rowTrends, row),
				new DirectionalArray (DirectionalArray.Direction.Row, rowFlucts, row));
		}
	}


	// Wrapper classes so we can share all of the wavelet code between 1d and
	// 2d wavelet operations by mapping 2d arrays to 1d array slices.
	private interface
	IWaveletArray
	{
		int Length {
			get;
		}

		double this[int n] {
			get;
			set;
		}
	}

	private class
	SimpleArray : IWaveletArray
	{
		private SimpleArray ()
		{
		}

		double[] arr;

		public SimpleArray (double[] arr)
		{
			this.arr = arr;
		}

		public int Length {
			get {
				return (arr.Length);
			}
		}

		public double this[int n] {
			get {
				return (arr[n]);
			}
			set {
				arr[n] = value;
			}
		}
	}

	private class
	DirectionalArray : IWaveletArray
	{
		private DirectionalArray ()
		{
		}

		public enum Direction {
			Row,
			Column,
		};

		Direction dir;
		double[,] arr;
		int resolveDim;

		public DirectionalArray (Direction dir, double[,] arr, int resolveDim)
		{
			this.dir = dir;
			this.arr = arr;
			this.resolveDim = resolveDim;
		}

		public int Length {
			get {
				if (dir == Direction.Row)
					return (arr.GetLength (1));
				else
					return (arr.GetLength (0));
			}
		}

		public double this[int n] {
			get {
				/*Console.WriteLine ("get: {0} of {1}, dimension {2}",
					n, (dir == Direction.Row) ? arr.GetLength (1) : arr.GetLength (0),
					dir);*/

				if (dir == Direction.Row)
					return (arr[resolveDim, n]);
				else
					return (arr[n, resolveDim]);
			}

			set {
				if (dir == Direction.Row)
					arr[resolveDim, n] = value;
				else
					arr[n, resolveDim] = value;
			}
		}
	}
}