// 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'.///* TestWavelet.cs * * Testing program for all kinds of supported wavelet transforms. */using Gtk;using GtkSharp;using System;public classTestWavelet{	public static void Main (string[] args)	{		if (args.Length < 1) {			Console.WriteLine ("usage: TestWavelet.exe image");			return;		}		Application.Init ();		Window win = new Window ("Test Wavelet");		win.DefaultSize = new Gdk.Size (400, 300);		win.DeleteEvent += new DeleteEventHandler (Window_Delete);		TestWavelet test = new TestWavelet ();		test.Test1D ();		Gtk.VBox vb = new Gtk.VBox ();		test.Test1DDualtree ();		//test.Test1DComplex ();		Gtk.Widget resNoise = test.AddNoise (args[0]);		vb.Add (resNoise);		Gtk.Widget resComplexPrint = test.PrintDualtree2DComplexWavelets ();		Gtk.Widget resComplex = test.TestDualtree2DComplex (args[0]);		vb.Add (resComplexPrint);		vb.Add (resComplex);		Gtk.Widget resRealPrint = test.PrintDualtree2DRealWavelets ();		Gtk.Widget resReal = test.TestDualtree2DReal (args[0]);		//Gtk.Widget res = test.Test2D (args[0]);		vb.Add (resRealPrint);		vb.Add (resReal);		win.Add (vb);		win.ShowAll ();		Application.Run ();	}	public Gtk.Widget AddNoise (string filename)	{		DisplayImage pic = new DisplayImage (filename);		ImageMap map = pic.ConvertToImageMap			(new DisplayImage.CanonicalPixelConverter ());		VBox vbLine = new VBox ();		HBox hb = new HBox ();		hb.Add (new Image (pic.Pbuf));		ImageMap noisy01 = (ImageMap) map.Clone ();		Denoise.InsertRandomNoise (noisy01, 0.1);		hb.Add (new Image (new DisplayImage (noisy01, 1.0).Pbuf));		ImageMap noisy02 = (ImageMap) map.Clone ();		Denoise.InsertRandomNoise (noisy02, 0.2);		hb.Add (new Image (new DisplayImage (noisy02, 1.0).Pbuf));		ImageMap noisy05 = (ImageMap) map.Clone ();		Denoise.InsertRandomNoise (noisy05, 0.5);		hb.Add (new Image (new DisplayImage (noisy05, 1.0).Pbuf));		// Denoise		double[] denoiseVal = { 0.01, 0.02, 0.05, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.4, };		int p = 0;		foreach (double thresh in denoiseVal) {			ImageMap denoised = Denoise.DenoiseSoft (noisy02, thresh);			Gtk.VBox vb = new Gtk.VBox ();			vb.Add (new Image (new DisplayImage (denoised, 1.0).Pbuf));			vb.Add (new Gtk.Label (String.Format ("denoised, thresh {0}", thresh)));			if (p % 4 == 0) {				vbLine.Add (hb);				hb = new HBox ();			}			p += 1;			hb.Add (vb);		}		vbLine.Add (hb);		return (vbLine);	}	public Gtk.Widget Test2D (string filename)	{		DisplayImage pic = new DisplayImage (filename);		ImageMap map = pic.ConvertToImageMap			(new DisplayImage.CanonicalPixelConverter ());		// Do the wavelet transform		double[,] picA = new double[map.YDim / 2, map.XDim / 2];		double[,] picH = new double[map.YDim / 2, map.XDim / 2];		double[,] picD = new double[map.YDim / 2, map.XDim / 2];		double[,] picV = new double[map.YDim / 2, map.XDim / 2];		IWavelet wave = new Farras10 ();		double[,] picArr = map.ValueArray;		WaveletProcessing wp = new WaveletProcessing ();		wp.Convolve2d (wave, picArr, picA, picH, picD, picV);		HBox hb = new HBox ();		hb.Add (new Image (pic.Pbuf));		hb.Add (new Image (new DisplayImage (map, 1.0).Pbuf));		/*		double[][,] vArr = { picA, picH, picD, picV };		foreach (double[,] arr in vArr) {			ImageMap subMap = new ImageMap (arr);			subMap.Normalize ();			DisplayImage subPic = new DisplayImage (subMap, 1.0);			hb.Add (new Image (subPic.Pbuf));		}		*/		double[,] comb = wp.tileCombine2d (picA, picH, picD, picV, false);		ImageMap sMap = new ImageMap (comb);		sMap.Normalize ();		DisplayImage sPic = new DisplayImage (sMap, 1.0);		hb.Add (new Image (sPic.Pbuf));		double[,] restored = new double[picArr.GetLength (0), picArr.GetLength (1)];		wp.Inverse2d (wave, restored, picA, picH, picD, picV);		hb.Add (new Image (new DisplayImage (new ImageMap (restored), 1.0).Pbuf));		return (hb);	}	// Display the six complex fluctuation wavelets of the 2d complex dualtree DWT	public Gtk.Widget PrintDualtree2DComplexWavelets ()	{		int J = 4;		int L = 3 * (1 << (J + 1));		int N = L / (1 << J);		WaveletProcessing wp = new WaveletProcessing ();		IWaveletDTCW wave1rst = new Farras10CFirst ();		IWaveletDTCW wave = new KingsburyQC10 ();		double[,] input = new double[3 * L, 6 * L];		Dualtree2dComplex[] dt = wp.DualtreeTransform2dComplex (wave1rst, wave,			input, J);		Dualtree2dComplex tree1 = dt[0];		Dualtree2dComplex tree2 = dt[1];		while (tree1.LastLevel == false) {			tree1 = tree1.Next;			tree2 = tree2.Next;		}		double[][,] flucts = {			tree2.GetDirectStore (0, 1), tree1.GetDirectStore (0, 2),			tree2.GetDirectStore (0, 0), tree1.GetDirectStore (0, 0),			tree2.GetDirectStore (0, 2), tree1.GetDirectStore (0, 1),			tree2.GetDirectStore (1, 1), tree1.GetDirectStore (1, 2),			tree2.GetDirectStore (1, 0), tree1.GetDirectStore (1, 0),			tree2.GetDirectStore (1, 2), tree1.GetDirectStore (1, 1),		};		int tNH = flucts.Length / 2;		for (int t = 0 ; t < flucts.Length ; ++t)			flucts[t][N/2 + (t/tNH) * N, N/2 + (t % tNH) * N] = 1.0;		double[,] restored = wp.DualtreeTransform2dComplexInverse (wave1rst, wave, dt);		// Third row: magnitude product of real and imaginary part		for (int y = 0 ; y < L ; ++y) {			for (int x = 0 ; x < (6 * L) ; ++x) {				restored[2 * L + y, x] = Math.Sqrt (					Math.Pow (restored[y, x], 2.0) + Math.Pow (restored[L + y, x], 2.0));			}		}		ImageMap map = new ImageMap (restored);		map.Normalize ();		HBox hb = new HBox ();		hb.Add (new Image (new DisplayImage (map, 1.0).Pbuf));		return (hb);	}	public Gtk.Widget TestDualtree2DComplex (string filename)	{		Console.WriteLine ("=== TestDualtree2DComplex BEGIN ===\n");		DisplayImage pic = new DisplayImage (filename);		ImageMap map = pic.ConvertToImageMap			(new DisplayImage.CanonicalPixelConverter ());		double[,] picArr = map.ValueArray;		WaveletProcessing wp = new WaveletProcessing ();		IWaveletDTCW wave1rst = new Farras10CFirst ();		IWaveletDTCW wave = new KingsburyQC10 ();		Dualtree2dComplex[] dt = wp.DualtreeTransform2dComplex (wave1rst, wave,			picArr, 2);		HBox hb = new HBox ();		hb.Add (new Image (pic.Pbuf));		hb.Add (new Image (new DisplayImage (map, 1.0).Pbuf));		for (int cn = 0 ; cn < 4 ; ++cn) {			double[,] comb = wp.tileCombine2d				(dt[cn / 2].GetTrend (cn % 2),				dt[cn / 2].GetDirectStore (cn % 2, 0),				dt[cn / 2].GetDirectStore (cn % 2, 1),				dt[cn / 2].GetDirectStore (cn % 2, 2), true);			ImageMap sMap = new ImageMap (comb);			DisplayImage sPic = new DisplayImage (sMap, 1.0);			hb.Add (new Image (sPic.Pbuf));		}		for (int cn = 0 ; cn < 4 ; ++cn) {			double[][,] subbands = {				dt[cn / 2].GetDirectStore (cn % 2, 0),				dt[cn / 2].GetDirectStore (cn % 2, 1),				dt[cn / 2].GetDirectStore (cn % 2, 2),			};			for (int band = 0 ; band < subbands.Length ; ++band) {				double sum = 0.0;				int count = 0;				double[,] plane = subbands[band];				for (int y = 0 ; y < plane.GetLength (0) ; ++y) {					for (int x = 0 ; x < plane.GetLength (1) ; ++x) {						count += 1;						sum += Math.Abs (plane[y, x]);						/*					if (Math.Abs (plane[y, x]) < 0.1)						plane[y, x] = 0.0;						*/					/*					double s = Math.Max (Math.Abs (plane[y, x] - 0.02), 0.0);					s = (s / (s + 0.02)) * plane[y, x];					plane[y, x] = s;					*/					}				}				Console.WriteLine ("band {0} average: {1}", band,					sum / ((double) count));			}		}		Console.WriteLine ("= Reconstructing...");		double[,] restored = wp.DualtreeTransform2dComplexInverse (wave1rst, wave, dt);		hb.Add (new Image (new DisplayImage (new ImageMap (restored), 1.0).Pbuf));		Console.WriteLine ("=== TestDualtree2DComplex END ===\n");		return (hb);	}	// Display the three fluctuation wavelets of the 2d real dualtree DWT	public Gtk.Widget PrintDualtree2DRealWavelets ()	{		int J = 4;		int L = 3 * (1 << (J + 1));		int N = L / (1 << J);		WaveletProcessing wp = new WaveletProcessing ();		IWaveletDTCW wave1rst = new Farras10CFirst ();		IWaveletDTCW wave = new KingsburyQC10 ();		double[,] input = new double[3 * L, 3 * L];		Dualtree2dReal[] dt = wp.DualtreeTransform2dReal (wave1rst, wave,			input, J);		Dualtree2dReal tree1 = dt[0];		Dualtree2dReal tree2 = dt[1];		while (tree1.LastLevel == false) {			tree1 = tree1.Next;			tree2 = tree2.Next;		}		double[,][,] flucts = {			{ (double[,]) tree1[1], (double[,]) tree1[2], (double[,]) tree1[3], },			{ (double[,]) tree2[1], (double[,]) tree2[2], (double[,]) tree2[3], },		};		for (int t = 0 ; t < 2 ; ++t)			for (int m = 0 ; m < 3 ; ++m)				flucts[t, m][N/2 + t * N, N/2 + m * N] = 1.0;		double[,] restored = wp.DualtreeTransform2dRealInverse (wave1rst, wave, dt);		for (int y = 0 ; y < L ; ++y) {			for (int x = 0 ; x < (3 * L) ; ++x) {				restored[2 * L + y, x] = Math.Sqrt (					Math.Pow (restored[y, x], 2.0) + Math.Pow (restored[L + y, x], 2.0));			}		}		ImageMap map = new ImageMap (restored);		map.Normalize ();		HBox hb = new HBox ();		hb.Add (new Image (new DisplayImage (map, 1.0).Pbuf));		return (hb);	}	public Gtk.Widget TestDualtree2DReal (string filename)	{		DisplayImage pic = new DisplayImage (filename);		ImageMap map = pic.ConvertToImageMap			(new DisplayImage.CanonicalPixelConverter ());		double[,] picArr = map.ValueArray;		WaveletProcessing wp = new WaveletProcessing ();		IWaveletDTCW wave1rst = new Farras10CFirst ();		IWaveletDTCW wave = new KingsburyQC10 ();		Dualtree2dReal[] dt = wp.DualtreeTransform2dReal (wave1rst, wave,			picArr, 2);		HBox hb = new HBox ();		hb.Add (new Image (pic.Pbuf));		hb.Add (new Image (new DisplayImage (map, 1.0).Pbuf));		double[,] comb1 = wp.tileCombine2d			((double[,]) dt[0][0], (double[,]) dt[0][1],			(double[,]) dt[0][2], (double[,]) dt[0][3], true);		double[,] comb2 = wp.tileCombine2d			((double[,]) dt[1][0], (double[,]) dt[1][1],			(double[,]) dt[1][2], (double[,]) dt[1][3], true);		ImageMap sMap = new ImageMap (comb1);		//sMap.Normalize ();		DisplayImage sPic = new DisplayImage (sMap, 1.0);		hb.Add (new Image (sPic.Pbuf));		sMap = new ImageMap (comb2);		//sMap.Normalize ();		sPic = new DisplayImage (sMap, 1.0);		hb.Add (new Image (sPic.Pbuf));		double[][,] flucts = { (double[,]) dt[0].Next[1],			(double[,]) dt[0].Next[2], (double[,]) dt[0].Next[3],			(double[,]) dt[1].Next[1], (double[,]) dt[1].Next[2], (double[,]) dt[1].Next[3], };		foreach (double[,] plane in flucts) {			double sum = 0.0;			int count = 0;			for (int y = 0 ; y < plane.GetLength (0) ; ++y) {				for (int x = 0 ; x < plane.GetLength (1) ; ++x) {					count += 1;					sum += Math.Abs (plane[y, x]);					/*if (Math.Abs (plane[y, x]) < 0.02)						plane[y, x] = 0.0;*/					/*					double s = Math.Max (Math.Abs (plane[y, x] - 0.02), 0.0);					s = (s / (s + 0.02)) * plane[y, x];					plane[y, x] = s;					*/				}			}			Console.WriteLine ("plane average: {0}", sum / ((double) count));		}		Console.WriteLine ("= Reconstructing...");		double[,] restored = wp.DualtreeTransform2dRealInverse (wave1rst, wave, dt);		hb.Add (new Image (new DisplayImage (new ImageMap (restored), 1.0).Pbuf));		return (hb);	}	public void Test1DDualtree ()	{		Console.WriteLine ("\n=== Test1DDualtree BEGIN ===");		double[] numBase = {			2.0, 3.0, 6.0, 1.0,			4.0, 6.0, 8.0, 10.0,			-4.0, -2.0, -1.0, 1.0,			2.0, 3.3, 3.4, 4.0,			0.0, 7.0, 5.0, 11.0,			14.0, -6.0, -7.0, 1.0,			-2.0, -2.0, 1.0, 3.0,			4.0, 3.0, 2.0, 1.0, };		double[] num = new double[numBase.Length * 8];		for (int n = 0 ; n < num.Length ; ++n)			num[n] = numBase[n % numBase.Length] * ((n / numBase.Length) + 1.0);		IWaveletDTCW wave1rst = new Farras10CFirst ();		IWaveletDTCW wave = new KingsburyQC10 ();		WaveletProcessing wp1d = new WaveletProcessing ();		Console.WriteLine ("num[] =");		foreach (double val in num)			Console.WriteLine ("\t{0}", val);		Dualtree[] dt = wp1d.DualtreeTransform1d (wave1rst, wave, num, 5);		for (int t = 0 ; t < 2 ; ++t) {			Console.WriteLine ("Tree {0}", t + 1);			Dualtree tree = dt[t];			int level = 0;			while (tree != null) {				Console.WriteLine ("    Level {0}", level);				for (int signal = 0 ; signal < 2 ; ++signal) {					Console.WriteLine ("        {0}",						signal == 0 ? "Trend" : "Fluctuation");					foreach (double val in ((double[]) tree[signal]))						Console.WriteLine ("            {0}", val);				}				tree = tree.Next;				level += 1;			}		}		Console.WriteLine ("\n  === Reconstruction ===\n");		double[] reconstruction = wp1d.DualtreeTransform1dInverse			(wave1rst, wave, dt);		Console.WriteLine ("reconstruction[] =");		foreach (double val in reconstruction)			Console.WriteLine ("\t{0}", val);		Console.WriteLine ("\n=== Test1DDualtree END ===");	}	public void Test1D ()	{		//double[] num = { 10.0, 8.0, 7.0, 3.0, -3.0, 5.0, 8.0, 9.0, };		double[] num = { 2.0, 3.0, 6.0, 1.0, 4.0, 6.0, 8.0, 10.0, -4.0, -2.0, -1.0, 1.0 };		IWavelet wave1d = new Farras10 ();		WaveletProcessing wp1d = new WaveletProcessing ();		double[] trend = new double[6];		double[] fluctuation = new double[6];		wp1d.Convolve1d (wave1d, ref num, ref trend, ref fluctuation);		Console.WriteLine ("Trend");		foreach (double val in trend)			Console.Write ("{0}, ", val);		Console.WriteLine ("\nFluctuation");		foreach (double val in fluctuation)			Console.Write ("{0}, ", val);		Console.WriteLine ("");		num[0] = -1.3;		wp1d.Inverse1d (wave1d, num, trend, fluctuation);		Console.WriteLine ("\nInversed function again:");		foreach (double val in num)			Console.Write ("{0}, ", val);		Console.WriteLine ("");	}	static void Window_Delete (object obj, DeleteEventArgs args)	{		Application.Quit ();		args.RetVal = true;	}}