denoise.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'.
//

using System;


public class
Denoise
{
	// 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);
	}
}