// 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'.//using System;public classDenoise{	// 0.0 <= noiseFactor <= 1.0	public static void InsertRandomNoise (ImageMap map, double noiseFactor)	{		Random rnd = new Random ();		for (int y = 0 ; y < map.YDim ; ++y) {			for (int x = 0 ; x < map.XDim ; ++x) {				map[x, y] += rnd.NextDouble () * noiseFactor;				if (map[x, y] > 1.0)					map[x, y] = 1.0;			}		}	}	public static ImageMap DenoiseSoft (ImageMap map, double threshhold)	{		IWaveletDTCW wave1rst = new Farras10CFirst ();		IWaveletDTCW wave = new KingsburyQC10 ();		WaveletProcessing wp = new WaveletProcessing ();		int smaller = map.XDim < map.YDim ? map.XDim : map.YDim;		int larger = map.XDim > map.YDim ? map.XDim : map.YDim;		int wavelen = wave1rst.SupportSize > wave.SupportSize ?			wave1rst.SupportSize : wave.SupportSize;		// Get the largest scale we could do.		int wscales = 0;		while ((1 << wscales) < wavelen)			wscales += 1;		int scales = 0;		while ((1 << scales) < larger)			scales += 1;		int dimension = 1 << scales;		scales -= wscales;		if (scales <= 0) {			throw (new ArgumentException (String.Format				("Input too small: size {0} is below wavelet length {1}.",				smaller, wavelen)));		}		Console.WriteLine ("Denoise: Using {0} scales.", scales);		int xDim = map.XDim;		int yDim = map.YDim;		ImageMap upscaled = map.EnlargeCanvas (dimension, dimension);		// Forward transform		Dualtree2dComplex[] dimage = wp.DualtreeTransform2dComplex (wave1rst, wave,			upscaled.ValueArray, scales);		// Soft denoising.		int soft = 0;		for (int treeNumber = 0 ; treeNumber < 2 ; ++treeNumber) {			Dualtree2dComplex tree = dimage[treeNumber];			int tLevel = 0;			while (tree != null) {				int subSoft = SoftThreshhold (tree, threshhold);				Console.WriteLine ("tree {0} level {1}, softed {2} coefficients",					treeNumber, tLevel, subSoft);				soft += subSoft;				tree = tree.Next;			}		}		double[,] restored = wp.DualtreeTransform2dComplexInverse (wave1rst,			wave, dimage);		ImageMap fused = new ImageMap (restored);		return (fused.ShrinkCanvas (xDim, yDim));	}	private static int SoftThreshhold (Dualtree2dComplex tree,		double threshhold)	{		int softed = 0;		for (int band = 0 ; band < 3 ; ++band) {			Dualtree2dComplex.ComplexArray tband = tree[band];			for (int y = 0 ; y < tband.GetLength (0) ; ++y) {				for (int x = 0 ; x < tband.GetLength (1) ; ++x) {					double realV = tband[y, x, 0];					double complexV = tband[y, x, 1];					double magnitude = Math.Sqrt (Math.Pow (realV, 2.0) +						Math.Pow (complexV, 2.0));					// Threshholding					if (magnitude < threshhold) {						softed += 1;						magnitude = 0.0;					}					realV = (magnitude / (magnitude + threshhold)) * realV;					complexV = (magnitude / (magnitude + threshhold)) * complexV;					tband[y, x, 0] = realV;					tband[y, x, 1] = complexV;				}			}		}		return (softed);	}}