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

/* DOFProcessing.cs - Extended Depth of Focus processing algorithms
 *
 * algorithms based on the short descriptions in the research paper
 * "Extended depth-of-focus for multi-channel microscopy images: a complex
 * wavelet approach" by Brigitte Forster, Dimitri Van De Ville, Jesse Berent,
 * Daniel Sage, Michael Unser.
 *
 * Available online at http://bigwww.epfl.ch/preprints/forster0401p.html
 *
 * NOTICE: so far, only real processing has been implemented, because of
 * outstanding bugs in the complex wavelet transform algorithm. The algorithm
 * in this file would require no internal changes (except C# type changes),
 * though. Do not compare the depth-of-focus fusion performance of this
 * program with the ones presented in the research paper, as this algorithm is
 * not as powerful yet.
 *
 * (C) Copyright 2004 -- Sebastian Nowozin ([e-mail])
 */

using System;
using System.Collections;


public class
DOFProcessing
{
	WaveletProcessing wp = new WaveletProcessing ();
	IWaveletDTCW wave1rst = new Farras10CFirst ();
	IWaveletDTCW wave = new KingsburyQC10 ();

	int currentCount = 0;
	ArrayList originals = new ArrayList ();
	public ArrayList Originals {
		get {
			return (originals);
		}
	}

	int scales = -1;
	int dimension = -1;

	Dualtree2dComplex[] dt = null;

	// All images have to have the same dimension.
	int xDim = -1;
	int yDim = -1;

	public DOFProcessing ()
	{
	}

	public DOFProcessing (int scales)
	{
		this.scales = scales;
	}

	public ImageMap ProduceImage ()
	{
		double[,] restored = wp.DualtreeTransform2dComplexInverse (wave1rst, wave, dt);

		ImageMap fused = new ImageMap (restored);
		fused = fused.ShrinkCanvas (xDim, yDim);

		return (fused);
	}

	// Fuse one image into the overall image
	// Return the number of wavelet coefficients replaced
	public int FuseOne (string filename)
	{
		DisplayImage pic = new DisplayImage (filename);
		ImageMap map = pic.ConvertToImageMap
			(new DisplayImage.CanonicalPixelConverter ());

		return (FuseOne (map));
	}

	public int FuseOne (ImageMap map)
	{
		currentCount += 1;
		originals.Add (map);

		if ((xDim != -1 && xDim != map.XDim) || (yDim != -1 && yDim != map.YDim)) {
			throw (new ArgumentException (String.Format
				("Image files have to be of the same dimension, required: {0}x{1}, found: {2}x{3}",
				xDim, yDim, map.XDim, map.YDim)));
		}

		if (dt == null) {
			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;

			scales = 0;
			while ((1 << scales) < larger)
				scales += 1;

			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)));
			}
			xDim = map.XDim;
			yDim = map.YDim;

			Console.WriteLine ("Images are ({0}x{1}), upscaled to ({2}x{3}), using {4} scales.",
				xDim, yDim, dimension, dimension, scales);
		}

		ImageMap upscaled = map.EnlargeCanvas (dimension, dimension);

		// Forward transform
		Dualtree2dComplex[] dimage = wp.DualtreeTransform2dComplex (wave1rst, wave,
			upscaled.ValueArray, scales);

		if (dt == null) {
			dt = dimage;

			return (-1);
		}

		// "The largest absolute value of the coefficients in the subbands
		// will correspond to sharper brightness changes and therefore to the
		// most salient features. A good integration rule consists in
		// selecting the slice with the largest absolute value of the wavelet
		// coefficients at each point." (Forster, see top of this file).
		//
		// It is however, not quite clear if each subband should be treated
		// independently, or one slice dominates all subbands. Later in the
		// paper, a "Subband consistency" check is introduced, from which I
		// deduce that most likely each subband shall be treated
		// independently.
		int replaced = 0;

		Console.WriteLine ("Merging image coefficients...");
		for (int treeNumber = 0 ; treeNumber < 2 ; ++treeNumber) {
			Dualtree2dComplex treeNew = dimage[treeNumber];
			Dualtree2dComplex treeBase = dt[treeNumber];

			int tLevel = 0;

			while (treeNew != null) {
				int subReplaced = SalienceComplexLevel (treeBase, treeNew);

				Console.WriteLine ("    tree {0} level {1}, replaced {2} coefficients",
					treeNumber, tLevel, subReplaced);
				replaced += subReplaced;

				treeNew = treeNew.Next;
				treeBase = treeBase.Next;
				tLevel += 1;
			}
		}

		Console.WriteLine ("=> TOTAL replaced {0} coefficients", replaced);

		return (replaced);
	}

	private int SalienceComplexLevel (Dualtree2dComplex treeBase,
		Dualtree2dComplex treeNew)
	{
		int replaced = 0;

		for (int band = 0 ; band < 3 ; ++band) {
			Dualtree2dComplex.ComplexArray baseBand = treeBase[band];
			Dualtree2dComplex.ComplexArray newBand = treeNew[band];

			for (int y = 0 ; y < baseBand.GetLength (0) ; ++y) {
				for (int x = 0 ; x < baseBand.GetLength (1) ; ++x) {
					double absNewBand = Math.Sqrt (
						Math.Pow (newBand[y, x, 0], 2.0) +
						Math.Pow (newBand[y, x, 1], 2.0));
					double absBaseBand = Math.Sqrt (
						Math.Pow (baseBand[y, x, 0], 2.0) +
						Math.Pow (baseBand[y, x, 1], 2.0));

					if (absNewBand > absBaseBand) {
						baseBand[y, x, 0] = newBand[y, x, 0];
						baseBand[y, x, 1] = newBand[y, x, 1];
						replaced += 1;
					}
				}
			}
		}

		return (replaced);
	}

	private int ChooseSalient (ImageMap p, ImageMap c)
	{
		int replaced = 0;

		for (int y = 0 ; y < p.YDim ; ++y) {
			for (int x = 0 ; x < p.XDim ; ++x) {
				if (Math.Abs (c[x, y]) > Math.Abs (p[x, y])) {
					p[x, y] = c[x, y];
					replaced += 1;
				}
			}
		}

		return (replaced);
	}
}


public class
MultiClipping
{
	public static int[,] ClippingMap (ArrayList originals, ImageMap fused)
	{
		ImageMap input = (ImageMap) originals[0];
		int maxX = input.XDim;
		int maxY = input.YDim;

		int[,] top = new int[maxY, maxX];

		for (int y = 0 ; y < maxY ; ++y) {
			for (int x = 0 ; x < maxX ; ++x) {
				double minDiff = Double.MaxValue;

				// Find the minimum delta original
				for (int n = 0 ; n < originals.Count ; ++n) {
					ImageMap map = (ImageMap) originals[n];

					double diff = Math.Abs (map[x, y] - fused[x, y]);
					if (diff < minDiff) {
						minDiff = diff;
						top[y, x] = n;
					}
				}
			}
		}

		return (top);
	}
}