// Waveblend - complex dualtree based image fusion// (C) Copyright 2004 -- Sebastian Nowozin <nowozin@cs.tu-berlin.de>//// 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 classWaveletProcessing : 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;