// 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'.///* 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 (nowozin@cs.tu-berlin.de) */using System;using System.Collections;public classDOFProcessing{	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 classMultiClipping{	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);	}}