public final class Tool {

    /** Forbit the construction of this class */
    private Tool() {}
    /** The tool implementation in use*/
    static private Tool.Logger currentLogger;
    static private boolean errorShown = false;

    /** Simple logger */
    public interface Logger {
	public void writeTrace(String message); 
	public void writeError(String message, boolean fatal);
	public void writeShortMessage(String message);
    }

    /** Inits a standard tool */
    static {
	// basic logger goes to System.out
	currentLogger = new Tool.Logger() {
	    public void writeTrace(String message) {
		System.out.println( "[fairSIM] "+message);
		System.out.flush();
	    }
	    
	    public void writeError(String message, boolean fatal) {
		String prefix = (fatal)?("[fsFATAL]"):("[fsERROR]");
		System.err.println( prefix+" "+message);
		System.err.flush();
	    }
    
	    public void writeShortMessage(String message) {
		System.out.println( "-fairSIM- "+message);
		System.out.flush();
	    }
	};
	// want to catch exceptions
	/*
	Thread.setDefaultUncaughtExceptionHandler( new Thread.UncaughtExceptionHandler () {
	    public void uncaughtException(Thread t, Throwable e) {
		Tool.trace("Problem, caugt exception: "+e);
		e.printStackTrace();
	    }
	});
	*/
    }

    /** Write a trace message */
    static public final void trace(String message) {
	if (currentLogger!=null)
	    currentLogger.writeTrace( message );
    }
    
    /** Output a short status / info message */
    static public final void tell(String message) {
	if (currentLogger!=null)
	    currentLogger.writeShortMessage( message);
    }
   
    /** Write an error message */
    static public final void error(final String message, boolean fatal ) {
        new Thread(new Runnable() {
            public void run() {
                if (!errorShown) {
                    errorShown = true;
                    JOptionPane.showMessageDialog(null, message, "fairSIM Error", JOptionPane.ERROR_MESSAGE);
                    try {
                        Thread.sleep(3000);
                    } catch (InterruptedException ex) {
                        Tool.trace("Tool: errorShown sleep interrupted, why?");
                    }
                    errorShown = false;
                }
            }
        }).start();
	if (currentLogger!=null)
	    currentLogger.writeError( message, fatal );
    }

    /** Creates an absolute file from path provided a string.
     *  This especially takes care of "~" to set the users home */
    static public File getFile(String path) {
	if (path.startsWith("~"+File.separator)){
	    path=System.getProperty("user.home")+path.substring(1);
	}
	return new File(path).getAbsoluteFile();
    }


    /** Decode a BCD timestamp (as used by PCO) 
     *	@param stamp input, typically image acquired by camera (first 16 entries used)
     *	@return The timestamp, in microseconds since epoch
     * */
    static public long decodeBcdTimestamp( short [] stamp ) {
	long stampNr= bcdDecode(stamp, 0, 4);

	int year = bcdDecode(stamp, 4, 6);
	int month= bcdDecode(stamp, 6, 7);
	int day  = bcdDecode(stamp, 7, 8);

	int h	  = bcdDecode(stamp,  8, 9);
	int min  = bcdDecode(stamp,  9,10);
	int sec  = bcdDecode(stamp, 10,11);
	int  us  = bcdDecode(stamp, 11,14);

	Calendar cld = Calendar.getInstance();
	cld.set( year, month-1, day, h, min, sec );

	long ret = cld.getTimeInMillis()/1000; 
	// TODO: for some reason, 'getTimeInMillis' is not a multiple of 1000
	
	ret = (ret*1000000) +us;

	return ret;
    }

    /** Decode double-packed (2 digit per byte) BCD-encoded values.
     *  Beware of overflows, even with long
     *  @param arr The input array to decode from 
     *	@param start start of range to decode
     *	@param end end of range to decode
     *	@return the decoded number 
     *  */
    public static int bcdDecode( short [] arr, int start, int end) {
	int ret=0;
	int count=0;
	
	for (int j=end-1; j>=start; j--) {
	    int val1 =  arr[j] & 0x000F ;
	    int val2 = (arr[j] & 0x00F0)>>4 ;
	    int mult = (int)Math.pow(10, count*2);
	    //IJ.log(""+mult);
	    ret += (val1+val2*10)*mult;
	    count++;
	}

	return ret;
    }




    /** Format a 'milliseconds since 1 Jan 1970' timestamp in ISO */ 
    static public String readableTimeStampMillis( long ms , boolean spaces ) {
	java.text.DateFormat df;
	if (!spaces) {
	    df = new java.text.SimpleDateFormat("yyyyMMdd'T'HHmmss'Z'");
	    df.setTimeZone(java.util.TimeZone.getTimeZone("UTC"));
	} else {
	    df = new java.text.SimpleDateFormat("yyyy-MM-dd' T 'HH:mm:ss '('Z')'");
	}

	String nowAsISO = df.format(new java.util.Date(ms));
	return nowAsISO;
    }

    /** Format a 'seconds since 1 Jan 1970' timestamp in ISO */ 
    static public String readableTimeStampSeconds( double seconds , boolean spaces) {
	long val = (long)(seconds*1000);
	return readableTimeStampMillis(val, spaces);
    }

    /** Implement and pass Tool.Logger to redirect log output,
     *  or set null to disable output completely */
    public static void setLogger( Tool.Logger t ) {
	currentLogger = t;
    }

    /** Shuts down all multi-threading pools */
    public static void shutdown() {
	SimpleMT.shutdown();
    }	

    /** Return a Tool.Timer, which is automatically started. */
    static public Timer getTimer() { return new Timer(); };
    
    /** A simple timer. TODO: The meaning of stop, pause, ... could
     *  be much clearer */
    public static class Timer {
	long start, stop, runtime, outtime;
	Timer() { 
	    start();
	}
	/** start the timer */
	public void start() { 
	    //start = System.currentTimeMillis(); 
	    start =  System.nanoTime(); 
	};
	/** stop the timer (next start resets it) */
	public void stop() { 
	    //stop = System.currentTimeMillis(); 
	    stop = System.nanoTime(); 
	    runtime += stop-start;
	    outtime=runtime;
	    runtime=0;
	    }
	/** pause the timer (next start continues) */
	public void hold(){
	    //stop = System.currentTimeMillis();
	    stop = System.nanoTime(); 
	    runtime += stop-start;
	    outtime  = runtime;
	    start =stop;
	}
	/** get the milliseconds on this timer */
	public double msElapsed() {
	    return outtime/1000000.;
	}

	/** output the amount of milliseconds counted */
	@Override public String toString(){ 
	    return String.format("%10.3f ms",(outtime/1000000.));
	}
    }
    
    public static String[] decodeArray(String encodedArray) {
        String[] split = encodedArray.split(";");
        String[] data = new String[split.length - 1];
        for (int i = 0; i < data.length; i++) {
            data[i] = split[i + 1];
        }
        return data;
    }
    
    public static int[] decodeIntArray(String encodedArray) {
        String[] stringArray = decodeArray(encodedArray);
        int len = stringArray.length;
        int[] intArray = new int[len];
        for (int i = 0; i < len; i++) {
            intArray[i] = Integer.parseInt(stringArray[i]);
        }
        return intArray;
    }

    private static String encodeArray(String prefix, String[] array) {
        int len = array.length;
        String output = prefix;
        if (len == 0) {
            return output;
        } else {
            for (int i = 0; i < len; i++) {
                output += ";" + array[i];
            }
            return output;
        }
    }
    
    public static String encodeArray(String prefix, int[] array) {
        int len = array.length;
        String[] sArray = new String[len];
        for (int i = 0; i < len; i++) {
            sArray[i] = Integer.toString(array[i]);
        }
        return encodeArray(prefix, sArray);
    }
    
    public static <T> String encodeArray(String prefix, T[] array) {
        if (array instanceof String[]) {
            return encodeArray(prefix, (String[])array);
        } else {
            int len = array.length;
            String[] sArray = new String[len];
            for (int i = 0; i < len; i++) {
                sArray[i] = array[i].toString();
            }
            return encodeArray(prefix, sArray);
        }
    }

    /** A generic callback interface */
    public static interface Callback<T> {
	public void callback(T a);
    }

    /** A generic tuple */
    public static class Tuple<F,S> {
	public final F first;
	public final S second;

	public Tuple (F first, S second) {
	    this.first  = first;
	    this.second = second;
	}
    }


    /* TODO: compare this to utils.Future and such, and maybe finish it
    public static class Errant<D, Tool.Callback<R>> {
	
	final D val;
	final Tool.Callback<R> iface;

	protected Errant( D val, Tool.Callback<R> iface) {

	}


	public returnResult(R) {

    } */
     

} public interface ImageSelector {

    /** Information returned about an open image */
    public class ImageInfo {
	/** Image width and height (pxl) */
	public final int width, height;
	/** Depth / length of stack */
	public final int depth;
	/** Size of pxl in microns (or -1 if not available) */
	public final double micronsPerPxl;
	/** "Name" of the image (or null if not available) */
	public final String name;
	/** ID, free to be set by implementors of ImageInfo */
	public int id;
		
	public ImageInfo( int width, int height, int depth,
	    double micronsPerPxl, String name, int id ) {
	    this.width		= width;
	    this.height		= height;
	    this.depth		= depth;
	    this.micronsPerPxl	= micronsPerPxl;
	    this.name		= name;
	    this.id		= id;
	}

	@Override
	public String toString() { return name; }

	/** Compares images on size and id only.
	 *  ImageInfo's are defined to be the same for
	 *  width, height, depth and id being equal, while
	 *  name and micronsPerPixel may differ */
	@Override
	public boolean equals( Object e ) {
	    // other types (and null): false
	    if (!(e instanceof ImageInfo ))
		return false;
	    // compare w,h,d,id
	    ImageInfo other = (ImageInfo)e;

	    if (    this.width  == other.width
		&&  this.height == other.height
		&&  this.depth	== other.depth
		&&  this.id	== other.id )
	    return true;
	    return false;
	}

    }


    /* see fiji/ImageOpener, currently, the callback-idea
     * to update the list seems to cause the annoying
     * deadlocks in the code... reverting to manual
     * updates of the image list */
    /*
    // Nofication that available images changed 
    public interface Callback {
	// gets called on change 
	public void call();
    }
    // register a callback 
    public void addCallback( Callback b );
    // de-register a callback 
    public void removeCallback( Callback b );
    */

    /** Returns the number of open images */
    public int getOpenImageCount();    
    
    /** Get list of open images */
    public ImageInfo [] getOpenImages();

    /** Returns the image referred to by an ImageInfo, at position pos
     * (or null if not found) */
    public Vec2d.Real getImage( ImageInfo which, int pos );
    
    /** Returns the image stack referred to by an ImageInfo, 
     * (or null if not found) */
    public Vec2d.Real [] getImages( ImageInfo which );


    /** Dummy implementation for testing (GUI) */
    public class Dummy implements ImageSelector {
	final int n;
	public Dummy(int i) { n=i; }
	final String [] names = { "Test Image 0", "Test Image 1",
	    "Some other image (w/o calibration)", 
	    "Very long image name just to test if GUI is wide enough",
	    "Image with 45 planes"};

	
	@Override
	public int getOpenImageCount() { return n; }
	@Override
	public ImageInfo [] getOpenImages() {
	    ImageInfo [] ret = new ImageInfo[ n ];
	    for (int i=0; i<n; i++)
		ret[i] = new ImageInfo(512, 512, idxToZ(i),
		    ((i!=2)?(0.082):(-1)),
		    names[i%names.length], i );
	    return ret;
	}
	@Override
	public Vec2d.Real getImage(ImageInfo w, int pos) {
	    return Vec2d.createReal( 512, 512 );
	}
	@Override
	public Vec2d.Real [] getImages(ImageInfo w) {
	    return Vec2d.createArrayReal( idxToZ(w.id) , 512, 512 );
	}

    
	private int idxToZ(int i) {
	    return (i!=1)?(15*(i+1)):(9*(i+1));
	}

    }




} public interface ImageDisplay {

    /** Set the i'th image to display the content of v.
     *  @param v Vector content to display (implementation must copy content)
     *  @param i Position of image to set
     *	@param label Label of image (may be 'null')
     *	@param m Markers to draw on the image (see {@link ImageDisplay.Marker}
     *	*/
    public abstract void setImage( Vec2d.Real v, int i, String label, Marker ... m); 

    /** Add an image to the end of the stack. 
     *  @param v Vector content to display (implementation must copy content)
     *	@param label Label of image (may be 'null')
     *	@param m Markers to draw on the image (see {@link ImageDisplay.Marker}
     *	*/
    public abstract void addImage( Vec2d.Real v, String label, Marker ... m );
    
    /** TODO: comment this */
    public void addImage( Vec3d.Real v, boolean project, String label, ImageDisplay.Marker ... m) ;

    /** Set that the display should be visible. Implementations (e.g.
     *  output to disk) may ignore this flag. */
    public void display();

    /** Called to tell a display it is not longer needed. Thus, it should close
     * and can drop all references, data etc. */
    public void drop();

    /** Returns the number of images that have been added
     *  to the display. */
    public int getCount();

    /** Register a listener */
    public void addListener( Notify n ); 
	    
    /** De-register a listener */
    public void removeListener( Notify n ); 

    /** Get image width */
    public int width();

    /** Get image height */
    public int height();

    // ------------------------------------------------------------------------


    /** Marker to draw onto an image */
    final class Marker {
	/** Coordinates x-pos, y-pos, width, height */
	public final double x,y,w,h;
	/** Draw an elipse instead of a rectangle */
	public final boolean drawElipse;
	
	/** Create a square marker
	 * @param x x-position center
	 * @param y y-position center
	 * @param w witdh 
	 * @param h height 
	 * @param elipse If true, draw a elipse instead of a rectangle */
	public Marker(double x, double y, double w, double h, boolean elipse) {
	    this.x=x; this.y=y; this.w=w; this.h=h; this.drawElipse = elipse;
	}
    }
    
    /** Position change in an ImageDisplay */
    public interface Notify {
	/** Displayed image changed position */
	public void newPosition( ImageDisplay e, int p );
    }


    /** Creation of ImageDisplays */
    public interface Factory {
	/** Create a new ImageDiplay.
	 *  @param w Width of the display
	 *  @param h Height of the display 
	 *  @param t Title of the display */
	public ImageDisplay create(int w, int h, String t);

    }


} public class OtfProvider {
    
    // --- internal parameters ----

    public enum APPROX_TYPE {
	EXPONENTIAL,
	SPHERICAL,
	NONE
    };

    // vals[band][idx], where idx = cycles / cyclesPerMicron
    private Cplx.Float [][] vals	=  null; 
    private Cplx.Float [][] valsAtt	=  null; 
    private float      [][] valsOnlyAtt =  null; 
    
    // physical units
    private double cyclesPerMicron;
    private double na, lambda, cutOff;
    
    // was the OTF set from estimate? 
    // Are there different bands? If so, how many?
    private boolean isEstimate, isMultiBand; 
    int maxBand=1; 
    double estimateAValue=.3;
    
    // for estimate: sample with how many points?
    // for data: this is overridden by data
    private int samplesLateral = 512;

    // for vector output: physical pixel size
    private int vecSize=-1;
    private double vecCyclesPerMicron=-1;

    // attenuation strength (0..1) and fhwm (in cycles/micron)
    private double attStrength = .99, attFWHM = 1.2;
    private boolean useAttenuation;


    /** For [0..cutoff] normalized to [0..1], return the ideal OTF. OTF of an ideal, 
     *  i.e. abberation-free lens system, resolution-limited by a circular pupil in
     *  the Fourier place,
     *  "Joseph W. Goodman, Introduction to Fourier Optics, 3. edition, page 145".
     *  @param dist Distance to cutoff, in normalized range 0..1. Values outsize of 0..1 return 0.
     */
    public static double valIdealOTF(double dist) {
	if ((dist<0)||(dist>1))
	    return 0;
	return (2/Math.PI)*(Math.acos(dist) - dist*Math.sqrt(1-dist*dist)); 
    }


    // ------ setup / access methods -------

    /** use factory methods instead */
    private OtfProvider() {};


    // TODO: deprecate and finally remove this function version
    public static OtfProvider fromEstimate(double na, double lambda, double a ) {
	return fromEstimate(na,lambda,a,APPROX_TYPE.EXPONENTIAL);
    }

    /** Create a new OTF from a (very basic) estimate.
     *  The curvature factor account for deviation of
     *  real-world OTFs from the theoretical optimum.
     *  a=1 yields an ideal OTF, a=0.2 .. 0.4 are more realistic,
     *  a has no effect on cutoff.
     *  @param na Objectives NA
     *  @param lambda Emission wavelength (nm)
     *  @param a curvature factor, a = [0..1]
     *  @param compType which compensation formula to use
     * */
    public static OtfProvider fromEstimate(double na, double lambda, double a, APPROX_TYPE compType ) {
	if ( (a<0)||(a>1) || (na<0.3) || (na>2.2) || (lambda<300) || (lambda>1500) )
	    throw new IllegalArgumentException("unphysical input parameters");

	OtfProvider ret = new OtfProvider();

	// compute / set parameters
	ret.na = na; ret.lambda = lambda;
	ret.cutOff = 1000 / ( lambda / na / 2 ); 
	ret.cyclesPerMicron = ret.cutOff / ret.samplesLateral ;
	ret.vals	= new Cplx.Float[1][ ret.samplesLateral ];
	ret.valsAtt	= new Cplx.Float[1][ ret.samplesLateral ];
	ret.valsOnlyAtt = new	   float[1][ ret.samplesLateral ];
	ret.isMultiBand = false;
	ret.isEstimate  = true;
	ret.estimateAValue = a;

	// sample some values up to cutoff
	for (int i=0; i< ret.samplesLateral ; i++) {
	    // v: normalize [0..cutoff] -> [0..1]
	    double v = i/(double)ret.samplesLateral ;
	    // get OTF at v, multiply empirical correction for curvature
	    float r =(float)valIdealOTF(v);
	    
	    if ( compType == APPROX_TYPE.EXPONENTIAL ) {
		r *= (float)Math.pow(a,v);
	    } 
	    if ( compType == APPROX_TYPE.SPHERICAL ) {
		r *= (float)(1 - ( (1-a)* (1 - 4*(Math.pow(v-.5,2)))  ));
	    }

	    ret.vals[0][i] = new Cplx.Float( r, 0 );
	}
	
	// initialize attenuation cache
	ret.setAttenuation( ret.attStrength, ret.attFWHM );
	ret.switchAttenuation( false );

	return ret;
    }

    /** Returns a short description (GUI display, etc) */
    public String printState(boolean html) {
	String ret ="";
	if (html)
	    ret = String.format("NA %4.2f, \u03bb<sub>em</sub> %4.0f, ", na, lambda);
	else 
	    ret = String.format("NA %4.2f, lambda %4.0f, ", na, lambda);
	if (isEstimate) 
	    ret+=String.format("(est., a=%4.2f)", estimateAValue);
	else
	    ret+=String.format("(from file)");
	return ret;
    }


    /** Return the OTF cutoff, unit is cycles/micron */
    public double getCutoff() {
	return cutOff;
    }

    /** Get the OTF value at 'cycl'.  
     *  @param band OTF band
     *  @param cycl Position in cycles/micron
     *  @param att  If true, return attenuated value (see {@link #setAttenuation})
     * */
    public Cplx.Float getOtfVal(int band, double cycl, boolean att) {
	// checks
	if ( !this.isMultiBand ) 
	    band=0;
	if ((band >= maxBand)||(band <0))
	    throw new IndexOutOfBoundsException("band idx too high or <0");
	if ( cycl < 0 )
	    throw new IndexOutOfBoundsException("cylc negative!");
	
	// out of support, return 0
	if ( cycl >= cutOff )
	    return Cplx.Float.zero();
	
	final double pos = cycl / cyclesPerMicron;
	if (Math.ceil(pos) >= samplesLateral)
	    return Cplx.Float.zero();

    
	// for now, linear interpolation, could be better with a nice cspline
	int lPos = (int)Math.floor( pos );	
	int hPos = (int)Math.ceil( pos );
	float f = (float)(pos - lPos);
    
	if (att) {
	    Cplx.Float retl = valsAtt[band][lPos].mult( 1-f );
	    Cplx.Float reth = valsAtt[band][hPos].mult( f );
	    return Cplx.add( retl, reth );
	} else {
	    Cplx.Float retl = vals[band][lPos].mult( 1-f );
	    Cplx.Float reth = vals[band][hPos].mult( f );
	    return Cplx.add( retl, reth );
	}

    }
   
    /** Sets pixel size, for output to vectors
     *	@param cyclesPerMicron Pixel size of output vector, in cycles/micron */
    public void setPixelSize( double cyclesPerMicron ) {
	if (cyclesPerMicron<=0)
	    throw new IllegalArgumentException("pxl size must be positive");
	vecCyclesPerMicron=cyclesPerMicron;
    }

    // ------ Attenuation -----

    /** Returns the attenuation value at dist.
     *  @param dist Distance to center in cycles/micron
     *  @param str  Strength of the attenuation
     *  @param fwhm FWHM of Attenuation, in cycles/micron  */
    public static float valAttenuation( final double dist, final double str, final double fwhm) {
	    return (float)(1 - str*Math.exp( -( Math.pow(dist,2) ) / (2*Math.pow( fwhm/2.355 ,2)) ));
    }

    /** Get the attenuation value at 'cycl'. 
     *  @param band OTF band
     *  @param cycl Position in cycles/micron
     * */
    public float getAttVal(int band, double cycl) {
	// checks (TODO: mult-band support for attenuation)
	if ( !this.isMultiBand ) 
	    band=0;
	if ((band >= maxBand)||(band <0))
	    throw new IndexOutOfBoundsException("band idx too high or <0");
	if ( cycl < 0 )
	    throw new IndexOutOfBoundsException("cylc negative!");
	
	// out of support, return 1
	if ( cycl >= cutOff )
	    return 1.f;
	
	double pos = cycl / cyclesPerMicron;
    
	// for now, linear interpolation, could be better with a nice cspline
	int lPos = (int)Math.floor( pos );	
	int hPos = (int)Math.floor( pos );
	float f = (float)(pos - lPos);
    
	float reth = valsOnlyAtt[band][lPos] * ( 1-f );
	float retl = valsOnlyAtt[band][hPos] * f ;
	return retl+reth;
    }
    
    /** Sets the OTF attenuation parameters.
     * Attenuation is also swichted on by this function.
     * @param strength Strength of attenuation, 0..1, usually 0.9 .. 0.99
     * @param fwhm     FWHM of the attenuation, in cycles / micron
     * */
    public void setAttenuation( double strength, double fwhm) {
	
	attStrength = strength; attFWHM = fwhm;
	
	// update cached attenuation values
	for ( int b = 0; b<vals.length; b++)
	for ( int v = 0; v<vals[b].length; v++) {
	    double dist = v * cyclesPerMicron;
	    valsOnlyAtt[b][v] = valAttenuation( dist, attStrength, attFWHM ) ;
	    valsAtt[b][v] = vals[b][v].mult( valsOnlyAtt[b][v] ); 
	}
    }

    /** Like {@link #setAttenuation} called with the current FWHM */
    public void setAttenuationStrength( double strength ) {
	setAttenuation( strength, attFWHM );
    }
    /** Like {@link #setAttenuation} called with the current strength */
    public void setAttenuationFWHM( double fwhm ) {
	setAttenuation( attStrength, fwhm );
    }
    


    /** Set if to apply attenuation to the OTF. */
    public void switchAttenuation( boolean on ) {
	useAttenuation = on;
    }
    
    /** Get if attenuation will be used. */
    public boolean isAttenuate() {
        return useAttenuation;
    }
    
    /** Get attenuation strength */
    public double getAttStr(int band) {
	return this.attStrength;
    }
    
    /** Get attenuation FWHM */
    public double getAttFWHM(int band) {
	return this.attFWHM;
    }

    // ------ applying OTF to vectors -------

    /** Multiplies / outputs OTF to a vector. Quite general function,
     *  some wrappers are provided for conveniece. 
     *  @param vec  Vector to write / multiply to
     *	@param band OTF band 
     *	@param kx OTF center position offset x
     *	@param ky OTF center position offset y
     *  @param useAtt if to use attenuation (independent of how {@link #switchAttenuation} is set)
     *  @param write  if set, vector is overridden instead of multiplied
     *  */
    public void otfToVector( final Vec2d.Cplx vec, final int band, 
	final double kx, final double ky,
	final boolean useAtt, final boolean write ) {
	
	// parameters
	if (vecCyclesPerMicron <=0)
	    throw new IllegalStateException("Vector pixel size not initialized");
	final int w = vec.vectorWidth(), h = vec.vectorHeight();

	// loop output vector
	new SimpleMT.StrPFor(0,h) {
	    public void at(int y) {
		for (int x=0; x<w; x++) {
		    // wrap to coordinates: x in [-w/2,w/2], y in [-h/2, h/2]
		    double xh = (x<w/2)?( x):(x-w);
		    double yh = (y<h/2)?(-y):(h-y);
		    // from these, calculate distance to kx,ky, convert to cycl/microns
		    double rad = MTool.fhypot( xh-kx, yh-ky );
		    double cycl = rad * vecCyclesPerMicron;
		    
		    // over cutoff? just set zero
		    if ( cycl > cutOff ) {
			vec.set(x,y, Cplx.Float.zero());
		    } 
		    
		    // within cutoff?
		    if ( cycl <= cutOff ) {
		    
			// get the OTF value
			Cplx.Float val = getOtfVal(band, cycl, useAtt);

			// multiply to vector or write to vector
			if (!write) {
			    vec.set(x, y, vec.get(x,y).mult( val.conj() ) );
			} else {
			    vec.set(x, y, val );
			}
		    }
		}
	    }
	}; 

    }

    // ------ Applying OTF to vectors ------
    
    /** Multiplied conjugated OTF to a vector.
     *  The desired vector pixel size has to be set (via {@link #setPixelSize}) first.
     *  Attenuation is used if set via {@link #switchAttenuation}.
     *  @param vec  Vector to write to
     *	@param band OTF band */
    public void applyOtf(Vec2d.Cplx vec, final int band) {
	otfToVector( vec, band, 0, 0, useAttenuation, false ) ; 
    }

    /** Multiplied conjugated OTF to a vector.
     *  The desired vector pixel size has to be set (via {@link #setPixelSize}) first.
     *  Attenuation is used if set via {@link #switchAttenuation}.
     *  @param vec  Vector to write to
     *	@param band OTF band 
     *	@param kx OTF center position offset x
     *	@param ky OTF center position offset y */
    public void applyOtf(Vec2d.Cplx vec, final int band, double kx, double ky ) {
	otfToVector( vec, band, kx, ky, useAttenuation, false ) ; 
    }

    /** Create a 2-dimension, radial symmetric vector of the OTF, centered at kx,ky.
     *  The desired vector pixel size has to be set (via {@link #setPixelSize}) first.
     *  Attenuation is not applied. 
     *  @param vec  Vector to write to
     *	@param band OTF band 
     *	@param kx Position / offset kx
     *	@param ky Poistion / offset ky
     *	*/
    public void writeOtfVector(final Vec2d.Cplx vec, final int band, 
	final double kx, final double ky) {
	otfToVector( vec, band, kx, ky, false, true ) ; 
    }
    
    /** Create a 2-dimension, radial symmetric vector of the OTF, centered at kx,ky.
     *  The desired vector pixel size has to be set (via {@link #setPixelSize}) first.
     *  Attenuation is applied if set via {@link #switchAttenuation}. 
     *  @param vec  Vector to write to
     *	@param band OTF band 
     *	@param kx Position / offset kx
     *	@param ky Poistion / offset ky
     *	*/
    public void writeOtfWithAttVector(final Vec2d.Cplx vec, final int band, 
	final double kx, final double ky) {
	otfToVector( vec, band, kx, ky, useAttenuation, true ) ; 
    }




    // ------ Other vectors ------

    /** Creates an apotization vector.
     *  This yields an ideal OTF, with some 'bend' (augmenting medium frequencies),
     *  and a new cutoff value. Usually, cutoff is set to 2x original cutoff.
     *  Bend is used as "apo = idealOtf^bend"
     *  @param vec  Vector to write to
     *  @param bend Bend to augment medium frequencies, 0..1, 1 for ideal OTF
     *  @param cutOff Cutoff, as factor to the OTF cutoff (so e.g. 2 for 2x lateral improvement) */
    public void writeApoVector(final Vec2d.Cplx vec, final double bend, final double cutOff ) {
	if (vecCyclesPerMicron <=0)
	    throw new IllegalStateException("Vector pixel size not initialized");
	final int w = vec.vectorWidth(), h = vec.vectorHeight();

	//for (int y=0; y<h; y++)
	new SimpleMT.PFor(0,h) {
	    public void at(int y) {
		for (int x=0; x<w; x++) {
		    // wrap to coordinates: x in [-w/2,w/2], y in [-h/2, h/2]
		    double xh = (x<w/2)?( x):(x-w);
		    double yh = (y<h/2)?(-y):(h-y);
		    // from these, calculate distance to 0, convert to phys. units
		    double rad = MTool.fhypot( xh, yh );
		    double cycl = rad * vecCyclesPerMicron;
		    // calculate fraction of cutoff, get idealOTF, augment with 'bend'
		    double frac = cycl / (getCutoff()*cutOff);
		    double val  = Math.pow( valIdealOTF( frac ), bend );
		    // set output to that value	
		    vec.set(x,y, new Cplx.Float((float)val));
		}
	    }
	}; 
    }

    /** Creates an attenuation vector (for optical sectioning).
     *  Pixel size of output has to be set via {@link #setPixelSize}.
     *  @param out  The vector that will be multiplied
     *  @param str  Strength of the attenuation
     *  @param fwhm FWHM of Attenuation, in cycles/micron 
     *  @param kx   x-pos of the attenuation
     *  @param ky   y-pos of the attenuation
     *  */
    public void writeAttenuationVector(final Vec2d.Real out, final double str, final double fwhm, 
	final double kx, final double ky) {
    
	if (vecCyclesPerMicron <=0)
	    throw new IllegalStateException("Vector pixel size not initialized");
       
	final int w = out.vectorWidth();
	final int h = out.vectorHeight();

	//for (int y= 0;y<h ;y++)
	new SimpleMT.PFor(0,h) {
	    public void at(int y) {
		for (int x= 0;x<w ;x++) {
		
		    // wrap to coordinates: x in [-w/2,w/2], y in [-h/2, h/2]
		    double xh = (x<w/2)?( x):(x-w);
		    double yh = (y<h/2)?(-y):(h-y);

		    // from these, calculate distance to kx,ky, convert to cycl/microns
		    double rad = MTool.fhypot( xh-kx, yh-ky );
		    double cycl = rad * vecCyclesPerMicron;
		    
		    out.set(x,y, valAttenuation(cycl,str,fwhm) );
		}
	    }
	};
    }

    // ------ Mask OTF region ------
   
    /** Set components not covered by OTF support (see {@link #getCutoff}) to zero.
     *  @param vec The vector to multiply with the conj. OTF
     *  @param kx   x-pos of the attenuation
     *  @param ky   y-pos of the attenuation
     * */
    public void maskOtf(final Vec2d.Cplx vec, final double kx, final double ky) {
	
	final int w = vec.vectorWidth(), h = vec.vectorHeight();
	final double otfSupport = getCutoff();

	//for (int y=0; y<h; y++) {{
	new SimpleMT.PFor(0, h) {
	    public void at(int y) {
		for (int x=0; x<w; x++) {
		    // wrap to coordinates: x in [-w,w], y in [-h, h]
		    double xh = (x<w/2)?( x):(x-w);
		    double yh = (y<h/2)?(-y):(h-y);
		    // from these, calculate distance to (kx,ky), convert to cycl/microns
		    double rad = MTool.fhypot( xh-kx, yh-ky ) * vecCyclesPerMicron;
		    // if outside of support, set zero
		    if (rad>otfSupport)
			vec.set(x,y, Cplx.Float.zero());
		}
	    }
	}; 

    }

    // ------ Load / Save operations ------

    /** Create an OTF stored in a config.
     *	@param cfg The config to load from
     *  */
    public static OtfProvider loadFromConfig( Conf cfg ) 
	throws Conf.EntryNotFoundException {

	Conf.Folder fld = cfg.r().cd("otf2d");
	OtfProvider ret = null;
	boolean estimate = (!fld.contains("data"));
	
	// Initialize as estimate
	if (estimate) {
	    ret = fromEstimate(
		fld.getDbl("NA").val(),
		fld.getInt("emission").val(),
		fld.getDbl("a-estimate").val()
		);
	}
	
	// Initialize from data	
	if (!estimate) {
	
	    // main parameters
	    ret = new OtfProvider();
	    ret.na	    = fld.getDbl("NA").val();
	    ret.lambda  = fld.getInt("emission").val();
	    ret.cutOff  = 1000 / (ret.lambda / ret.na /2);
	    ret.isEstimate = false;
	    
	    // copy parameters
	    Conf.Folder data = fld.cd("data");
	    ret.maxBand = data.getInt("bands").val();
	    ret.isMultiBand = (ret.maxBand>1);
	    ret.isEstimate  = false;
	    ret.samplesLateral =data.getInt("samples").val(); 

	    ret.cyclesPerMicron = data.getDbl("cycles").val();
	    
	    // init bands
	    ret.vals	= new Cplx.Float[ret.maxBand][ ret.samplesLateral ];
	    ret.valsAtt	= new Cplx.Float[ret.maxBand][ ret.samplesLateral ];
	    ret.valsOnlyAtt = new  float[ret.maxBand][ ret.samplesLateral ];

	    // copy bands
	    for (int b=0; b<ret.maxBand; b++) {
		byte  [] bytes = data.getData(String.format("band-%d",b)).val();
		float [] val   = Conf.fromByte( bytes );
		
		for (int i=0; i< ret.samplesLateral ; i++) {
		    ret.vals[b][i] = new Cplx.Float( val[2*i], val[2*i+1]);	    
		    
		}
	    }
	
	    ret.setAttenuation( ret.attStrength, ret.attFWHM );
	    ret.switchAttenuation( false );
	}

	// attenuation, if set
	if ( fld.contains("attenuation") ) {
	    Conf.Folder att = fld.cd("attenuation");
	    ret.setAttenuation( att.getDbl("strength").val(), att.getDbl("FWHM").val());
	    ret.switchAttenuation( true );
	}
	
	return ret;
    }

    /** Write out an OTF to a config */
    public void saveConfig( Conf cfg ) {

	Conf.Folder fld  = cfg.r().mk("otf2d");
	fld.newDbl("NA").setVal( this.na );
	fld.newInt("emission").setVal( (int)this.lambda );
	
	// if this is an estimate
	if ( this.isEstimate ) {
	    fld.newDbl("a-estimate").setVal( this.estimateAValue );
	}

	// if this was read from file
	if ( !this.isEstimate ) {
	    Conf.Folder data = fld.mk("data");
	    data.newInt("bands").setVal( this.maxBand );
	    data.newInt("samples").setVal( this.samplesLateral );
	    data.newDbl("cycles").setVal( this.cyclesPerMicron );
	    
	    // store bands
	    for (int b=0; b< this.maxBand; b++) {
		float [] tmp = new float[ this.samplesLateral * 2 ];
		for ( int i=0; i<tmp.length/2; i++) {
		    tmp[2*i+0] = this.vals[b][i].re;
		    tmp[2*i+1] = this.vals[b][i].im;
		}
		byte [] outData = Conf.toByte( tmp );
		data.newData(String.format("band-%d",b)).setVal( outData );
	    }
	}

	// attenuation
	if ( this.useAttenuation ) {
	    Conf.Folder att = fld.mk("attenuation");
	    att.newDbl("strength").setVal( attStrength ); 
	    att.newDbl("FWHM").setVal( attFWHM ); 
	}




    }





    // ------ Testing ------

    /** For testing, outputs some OTF to stdout */
    public static void main( String args [] ) throws Exception {

	if (args.length==0) {
	    System.out.println("Use: t - Timing, o - Output, i - Input ");
	    return;
	}

	OtfProvider otf = OtfProvider.fromEstimate( 1.4, 515, .241);
	
	// output
	if (args[0].equals("o")) {

	    for (float cycl = 0; cycl < 6; cycl +=0.025 )
		System.out.println(String.format(" %5.3f %6.4f A", 
		    cycl, otf.getOtfVal(0, cycl, false).re ));
	    for (int i=0;i<512;i++)
		System.out.println(String.format(" %d %6.4f B",  
		    i, otf.getOtfVal(0, i/(512*0.082), false ).re )); 
	    
	    // save as test 
	    if ( args.length >= 2 ) {
		Conf cfgOut = new Conf("fairsim");
		otf.saveConfig(cfgOut);
		cfgOut.saveFile( args[1]);
	    }
	}
	
	// timing
	if (args[0].equals("t")) {

	    otf.setPixelSize(0.023);
	    //otf.switchAttenuation( true );
	    Vec2d.Cplx [] test = Vec2d.createArrayCplx(15,512,512);

	    for (int loop=0; loop<5; loop++) {

		// Test1: Create the OTF once, cache, apply from cache
		Tool.Timer t1 = Tool.getTimer();
		Vec2d.Cplx otfCache = Vec2d.createCplx( 512,512 );
		otf.writeOtfVector( otfCache, 0, 20, 50);
		t1.start();
		for (int i=0;i<15; i++) 
		    test[i].timesConj( otfCache );
		t1.stop();
		
		// Test1: Apply OTF on the fly
		Tool.Timer t2 = Tool.getTimer();
		t2.start();
		for (int i=0;i<15; i++) 
		    otf.applyOtf( test[i], 0 ) ;
		t2.stop();
		
		Tool.trace("OTF cached: "+t1);
		Tool.trace("OTF on-fly: "+t2);
	    }
	}
	
	// input
	if (args[0].equals("i")) {

	    Conf cfg = Conf.loadFile( args[1] );

	    OtfProvider otfl = OtfProvider.loadFromConfig( cfg ); 

	    for (float cycl = 0; cycl < 6; cycl +=0.025 )
		System.out.println(String.format(" %5.3f %6.4f A", 
		    cycl, otfl.getOtfVal(0, cycl, false).re ));
	    
	    // save as test 
	    if ( args.length >= 3 ) {
		Conf cfgOut = new Conf("fairsim");
		otfl.saveConfig(cfgOut);
		cfgOut.saveFile( args[2]);
	    }

	}



	Tool.shutdown();
    }

} public class SimParam implements Vec2d.Size {

    public enum FilterStyle {
	Wiener("Wiener filter"),	
	RLin("RL on input"), 
	RLout("RL on output"),
	RLboth("RL on both");

	final String name;
	FilterStyle(String a) { name = a; };

	@Override public String toString() { return name; };
    }

    /** number of angles/pattern directions */
    final protected int nrDirs;	
    final protected int nrBands;	
    final protected int nrPhases;
    final protected Dir [] directions ;

    // Image parameters
    private int imgSize=-1;			    // image size in pxl
    private double micronsPerPixel=-1;		    // spatial extent of pixel
    private double cyclesPerMicron=-1;		    // freq extent of FFT pixel

    private int stackSize=-1;			    // z-stack length in #slices
    private double micronsPerSlice=-1;		    // spatial extent of the z-slice stage shift
    private double cyclesPerMicronInZ=-1;	    // freq extent of FFT in z-slices

    private IMGSEQ imgSeq = IMGSEQ.PAZ;		    // order of images in input
    private CLIPSCALE clipScaleMode 
	= CLIPSCALE.NONE;			    // clip&scale of output


    private FilterStyle filterStyle = 
		FilterStyle.Wiener;		    // which filter to use

    private double wienerFilterParameter = 0.05;    // Wiener filter parameter
    private double apoCutOff = 2;		    // Apo cutoff parameter
    private double apoBend   = 0.9;		    // Apo bend parameter

    private int rlIterations = 5;		    // number of Richardson-Lucy iterations
    

    double modLowLimit = 0.4, modHighLimit = 1.1, defaultModulation = 0.65;

    // Transfer function, OTF attenuation, Apotization
    private OtfProvider   currentOtf2D=null;

    private long runtimeTimestamp = 0;


    /** Use factory method {@link #create} to obtain object */
    protected SimParam(int bands, int dirs, int phases, boolean threeD) { 
	if (dirs<1)
	    throw new RuntimeException("dirs < 1, not useful");
	
	nrDirs   = dirs;
	nrBands  = bands;
	nrPhases = phases;

	directions = new Dir[ nrDirs ];
	for (int i=0; i<nrDirs; i++)
	    directions[i] = new Dir( bands, phases, i );
	
    }
    
    /** New set of SIM reconstruction parameters for 2D reconstruction. 
     *	@param bands  Number of bands
     *	@param dirs   Number of pattern orientations
     *	@param phases Number of phases
     * */
    public static SimParam create(int bands, int dirs, int phases ) {
	return new SimParam( bands, dirs, phases, false);
    }
    
    /** New set of SIM reconstruction parameters for 2D reconstruction. 
     *	@param bands  Number of bands
     *	@param dirs   Number of pattern orientations
     *	@param phases Number of phases
     *	@param size Size in pixels (width or height) of raw images
     *	@param micronsPerPxl Size of a pixel in micrometers 
     * */
    public static SimParam create(int bands, int dirs, int phases, 
	int size, double micronsPerPxl, OtfProvider otf ) {
	SimParam sp = new SimParam( bands, dirs, phases, false);
	sp.setPxlSize( size, micronsPerPxl );
	sp.otf( otf );
	return sp;
    }
   
    
    /** Get parameter subset for pattern direction 'i' */
    public Dir dir(int i) {
	return directions[i];
    }

    /** Get the number of directions */
    public int nrDir() {
	return nrDirs;
    }
    
    /** Get the number of bands */
    public int nrBand() {
	return nrBands;
    }

    /** Get cycles / micron pxl size */
    public double pxlSizeCyclesMicron() {
	return cyclesPerMicron;
    }
    
   
    
    /** Get the image ordering */
    public IMGSEQ getImgSeq() {
	return imgSeq;
    }
    /** Set the image ordering */
    public SimParam setImgSeq( IMGSEQ i) {
	imgSeq = i;
	return this;
    }

    /** Get the image ordering */
    public CLIPSCALE getClipScale() {
	return clipScaleMode;
    }
    
    /** Set the image ordering */
    public SimParam setClipScale( CLIPSCALE i) {
	clipScaleMode = i;
	return this;
    }

    /** Get the number of images in one z plane */
    public int getImgPerZ() {
	int ret=0;
	for ( Dir d : directions )
	    ret+= d.nrPha();
	return ret;
    }
	
    /** Get number of phases */
    public int nrPha() {
        return nrPhases;
    }


    /** Set the image size
     *  @param pxl Number of pixels
     *  @param microns Number of microns per pxl */
    public SimParam setPxlSize( int pxl, double microns ) {
	imgSize  = pxl;
	micronsPerPixel = microns;
	cyclesPerMicron = 1/(pxl*microns);
	this.otf( currentOtf2D );	// propagate size to OTF
	return this;
    }



    /** Set the filter type to use */
    public void setFilterStyle( FilterStyle s ) {
	filterStyle = s;
    }

    /** Get the filter type to use */
    public FilterStyle getFilterStyle() {
	return filterStyle;
    }

    /** Determine if Wiener-filtering is used on output data*/
    public boolean useWienerFilter() {
	return ( filterStyle == FilterStyle.Wiener ||
		 filterStyle == FilterStyle.RLin );
    }



    /** Determine if RL-filtering is used on input data*/
    public boolean useRLonInput() {
	return ( filterStyle == FilterStyle.RLboth || 
		 filterStyle == FilterStyle.RLin   );
    }
    
    /** Determine if RL-filtering is used on output data*/
    public boolean useRLonOutput() {
	return ( filterStyle == FilterStyle.RLboth || 
		 filterStyle == FilterStyle.RLout   );
    }


    /** Set the number of RL iterations */
    public void setRLiterations( int n ) {
	rlIterations = n;
    }

    /** Get the number of RL iterations */
    public int getRLiterations() {
	return rlIterations;
    }


    /** Set the Wiener Filter parameter */
    public SimParam setWienerFilter( double wf ) {
	wienerFilterParameter = wf;
	return this;
    }
    
    /** Get the Wiener Filter parameter */
    public double getWienerFilter() {
	return wienerFilterParameter;
    }
    
    /** Set the APO cutoff factor (cutoff in relation to OTF cutoff) */
    public SimParam setApoCutoff( double af ) {
	apoCutOff = af;
	return this;
    }
   

    /** Set the APO bend parameter (curvature of the APO) */
    public SimParam setApoBend( double ab ) {
	apoBend = ab;
	return this;
    }
   


    /** Get the APO cutoff factor */
    public double getApoCutoff() {
	return apoCutOff;
    }
    
    /** Get the APO cutoff factor */
    public double getApoBend() {
	return apoBend;
    }



    /** Get the current otf */
    public OtfProvider otf() {
	return currentOtf2D;
    }
    

    /** Set a new OTF */
    public void otf(OtfProvider otf) {

	if (otf!=null) { 
	    currentOtf2D=otf;
	    otf.setPixelSize( cyclesPerMicron );
	}
    }


    // ----------------------------------------------------------------------------------

    /** Parameters for each pattern orientation. */
    public class Dir {
	
	private final int nrBands;  // how many bands
	private final int nrPhases; // how many phases
	private double pX, pY;	    // shift vector (band1)
	private double phaOff;	    // global phase offsets
	private double [] phases;   // phases
	private double [] modul;    // modulation 
	private int thisBand;	    // which band is this
	private boolean hasIndividualPhases;	// if non-equidist. phases are set

	private double angleIntensityFactor=1;   // the intensity factor to apply to this angle
	private double [] phaseIntensityFactors; // the intensity factors to apply phase-by-phase


	Dir(int nBands, int nPha, int thisBand) {
	    
	    if (nBands<2) 
		throw new RuntimeException("bands < 2, not useful");
	    if (nPha < (nBands*2) -1)
		throw new RuntimeException("not enough phases for bands:"+nBands+" "+nPha);
	    
	    nrBands = nBands;
	    nrPhases = nPha;
	    phases = new double[ nrPhases ];
	    resetPhases();  // inits the phases to equidistant
	    modul=new double[ nrBands ];
	    this.thisBand = thisBand;

	    for (int i=0;i<nrBands;i++)
		modul[i] = 1.0;
	
		phaseIntensityFactors = new double[nrPhases];
	    for (int i=0;i<nrPhases;i++)
		phaseIntensityFactors[i] = 1.0;

	}

	void failBand(int i) {
	    if ((i<0)||(i>=nrBands))
		throw new RuntimeException("Wrong band index");
	}

	/** Return number of bands */
	public int nrBand() { return nrBands; }
	
	/** Return number of components (2*bands-1) */
	public int nrComp() { return nrBands*2-1; }


	// --- Phases ---

	/** Set a global phase offset.  */
	public void setPhaOff( double pha ) {
	    phaOff = pha;
	}

	/** Set individual phases.
	 *  @param pha Phases to set
	 *  @param reset If true, resets the phase offset to 0
	 * */
	public void setPhases( double [] pha, boolean reset ) {
	    if (pha.length != nrPhases )
		throw new RuntimeException("Length mismatch");
	    if (reset)
		phaOff = 0;
	    for (int i=0; i<nrPhases; i++)
		phases[i] = pha[i];
	    hasIndividualPhases = true;
	}

	/** Set one individual phase.
	 * @param pha Phase to set
	 * @param phaIdx Index of phase to set
	 */
	public void setSinglePhase( double pha, int phaIdx ) {
	    phases[phaIdx] = pha;
	    hasIndividualPhases = true;
	}


	/** Reset individual phases back to equidistant. */
	public void resetPhases() {
	    resetPhases(1.);
	}

	/** Reset individual phases back to equidistant. */
	public void resetPhases(double multiplier) {
	    for (int i=0; i<nrPhases; i++) 
		phases[i] = ( 2*Math.PI*multiplier / nrPhases ) * i;
	    hasIndividualPhases = false;
	}
	
	/** Get absolute phases. I.e. phases with the offset applied.
	 *  This can directly be passed as 
	 *  'phases' to {@link BandSeparation#separateBands}. */
	public double [] getPhases() {
	    double [] ret = new double[nrPhases];
	    for (int i=0; i<nrPhases; i++)
		ret[i] = phases[i] + phaOff;
	    return ret;
	}

	/** Get number of phases */
	public int nrPha() {
	    return nrPhases;
	}

	/** Get phase #i, without offset applied. */
	public double getPhase(int i) {
	    return phases[i];
	}
	/** Get the phase offset */
	public double getPhaOff(){
	    return phaOff;
	}


	// --- Shifts ---

	/** Set shift px, py (from estimate for highest band). 
	 *  Lower bands are set equidistant in between. 
	 *  @param cor Either an array (x[0],y[1]) or 2 numbers: x,y
	 *  */
	public void setPxPy(double inPx, double inPy) { 
	    pX = inPx/(nrBands-1); 
	    pY = inPy/(nrBands-1); 
	}

	/** Return x for shift of band n. 
	 *  Bands count from 0.*/
	public double px(int band) {
	    failBand(band);
	    return pX*band;
	}
	
	/** Return y for shift of band n. 
	 *  Bands count from 0.*/
	public double py(int band) {
	    failBand(band);
	    return pY*band;
	}
    
	/** Return shift of band n. 
	 *  Bands count from 0.*/
	public double [] getPxPy(int band) {
	    failBand(band);
	    return new double [] { pX*band, pY*band };
	}

	/** Return angle of shift of band n (atan2). 
	 *  Bands count from 0.*/
	public double getPxPyAngle(int band) {
	    failBand(band);
	    return Math.atan2( pY*band, pX*band );
	}
	
	/** Return length of shift of band n (hypot). 
	 *  Bands count from 0.*/
	public double getPxPyLen(int band) {
	    failBand(band);
	    return Math.hypot( pY*band, pX*band );
	}

	/** Compute the estimated resolution improvement */
	public double getEstResImprovement() {

	    if (currentOtf2D == null) {
		throw new RuntimeException("No OTF set for SIM parameter instance");
	    }

	    double shift	= getPxPyLen( nrBands -1) * cyclesPerMicron;
	    double otfCutoff    = currentOtf2D.getCutoff();

	    return (shift+otfCutoff) / otfCutoff;

	}



	// --- modulation ---
	/** Set the modulation of band n. */
	public boolean setModulation(int b, double m) {
	    failBand(b);
	    modul[b]=m;
	    return ((m>modLowLimit)&(m<modHighLimit));
	}

	/** Get modulations */
	public double [] getRawModulations() {
	    double [] ret = new double [nrBands];
	    System.arraycopy( modul, 0, ret, 0, nrBands );
	    return ret;
	}
     
	/** Get modulations (limited to 'reasonable' range) */
	public double [] getModulations() {
	    double [] ret = new double [nrBands];
	    System.arraycopy( modul, 0, ret, 0, nrBands );
	    // if any modulation is outside of limits, return limit
	    for (int b=0; b< nrBands; b++) {
		if (ret[b]<modLowLimit)  ret[b] = defaultModulation;
		if (ret[b]>modHighLimit) ret[b] = defaultModulation;
	    }
	    return ret;
	}
     


	// -- correction factors ---
	
	/** Set an intensity correction factor for this angle */
	public void setAngleIntensityFactor( double f ) {
	    angleIntensityFactor = f;
	}

	/** Set a phase intensity correction factor for this angle and phase i */
	public void setPhaseIntensityFactor( int p, double f ) {
	    phaseIntensityFactors[p] = f;
	}

	/** Get the angle intensity factor */
	public double getAngleIntensityFactor() {
	    return angleIntensityFactor;
	}

	/** Get a phase intensity factor */
	public double getPhaseIntensityFactor(int p) {
	    return phaseIntensityFactors[p];
	}

	/** Returns a combined correction quotient: 1/(angleIntesity*phaseIntensity) */
	public double getIntensityQuotient(int p) {
	    return 1/(angleIntensityFactor*phaseIntensityFactors[p]);
	}



    }

    /** check if index belongs to a band */
    void failBand(int i) {
	if ((i<0)||(i>=nrBands))
	    throw new RuntimeException("Wrong band index");
    }

    // ----------------------------------------------------------------------------------
    
    /** Save the parameters to a configuration folder.
	Creates sub-folders 'sim-param'  */
    public void saveConfig( Conf.Folder cfg ) {
    
	Conf.Folder fd = cfg.mk("sim-param");

	fd.newInt("nr-angles").setVal(nrDirs);
	fd.newInt("nr-bands").setVal(nrBands);
	fd.newInt("nr-phases").setVal(nrPhases);
	fd.newStr("img-seq").val( imgSeq.name());
	fd.newInt("img-size-pxl").setVal(imgSize);
	fd.newDbl("microns-per-pxl").setVal(micronsPerPixel);
	fd.newDbl("wiener-parameter").setVal( wienerFilterParameter );
	fd.newDbl("apodization-cutoff").setVal( apoCutOff );
	fd.newDbl("apodization-bend").setVal( apoBend );
    
	for ( int d=0; d < nrDirs; d++ ) {
	    Conf.Folder df = fd.mk(String.format("dir-%d",d));
	    SimParam.Dir dir = this.dir(d);

	    //df.newInt("nr-bands").setVal(dir.nrBands);
	    //df.newInt("nr-phases").setVal(dir.nrPhases);
	    df.newDbl("shift").setVal( new double [] 
		{dir.pX *(dir.nrBands-1) , dir.pY * (dir.nrBands-1)} );
	    df.newDbl("phase-offset").setVal( dir.phaOff );
	    df.newDbl("modulations").setVal( dir.modul );
	    if ( dir.hasIndividualPhases )
		df.newDbl("phases").setVal( dir.phases );

	}
    }

    /** Initialize a parameter set from a configuraion folder */
    public static SimParam loadConfig( Conf.Folder cfg ) 
	throws Conf.EntryNotFoundException {
	
	Conf.Folder fd = cfg.cd("sim-param");
	
	// basics
	int nDirs = fd.getInt("nr-angles").val();
	int nBand = fd.getInt("nr-bands").val();
	int nPhas = fd.getInt("nr-phases").val();
	SimParam ret = new SimParam( nBand, nDirs, nPhas, false );
	
	// pixel size
	ret.setPxlSize( fd.getInt("img-size-pxl").val(),
			fd.getDbl("microns-per-pxl").val() );

	// image type
	IMGSEQ isq = IMGSEQ.fromName( fd.getStr("img-seq").val() );
	ret.setImgSeq(isq);
	
	// filter settings
	ret.setWienerFilter( fd.getDbl("wiener-parameter").val() );
	ret.setApoCutoff( fd.getDbl("apodization-cutoff").val() );
	
	// optional parameters that might have been stored
	if ( fd.contains("apodization-bend")) {
	    ret.setApoBend( fd.getDbl("apodization-cutoff").val());
	};


	// for each pattern direction ...
	for ( int d=0; d < ret.nrDirs; d++ ) {
	    Conf.Folder df = fd.cd(String.format("dir-%d",d));
	    SimParam.Dir dir = ret.dir(d);
	    
	    dir.setPxPy( df.getDbl("shift").vals()[0], df.getDbl("shift").vals()[1] );
	    dir.setPhaOff( df.getDbl("phase-offset").val());
	    double [] mod = df.getDbl("modulations").vals();
	    for ( int i=0; i< mod.length ; i++)
		dir.setModulation( i, mod[i]);
	    if ( df.contains("phases"))
		dir.setPhases( 	df.getDbl("phases").vals(), false );

	}



	// return the result
	return ret;
    }





   
    // ----------------------------------------------------------------------------------

    public static enum CLIPSCALE {
	BOTH("clip&scale"),
	CLIP("clip zeros"),
	NONE("raw values");

	final String name;

	CLIPSCALE(String n) {
	    name=n;
	}

	@Override 
	public String toString() {
	    return name;
	}

    }

    public static enum IMGSEQ {
	PAZ("p,a,z (def)"), 
	PZA("p,z,a (OMX)"), 
	ZAP("z,a,p (Zeiss)");
	
	final String name;
    
	IMGSEQ(String n) {
	    name=n;
	}
	
	@Override 
	public String toString() {
	    return name;
	}

	/** Calculate position in (raw) z stack.
	 *  @param ang Index of angle
	 *  @param pha Index of phase
	 *  @param z   Index of z
	 *  @param angMax   Number of angles
	 *  @param phaMax   Number of phases
	 *  @param zMax	    Number of z slices */
	public int calcPos( int ang, int pha, int z, int angMax, int phaMax, int zMax ) {
	    if (( ang>=angMax ) || ( pha >= phaMax ) || ( z >= zMax ) || 
		( ang <0 ) || ( pha < 0) || ( z < 0 ))
		throw new RuntimeException("Parameter wrong!");
	    
	    switch(this) {
		case PZA: return ( pha +   z*phaMax  + ang*phaMax*zMax );
		case ZAP: return (  z  + ang*zMax    + pha*zMax*angMax );
		case PAZ: 
		default:
			  return ( pha + ang*phaMax  + z*phaMax*angMax );
	    }
	}
	
	/** Calculate position in (raw) z stack.
	 *  @param ang Index of angle
	 *  @param pha Index of phase
	 *  @param z   Index of z
	 *  @param t   Index of time
	 *  @param angMax   Number of angles
	 *  @param phaMax   Number of phases
	 *  @param zMax	    Number of z slices */
	public int calcPosWithTime( int ang, int pha, int z, int t, int angMax, int phaMax, int zMax ) {
	    if (( ang>=angMax ) || ( pha >= phaMax ) || ( z >= zMax ) || 
		( ang <0 ) || ( pha < 0) || ( z < 0 ))
		throw new RuntimeException("Parameter wrong!");
	    
	    switch(this) {
		case PZA: return ( pha +   z*phaMax  + ang*phaMax*zMax + 
		    phaMax*zMax*angMax * t );
		case ZAP: return (  z  + ang*zMax    + pha*zMax*angMax +
		    phaMax*zMax*angMax * t );
		case PAZ: 
		default:
			  return ( pha + ang*phaMax  + z*phaMax*angMax +
		    phaMax*zMax*angMax * t );
	    }
	}

	/** Return a default SIM parameter set */
	public SimParam getParam() {
	    switch(this) {
		case PZA: return SimParam.create(3,3,5,512,0.082,null);
		case ZAP: return SimParam.create(3,5,5,1024,0.064,null);
		default:  return null;
	    }
	}
	
	/** Return the ENUM from name */
	public static IMGSEQ fromName( String n ){
	    for ( IMGSEQ i : IMGSEQ.values() )
		if ( n.trim().equals( i.name()) )
		    return i;
	    return null;
	}



    }

    // ----------------------------------------------------------------------------------

    /** update the internal timestamp to signal changes, return the timestamp */
    public long signalRuntimeChange() {
	runtimeTimestamp = System.currentTimeMillis();
	return runtimeTimestamp;
    }

    /** returns the last time 'signalRuntimeChange' was called, as
     * System.currentTimeMillis */
    public long getRuntimeTimestamp() {
	return runtimeTimestamp;
    }
    
    /** returns true if the timestamp provided is older than the stamp generated
     * by the last call to singalRuntimeChange. */ 
    public boolean compareRuntimeTimestamp( long timestamp ) {
	return (( timestamp - runtimeTimestamp ) < 0);
    }


    // ----------------------------------------------------------------------------------
    
    /** Return width */
    @Override public int vectorWidth() { return imgSize; }
    
    /** Return height */
    @Override public int vectorHeight() { return imgSize; }
    

    /** Returns a multi-line, human-readable output of parameters */
    public String prettyPrint(boolean phaInDeg) {
	final double pf = ((phaInDeg)?(180/Math.PI):(1));
	String ret = "#--- SIM parameter summary ---\n";
	ret += String.format("# pxl %7.5f microns (freq pxl %8.6f cycles/micron)\n ----\n",
	     micronsPerPixel, cyclesPerMicron);
	for ( Dir d : directions ) {
	    int b = d.nrBand()-1; 
	    ret += String.format("Nr bands: %d, shift of outer-most band:\n", b);
	    ret += String.format("px: %6.3f py: %6.3f (len %6.3f) \n",
		d.px(b), d.py(b), Math.sqrt(d.px(b)*d.px(b)+d.py(b)*d.py(b)));
	    ret += String.format("phases (%s) [ ",((phaInDeg)?("deg"):("rad")) );
	    for (double p : d.getPhases() )
		ret+=String.format("%6.3f ", p*pf);
	    ret += " ]\n Modulation: [ ";
	    for (double m: d.getModulations() ) {
		ret+=String.format("%6.3f ",m);
	    }
	    ret += " ]\n ----\n";
	}
	return ret;
    }
    
    // ----------------------------------------------------------------------------------

    /** For testing input / output to CONF */
    public static void main( String args [] ) throws Exception {

	if (args.length<2) {
	    System.out.println("Use: i - Input, o - Output [filename] ");
	    return;
	}

	SimParam foo = SimParam.create(3,4,5,512,0.082, null);

	// output
	if (args[0].equals("o")) {
	    Conf cfg = new Conf("fairsim");
	    foo.saveConfig( cfg.r() );
	    cfg.saveFile(args[1]);
	}
	
	// input
	if (args[0].equals("i")) {
	    Conf cfg = Conf.loadFile(args[1]);
	    SimParam bar = loadConfig(cfg.r());
	    System.out.println( bar.prettyPrint(true));
	}
	

	Tool.shutdown();
    }

} public class SimUtils {


    /** Places a 2d freq space vector into a vector of double size.
     *  @param inV Input vector
     *  @param outV Output vector, of doubled size
     *  */
    static public void placeFreq( Vec2d.Cplx inV, Vec2d.Cplx outV ) {
	outV.pasteFreq(inV);
    }

    /** Moves freq space data to kx, ky with subpixel precision, by phase
     *  multiplication in real space. Basically, just ifft, multiply phases, fft...
     *  @param inV Vector to move
     *  @param kx x-coord to move to
     *  @param ky y-coord to move to
     *  */
    static public void fourierShift( Vec2d.Cplx inV,
	final double kx, final double ky ) {

	// tranform to real space
	Transforms.fft2d( inV , true );
	Transforms.timesShiftVector( inV, kx, -ky );
	Transforms.fft2d( inV , false );
    }

    /** Moves freq space data to kx, ky with subpixel precision, by phase
     *  multiplication in real space. 
     *  In standard mode, input vector is copied into output, then moved
     *  by Fourier shift on the output vector.
     *  In fast mode, the subpixel move is carried out on the input vector,
     *  result than is copied to the correct (multiple integer) position in
     *  output vector.
     *
     *  @param inV Input vector, w x h (vector will be changed in fast mode!)
     *  @param outV Output vector, 2*w x 2*h
     *  @param kx x-coord to move to
     *  @param ky y-coord to move to
     *  */
    static public void pasteAndFourierShift( Vec2d.Cplx inV, Vec2d.Cplx outV,
	final double kx, final double ky, boolean doFastShift ) {

	if (!doFastShift) {
	    // standard, slow shift (two FFTs on 2w,2h )
	    outV.pasteFreq( inV, 0, 0);
	    outV.fft2d( true );
	    outV.fourierShift( kx, -ky );
	    outV.fft2d( false );
	} else {

	    // do the subpixel part on input
	    double xFrac = kx - Math.floor(kx);
	    double yFrac = -ky - Math.floor(-ky);
	    inV.fft2d( true );
	    inV.fourierShift( xFrac, yFrac);
	    inV.fft2d( false );
	    
	    // do the integer shift in paste
	    int xInt = (int)Math.floor(  kx );
	    int yInt = (int)Math.floor( -ky );
	    outV.pasteFreq( inV, xInt, yInt );

	    Tool.trace(String.format(" %7.4f %7.4f %7.4f %7.4f" ,
		xFrac, yFrac, (double)xInt, (double)yInt));
	}
    }


    /** Return the power spectrum (as a new vector) of 'in' */
    public static Vec2d.Real pwSpec(Vec2d.Cplx in , boolean swap) {
	Vec2d.Real pw = Vec2d.createReal(in);
	Transforms.computePowerSpectrum( in, pw, swap);
	return pw;
    }
    
    /** Return the power spectrum (as a new vector) of 'in' */
    public static Vec2d.Real pwSpec(Vec3d.Cplx in,  boolean swap) {
	Vec2d.Cplx in2d = Vec2d.createCplx( in.vectorWidth(), in.vectorHeight());
	in2d.project( in );
	Vec2d.Real pw   = Vec2d.createReal(in2d);
	Transforms.computePowerSpectrum( in2d, pw, swap);
	return pw;
    }


    /** Return the power spectrum (as a new vector) of 'in' */
    public static Vec2d.Real pwSpec(Vec2d.Cplx in ) {
	return pwSpec(in, true );
    }
    
    /** Return the power spectrum (as a new vector) of 'in' */
    public static Vec2d.Real pwSpec(Vec3d.Cplx in ) {
	return pwSpec(in, true );
    }


    
    /** Create an image by FFTin back 'in' to spatial */
    public static Vec2d.Real spatial(Vec2d.Cplx in ) {
	Vec2d.Real pw = Vec2d.createReal(in);
	Vec2d.Cplx tmp = in.duplicate();
	Transforms.fft2d( tmp, true );
	pw.copy( tmp );
	return pw;
    }
    
    /** Create an image by FFTin back 'in' to spatial */
    public static Vec2d.Real spatial(Vec3d.Cplx in ) {
	Vec2d.Real pw  = Vec2d.createReal(in.vectorWidth(), in.vectorHeight());
	Vec2d.Cplx tmp = Vec2d.createCplx( pw );
	tmp.project( in );
	Transforms.fft2d( tmp, true );
	pw.copy( tmp );
	return pw;
    }
    
    
    
    /** Create an image by FFTin back 'in' to spatial,
     *  clipping and scaling it afterwards */
    static Vec2d.Real spatial(Vec2d.Cplx in, SimParam.CLIPSCALE clipScale ) {
	Vec2d.Real pw = spatial(in);
	if (clipScale==SimParam.CLIPSCALE.CLIP)
	    clipAndScale( pw, true, false ); 
	if (clipScale==SimParam.CLIPSCALE.BOTH)
	    clipAndScale( pw, true, true ); 
	return pw;
    }


    /** Fades borders (sizes px) of the input to zero.
     *  Done by multiplying sin^2(x) mapped [0..px] to [pi/2..0].
     *  Good step before zero-padding data for high-res FFT. 
     *	@param img Vector to fade borders
     *	@param px  Size of faded region
     *  */
    public static void fadeBorderCos( Vec2d.Real img, int px ) {
	int w=img.vectorWidth(), h=img.vectorHeight();
	float [] dat = img.vectorData();
	final double fac = 1./px * Math.PI/2.;

	// top
	for (int y=0; y<px; y++)
	    for (int x=0; x<w ; x++)
		dat[x + y * w] *= Math.pow( Math.sin( y * fac ) , 2 );
	// bottom
	for (int y=h-px; y<h; y++)
	    for (int x=0; x<w ; x++)
		dat[x + y * w] *= Math.pow( Math.sin( (h-y-1) * fac ) , 2 );
	// left
	for (int y=0; y<h; y++)
	    for (int x=0; x<px ; x++)
		dat[x + y * w] *= Math.pow( Math.sin( x * fac ) , 2 );
	// right
	for (int y=0; y<h; y++)
	    for (int x=w-px; x<w ; x++)
		dat[x + y * w] *= Math.pow( Math.sin( (w-x-1) * fac ) , 2 );
	
	img.syncBuffer();
    }

    /** Fades borders (sizes px) of the input to zero.
     *  Done by multiplying sin^2(x) mapped [0..px] to [pi/2..0].
     *  Good step before zero-padding data for high-res FFT. 
     *	@param img Vector to fade borders
     *	@param px  Size of faded region
     *  */
    public static void fadeBorderCos( Vec2d.Cplx img, int px ) {
	int w=img.vectorWidth(), h=img.vectorHeight();
	float [] dat = img.vectorData();
	final double fac = 1./px * Math.PI/2.;

	// top
	for (int y=0; y<px; y++) {
	    for (int x=0; x<w ; x++) {
		double f = Math.pow( Math.sin( y * fac ) , 2 );
		dat[2*(x + y * w)+0] *= f;
		dat[2*(x + y * w)+1] *= f;
	    }
	}
	// bottom
	for (int y=h-px; y<h; y++) {
	    for (int x=0; x<w ; x++) {
		double f = Math.pow( Math.sin( (h-y-1) * fac ) , 2 );
		dat[2*(x + y * w)+0] *= f;
		dat[2*(x + y * w)+1] *= f;
	    }
	}
	// left
	for (int y=0; y<h; y++) {
	    for (int x=0; x<px ; x++) {
		double f = Math.pow( Math.sin( x * fac ) , 2 );
		dat[2*(x + y * w)+0] *= f;
		dat[2*(x + y * w)+1] *= f;
	    }
	}
	// right
	for (int y=0; y<h; y++) {
	    for (int x=w-px; x<w ; x++) {
		double f = Math.pow( Math.sin( (w-x-1) * fac ) , 2 );
		dat[2*(x + y * w)+0] *= f;
		dat[2*(x + y * w)+1] *= f;
	    }
	}
	
	img.syncBuffer();
    }



    /** Multiplies / overlays a region in the input image 'cntrl' with
     *  a sin pattern defined by kx, ky, pha. This function
     *  can be used to visualize the fitted pattern. 
     *	@param cntrl Vector to output to / modify
     *	@param kx x-coordinate of shift vector
     *	@param ky y-coordinate of shift vector
     *	@param pha Phase of pattern
     *  */
    public static void visualizePattern( Vec2d.Real cntrl, 
	    final double kx, final double ky, final double pha ) {
	    
	    if (cntrl.vectorWidth()!=cntrl.vectorHeight()) throw new
		RuntimeException("Please use a square image, w==h");
	    final int len = cntrl.vectorWidth();
	    float [] dat = cntrl.vectorData();

	    for (int y=len/2-50;y<len/2+150; y++)
	    for (int x=len/2-50;x<len/2+150; x++) {
		double v = Math.cos( 2*Math.PI* (-kx*x + ky*y)/(len) + pha );
		dat[ x + y*len ] = dat[ x + y*len]*(float)(v+2)*0.6f ;
	    }

	    cntrl.syncBuffer();
    }


    /** Subtracts a constant value (background) from each
     *  pixel.
     *  @param vec Image to work on
     *  @param bgr Value to subtract
     *  @return Percentage of clipped (lower-than-zero) values */
    public static double subtractBackground( Vec.Real vec,
	final double bgr ) {
	    
	int zCount=0;
	float [] dat = vec.vectorData();

	for (int i=0; i<dat.length; i++) {
	    dat[i]-=(float)bgr;
	    if (dat[i] <0) {
		dat[i]=0;
		zCount++;
	    }
	}
	
	vec.syncBuffer();

	return ((double)zCount) / dat.length;

    }

    /** Subtracts a constant value (background) from each
     *  pixel. Will set any imaginary component to zero!
     *  @param vec Image to work on
     *  @param bgr Value to subtract
     *  @return Percentage of clipped (lower-than-zero) values */
    public static double subtractBackground( Vec.Cplx vec,
	final double bgr ) {
	    
	int zCount=0;
	float [] dat = vec.vectorData();

	for (int i=0; i<vec.vectorSize(); i++) {
	    dat[2*i]-=(float)bgr;
	    dat[2*i+1]=0;
	    if (dat[2*i] <0) {
		dat[2*i]=0;
		zCount++;
	    }
	}
	
	vec.syncBuffer();

	return ((double)zCount) / dat.length;

    }


    /** Zero-clippes and scales an image.
     *  @param vec Image to work on
     *  @param clip Set values lower than 0 to 0
     *  @param scale Scale 0...max to 0..255 */
    public static void clipAndScale( Vec.Real vec, boolean clip, boolean scale ) {
	
	float [] dat = vec.vectorData();
	float max = 0;
	
	if ( scale ) {
	    // find maximum
	    for (int i=0; i<dat.length; i++) 
		if (dat[i]>max) max = dat[i];
	    // scale
	    for (int i=0; i<dat.length; i++) 
		dat[i] *= 255.f/max;
	}	

	if ( clip ) {
	    for (int i=0; i<dat.length; i++) 
		if (dat[i]<0) dat[i]=0;
	}

	vec.syncBuffer();
    }



} public class SimAlgorithm {

    /** Run the SIM parameter estimation.
     * Default to override any phase information set in SimParam.
     * @param param  The SIM parameter instance to work on
     * @param inFFT  The input images (in Fourier space)
     * @param fitBand On which band to perform the kx,ky fit
     * @param fitExclude How much (in fraction of OTF support) to exclude from fit
     * @param idf    ImageDisplayFactory for intermediate output (may be null)
     * @param visualFeedback Feedback Amount of visual feedback, 0..4
     * @param tEst   Runtime measurement (may be null) */
    public static void estimateParameters( final SimParam param, 
	Vec2d.Cplx [][] inFFT, final int fitBand, final double fitExclude,
	final ImageDisplay.Factory idf, 
	int visualFeedback, Tool.Timer tEst ) {
	estimateParameters(param, inFFT, fitBand, fitExclude, idf,
	    visualFeedback, tEst, false);
    }

    /** Run the SIM parameter estimation 
     * @param param  The SIM parameter instance to work on
     * @param inFFT  The input images (in Fourier space)
     * @param fitBand On which band to perform the kx,ky fit
     * @param fitExclude How much (in fraction of OTF support) to exclude from fit
     * @param idf    ImageDisplayFactory for intermediate output (may be null)
     * @param visualFeedback Feedback Amount of visual feedback, 0..4
     * @param tEst   Runtime measurement (may be null) 
     * @param keepPhases If true, phase information from SimParam will be used in band separation */
    public static void estimateParameters( final SimParam param, 
	Vec2d.Cplx [][] inFFT, final int fitBand, final double fitExclude,
	final ImageDisplay.Factory idf, 
	int visualFeedback, Tool.Timer tEst, boolean keepPhases ) {

	final int w = inFFT[0][0].vectorWidth(), h = inFFT[0][0].vectorHeight();
	final OtfProvider otfPr = param.otf();

	ImageDisplay pwSt=null,pwSt2=null, spSt=null, spSt2=null;
	if ((visualFeedback>0)&&(idf!=null)) {
	    pwSt  = idf.create(w,h, "Power Spectra" );
	    spSt  = idf.create(w,h, "Spatial images");
	    pwSt2 = idf.create(2*w,2*h, "Power Spectra" );
	    spSt2 = idf.create(2*w,2*h, "Spatial images");
	}

	// if no ImageFactory, do not generate feedback
	if (idf==null) {
	    visualFeedback = 0;
	}

	if (fitBand!=1 && fitBand!=2 ) throw new RuntimeException(
	    "Fitband neither 1 nor 2");

	if (tEst!=null) tEst.start();
    
	// The attenuation vector helps well to fade out the DC component,
	// which is uninteresting for the correlation anyway
	Vec2d.Real otfAtt = Vec2d.createReal( param );
	otfPr.writeAttenuationVector( otfAtt, .99, 0.15*otfPr.getCutoff(), 0, 0  ); 
	
	// loop through pattern directions
	for (int angIdx=0; angIdx<param.nrDir(); angIdx++) {
	
	    final SimParam.Dir dir = param.dir(angIdx);

	    // idx of low band (phase detection) and high band (shift vector detection)
	    // will be the same for two-beam
	    final int lb = 1;
	    final int hb = (param.dir(angIdx).nrBand()==3)?(3):(1);
	    final int fb = (fitBand==1)?(lb):(hb);

	    // compute band separation
	    Vec2d.Cplx [] separate = Vec2d.createArrayCplx( dir.nrComp(), w, h);
	    Vec2d.Cplx [] inputCpy = Vec2d.createArrayCplx( dir.nrPha(), w, h );
	    
	    for (int pha = 0 ; pha<dir.nrPha(); pha++) {
		inputCpy[pha] = inFFT[angIdx][pha].duplicate();
		inputCpy[pha].scal( new Cplx.Float( (float)dir.getIntensityQuotient(pha)) );
	    }
	    
	   
	    if (!keepPhases) {
		// reset any previously set phases
		dir.resetPhases();

		BandSeparation.separateBands( inputCpy , separate , 
		    0, dir.nrBand(), null);
	    } else {
		Tool.trace("Fitting with preset relative phases");
		BandSeparation.separateBands( inputCpy , separate , 
		    dir.getPhases(), dir.nrBand(), null);
	    }


	    // duplicate input, as we will modify it for coarse peak finding
	    // TODO: this could be avoided it neither coarse peak finding
	    // NOR visual feedback is requested
	    Vec2d.Cplx c0 = separate[0].duplicate();
	    Vec2d.Cplx c1 = separate[lb].duplicate();
	    Vec2d.Cplx c2 = separate[hb].duplicate();

	    // dampen region around DC 
	    c0.times( otfAtt );
	    c1.times( otfAtt );
	    c2.times( otfAtt ); 
	    
	    // compute correlation: ifft, mult. in spatial, fft back
	    Transforms.fft2d( c0, true);
	    Transforms.fft2d( c1, true);
	    Transforms.fft2d( c2, true);
	    c1.timesConj( c0 );
	    c2.timesConj( c0 );
	    Transforms.fft2d( c1, false);
	    Transforms.fft2d( c2, false);
	    
	    double [] peak ;
	    double minDist =2;
	   
	    if (fitExclude>0) {
		// find the highest peak in corr of band0 to highest band 
		// with minDist of otfCutoff from origin, store in 'param'
		minDist = fitExclude * otfPr.getCutoff() / param.pxlSizeCyclesMicron();
		peak = Correlation.locatePeak(  (fitBand==1)?(c1):(c2) , minDist );
		
		Tool.trace(String.format("Peak: (dir %1d) located (min %4.0f) at x %5.0f y %5.0f",
		    angIdx, minDist, peak[0], peak[1]));
	    } else {
		// just go with the values stored in param for the coarse peak estimate
		peak = new double[] { param.dir(angIdx).px(fitBand), param.dir(angIdx).py(fitBand) };
		Tool.trace(String.format("Peak: (dir %1d) from sim param preset at x %5.0f y %5.0f",
		    angIdx, peak[0], peak[1]));
	    }

	    // fit the peak to sub-pixel precision by cross-correlation of
	    // Fourier-shifted components
	    Vec2d.Real cntrl    = Vec2d.createReal(30,10);
	    peak = Correlation.fitPeak( separate[0], separate[fb], 0, fitBand, otfPr,
		-peak[0], -peak[1], 0.05, 2.5, cntrl );

	    // Now, either three beam / 3 bands ...
	    if (lb!=hb) {
		
		// peak should contain the shift band0<->band2, so if band0<->band1
		// was fitted, multiply by 2
		if (fitBand==1) {
		    peak[0]*=2; peak[1]*=2;
		}

		// At the peak position found, extract phase and modulation from band0 <-> band 1
		Cplx.Double p1 = Correlation.getPeak( separate[0], separate[lb], 
		    0, 1, otfPr, peak[0]/2, peak[1]/2, 0.05 );

		// TODO: this is a quick fix, this should be done properly by looking at the OTF
		// overlap, to figure out if parameter extraction is possible. This just catches
		// the 'definitely not possible' case
		Cplx.Double p2;

		if ( Math.hypot( peak[0], peak[1] ) * param.pxlSizeCyclesMicron() > otfPr.getCutoff() ) {
		    Tool.trace( "Second peak outside of OTF overlap (non-linear SIM?)" );
		    Tool.trace( "-> copying modulation depth for now...");
		    p2 = p1;   
		} else {
		    // Extract modulation from band0 <-> band 2
		    p2 = Correlation.getPeak( separate[0], separate[hb], 
			0, 2, otfPr, peak[0], peak[1], 0.05 );
		}

		Tool.trace(
		    String.format("Peak: (dir %1d): fitted --> x %7.3f y %7.3f p %7.3f (m %7.3f, %7.3f)", 
		    angIdx, peak[0], peak[1], p1.phase(), p1.hypot(), p2.hypot() ));
	
		// store the result
		param.dir(angIdx).setPxPy(   -peak[0], -peak[1] );
		param.dir(angIdx).setPhaOff( p1.phase() );
		param.dir(angIdx).setModulation( 1, p1.hypot() );
		param.dir(angIdx).setModulation( 2, p2.hypot() );
	    }
	    
	    // ... or two-beam / 2 bands
	    if (lb==hb) {
		// get everything from one correlation band0 to band1
		Cplx.Double p1 = Correlation.getPeak( separate[0], separate[1], 
		    0, 1, otfPr, peak[0], peak[1], 0.05 );

		Tool.trace(
		    String.format("Peak: (dir %1d): fitted --> x %7.3f y %7.3f p %7.3f (m %7.3f)", 
		    angIdx, peak[0], peak[1], p1.phase(), p1.hypot() ));
	
		// store the result
		param.dir(angIdx).setPxPy(   -peak[0], -peak[1] );
		param.dir(angIdx).setPhaOff( p1.phase() );
		param.dir(angIdx).setModulation( 1, p1.hypot() );
	    }



	    // --- output visual feedback of peak fit ---
	    if (visualFeedback>0) {
		
		// mark region excluded from peak finder
		// output the peaks found, with circles marking them, and the fit result in
		// the top corner for the correlation band0<->band2
		ImageDisplay.Marker excludedDC = 
		    new ImageDisplay.Marker(w/2,h/2,minDist*2,minDist*2,true);
		
		Vec2d.Real fittedPeak = SimUtils.pwSpec( (fitBand==1)?(c1):(c2) );
		fittedPeak.paste( cntrl, 0, 0, false );
		
		double f=(lb!=hb)?(fitBand/2.):(1);

		pwSt.addImage( fittedPeak, "dir "+angIdx+" c-corr band 0<>band "+fitBand,
		    new ImageDisplay.Marker( w/2-peak[0]*f, h/2+peak[1]*f, 10, 10, true),
		    excludedDC);
		
		// if there is a low band, also add it
		if ((visualFeedback>1)&&(lb!=hb))
		    pwSt.addImage( SimUtils.pwSpec( c1 ), "dir "+angIdx+" c-corr band 0<>low",
			new ImageDisplay.Marker( w/2-peak[0]/2, h/2+peak[1]/2, 10, 10, true));
	    }
		

	    // --- output visual feedback of overlapping regions (for all bands) ---
	    if (visualFeedback>1)  
	    for (int b=1; b<param.nrBand(); b++) {	
	    
		SimParam.Dir par = param.dir(angIdx);

		// find common regions in low and high band
		Vec2d.Cplx b0 = separate[0  ].duplicate();
		Vec2d.Cplx b1 = separate[2*b].duplicate();
	    
		Correlation.commonRegion( b0, b1, 0, b, otfPr,  
		    par.px(b), par.py(b), 0.15, (b==1)?(.2):(.05), true);

		// move the high band to its correct position
		b1.fft2d( true );
		b1.fourierShift( par.px(b), -par.py(b));
		b1.fft2d( false );
	
		// apply phase correction
		b1.scal( Cplx.Float.fromPhase( param.dir(angIdx).getPhaOff()*b ));
       
		// output the full shifted bands
		if ( visualFeedback>2 )  {
		    // only add band0 once	
		    if ( b==1 ) {
			Vec2d.Cplx btmp = separate[0].duplicate();
			otfPr.maskOtf( btmp, 0, 0);
			pwSt.addImage(SimUtils.pwSpec( btmp ), String.format(
			    "a%1d: full band0", angIdx, b ));
		    }

		    // add band1, band2, ...
		    Vec2d.Cplx btmp = separate[2*b].duplicate();
		    btmp.fft2d( true );
		    btmp.fourierShift( par.px(b), -par.py(b) );
		    btmp.fft2d( false );
		    otfPr.maskOtf( btmp, par.px(b), par.py(b));

		    pwSt.addImage(SimUtils.pwSpec( btmp ), String.format( 
			"a%1d: full band%1d (shifted %7.3f %7.3f)",
			angIdx,  b, par.px(b), par.py(b))); 
		}

		// output power spectra of common region
		pwSt.addImage(SimUtils.pwSpec( b0 ), String.format(
		    "a%1d: common region b0<>b%1d, band0", angIdx, b )); 
		pwSt.addImage(SimUtils.pwSpec( b1 ), String.format( 
		    "a%1d: common region b0<>b%1d, band%1d",angIdx, b, b)); 

		// output spatial representation of common region
		spSt.addImage(SimUtils.spatial( b0 ), String.format(
		    "a%1d: common region b0<>b%1d, band0", angIdx, b )); 
		spSt.addImage(SimUtils.spatial( b1 ), String.format( 
		    "a%1d: common region b0<>b%1d, band%1d",angIdx, b, b)); 
	    }
	
	}
    
	if (tEst!=null) tEst.stop();

	if (idf!=null) {
	    pwSt.display();
	    spSt.display();
	    pwSt2.display();
	    spSt2.display();
	}


    }


    /** Estimate absolute phases for each raw data frame.
     *  This follows the non-iterative phase estimation based on auto-correlation by Wicker et al.
     *  It requires(!) the k-Vectors to be set correctly (usually doable by running the parameter
     *  estimation with roughly correct relative phases)
     * @param param  The SIM parameter instance to work on
     * @param inFFT  The input images (in Fourier space)
     * @param fitExclude How much (in fraction of OTF support) to exclude from fit
     * @param tEst   Runtime measurement (may be null) 
    */
    public static void estimateAbsolutePhases( final SimParam param, 
	Vec2d.Cplx [][] inFFT, 
	Tool.Timer tEst ) {

	Tool.tell("running individual phase estimator");

	if (tEst != null) tEst.start();
   
	final OtfProvider otfPr = param.otf();

	// loop all directions
	for (int angIdx = 0; angIdx < param.nrDir(); angIdx ++ ) {
	    final SimParam.Dir par = param.dir(angIdx);
	    
	    // run Kai Wicker auto-correlation
	    double [] pha = new double[par.nrPha()];
	    for (int i=0; i < par.nrPha() ; i++) {
	       
		// copy input, weight with otf
		Vec2d.Cplx ac =  inFFT[angIdx][i].duplicate();
		otfPr.applyOtf( ac, 0); 
		
		// compute auto-correlation at px,py (shift of band1)
		Cplx.Double corr = Correlation.autoCorrelation( 
		    ac, par.px(1), par.py(1) );

		Tool.trace(String.format("a%1d img %1d, Phase(Wicker et. al.) : %5.3f  ",
		    angIdx, i, corr.phase()));

		pha[i] = corr.phase();
	    }
	    par.setPhases( pha, true );	
	    
	}

	if (tEst != null) tEst.stop();
    }


    // TODO: deprecate this function at some point
    // it is only here to be compatible with existing scripts / examples,
    // use more complete version below

    /** Run the SIM reconstruction 
     * @param param  The SIM parameter instance to work on
     * @param inFFT  The input images (in Fourier space)
     * @param idf    ImageDisplayFactory for intermediate output (may be null)
     * @param visualFeedback Feedback Amount of visual feedback, -1 (off) ..4
     * @param otfBeforeShift Apply the OTF before shifting bands
     * @param imgClipScale Clip zero values and scale (0..255) output images?
     * @param tRec   Runtime measurement (may be null) 
     * @return The reconstructed image */
    public static Vec2d.Real runReconstruction( final SimParam param, 
	Vec2d.Cplx [][] inFFT, ImageDisplay.Factory idf, int visualFeedback, 
	final boolean otfBeforeShift, final SimParam.CLIPSCALE imgClipScale, 
	Tool.Timer tRec ) {
	return runReconstruction( param, inFFT, idf, visualFeedback, otfBeforeShift, imgClipScale, null, null, tRec);
    }




    /** Run the SIM reconstruction 
     * @param param  The SIM parameter instance to work on
     * @param inFFT  The input images (in Fourier space)
     * @param idf    ImageDisplayFactory for intermediate output (may be null)
     * @param visualFeedback Feedback Amount of visual feedback, -1 (off) ..4
     * @param otfBeforeShift Apply the OTF before shifting bands
     * @param imgClipScale Clip zero values and scale (0..255) output images?
     * @param tRec   Runtime measurement (may be null) 
     * @param widefieldResult The widefield image computed during the reconstruction (may be null)
     * @param filteredWidefieldResult The filtered widefield image (may be null)
     * @return The reconstructed image */ 
    public static Vec2d.Real runReconstruction( final SimParam param, 
	Vec2d.Cplx [][] inFFT, ImageDisplay.Factory idf, int visualFeedback, 
	final boolean otfBeforeShift, final SimParam.CLIPSCALE imgClipScale, 
	Vec2d.Real widefieldResult,
	Vec2d.Real filteredWidefieldResult,
	Tool.Timer tRec ) {

	if (tRec != null) tRec.start();	

	final int w = inFFT[0][0].vectorWidth(), h = inFFT[0][0].vectorHeight();
    
	final double apoB=param.getApoBend(), apoF=param.getApoCutoff(); // Bend and mag. factor of APO
	
	final OtfProvider otfPr = param.otf();

	ImageDisplay pwSt=null,pwSt2=null, spSt=null, spSt2=null;
	
	if (idf!=null) {
	    pwSt  = idf.create(w,h, "Power Spectra" );
	    spSt  = idf.create(w,h, "Spatial images");
	    pwSt2 = idf.create(2*w,2*h, "Power Spectra" );
	    spSt2 = idf.create(2*w,2*h, "Spatial images");
	} else {
	    if (visualFeedback>=0) {
		Tool.trace("Reconstruction: No image display, turning off "+
		"all intermediate output");
		visualFeedback = -1;
	    }
	}
	
	// setup WienerFilter
	Tool.tell("Setting up Wiener filter");
	WienerFilter wFilter = new WienerFilter( param );
	double wienParam     = param.getWienerFilter();

	if (visualFeedback>0) {
	    Vec2d.Real wd = wFilter.getDenominator(wienParam);
	    wd.reciproc();
	    wd.normalize();
	    Transforms.swapQuadrant( wd );
	    pwSt2.addImage(wd, "Wiener denominator");
	}
	
	// vectors to store the result
	Vec2d.Cplx fullResult    = Vec2d.createCplx( param, 2);
	
	// zero-band OTF (for RL filtering)
	Vec2d.Cplx inputOtf = null;
	if ( param.useRLonInput() ) {
	    inputOtf = Vec2d.createCplx(param);
	    otfPr.writeOtfVector( inputOtf, 0,0,0);
	    if ( visualFeedback>1 )  {
		Vec2d.Real tmp  = Vec2d.createReal(param);
		tmp.copy( inputOtf );
		pwSt.addImage( tmp, "Widefield OTF" );
	    }
	}
    
	// loop all pattern directions
	for (int angIdx = 0; angIdx < param.nrDir(); angIdx ++ ) 
	{
	    final SimParam.Dir par = param.dir(angIdx);
	    Tool.tell("Reconstr. for angle "+(angIdx+1)+"/"+param.nrDir());
	    
	    Vec2d.Cplx [] separate  = Vec2d.createArrayCplx( par.nrComp(), w, h);
	  
	    // copy into temp. array (to not override input data) 
	    // and apply correction factor
	    Vec2d.Cplx [] tmpArray = Vec2d.createArrayCplx( par.nrPha(), w, h);
	    for (int i=0; i< par.nrPha() ;i++) {
		tmpArray[i].copy( inFFT[angIdx][i] );
		tmpArray[i].scal( new Cplx.Float( (float)par.getIntensityQuotient(i) ) );
		Tool.trace(String.format(" Input data intensity corrected: a%1d p%1d --> %7.5f",
		    angIdx, i, par.getIntensityQuotient(i)));
	    }

	    

	    // ---- Richardson-Lucy: Deconvolve input data here ----
	    if ( param.useRLonInput() ) {
		
			// deconvolve the input data
		for (int i=0; i< par.nrPha() ;i++) { 

		    if (visualFeedback>1) {
			spSt.addImage( SimUtils.spatial( tmpArray[i]), 
			    "input before deconv., ang "+angIdx+", phase "+i);
		    }

		    RLDeconvolution.deconvolve( tmpArray[i], inputOtf, 
			param.getRLiterations(), true);
		    
		    if (visualFeedback>0) {
			spSt.addImage( SimUtils.spatial( tmpArray[i]),
			    "Deconvolved input, ang "+angIdx+", phase "+i);
		    }
		
		}
	    }

	    // use the temp array as input for the band separation
	    BandSeparation.separateBands( tmpArray , separate , 
		par.getPhases(), par.nrBand(), par.getModulations());
	    

	    // Wiener filter: Apply OTF here
	    if (otfBeforeShift && param.useWienerFilter() )
		for (int i=0; i<(par.nrComp()) ;i++)  
		    otfPr.applyOtf( separate[i], (i+1)/2);

	    // ------- Shifts to correct position ----------
	    Vec2d.Cplx [] shifted		= Vec2d.createArrayCplx(5, 2*w, 2*h);

	    // band 0 is DC, so does not need shifting, only a bigger vector
	    SimUtils.placeFreq( separate[0],  shifted[0]);
	    
	    // higher bands need shifting
	    for ( int b=1; b<par.nrBand(); b++) {
		
		Tool.trace("reconstr.: dir "+angIdx+": shift band: "+b+" to: "+par.px(b)+" "+par.py(b));
		
		// first, copy to larger vectors
		int pos = b*2, neg = (b*2)-1;	// pos/neg contr. to band
		SimUtils.placeFreq( separate[pos] , shifted[pos]);
		SimUtils.placeFreq( separate[neg] , shifted[neg]);

		// then, fourier shift
		SimUtils.fourierShift( shifted[pos] ,  par.px(b),  par.py(b) );
		SimUtils.fourierShift( shifted[neg] , -par.px(b), -par.py(b) );
	    }
	   
	    // ------ OTF multiplication or masking ------
	   
	    if ( param.useWienerFilter() ) {
		if (!otfBeforeShift) {
		    // multiply with shifted OTF
		    otfPr.applyOtf( shifted[0], 0 );
		    for (int b=1; b<par.nrBand(); b++) {
			int pos = b*2, neg = (b*2)-1;	// pos/neg contr. to band
			otfPr.applyOtf( shifted[pos], b,  par.px(b),  par.py(b) );
			otfPr.applyOtf( shifted[neg], b, -par.px(b), -par.py(b) );
		    }
		} else {
		    // or mask for OTF support
		    for (int b=1; b<par.nrBand(); b++) {
			int pos = b*2, neg = (b*2)-1;	// pos/neg contr. to band
			otfPr.maskOtf( shifted[pos],  par.px(b),  par.py(b) );
			otfPr.maskOtf( shifted[neg], -par.px(b), -par.py(b) );
		    }
		}
	    }
	    // ------ Sum up result ------
	    
	    if (param.useRLonOutput()) {
		    shifted[0].scal( 1.f / param.nrDir());
	    }
	    
	    for (int i=0;i<par.nrBand()*2-1;i++) { 
		    fullResult.add( shifted[i] ); 
	    }
	
	    
	    // ------ Output intermediate results ------
	    
	    if (visualFeedback>0) {
	
		Tool.tell("Computing interm. results");

		// per-direction results
		Vec2d.Cplx result = Vec2d.createCplx(2*w,2*h);
		for (int i=0;i<par.nrBand()*2-1;i++)  
		    result.add( shifted[i] ); 

		// loop bands in this direction
		for (int i=0;i<par.nrBand();i++) {     

		    // TODO: All these should also use RL-filtering if set!
		    // get wiener denominator for (direction, band), add to full denom for this band
		    Vec2d.Real denom = wFilter.getIntermediateDenominator( angIdx, i, wienParam);
		
		    // add up +- shift for this band
		    Vec2d.Cplx thisband   = shifted[i*2];
		    if (i!=0)
			thisband.add( shifted[i*2-1] );
	
		    // output the wiener denominator
		    if (visualFeedback>1) {
			Vec2d.Real wd = denom.duplicate();
			wd.reciproc();
			wd.normalize();
			Transforms.swapQuadrant( wd );
			pwSt2.addImage( wd, String.format(
			    "a%1d: OTF/Wiener band %1d",angIdx,(i/2) ));
		    }
		    
		    // apply filter and output result
		    thisband.times( denom );
		    
		    pwSt2.addImage( SimUtils.pwSpec( thisband ) ,String.format(
			"a%1d: band %1d",angIdx,i));
		    spSt2.addImage( SimUtils.spatial( thisband, imgClipScale ) ,String.format(
			"a%1d: band %1d",angIdx,i));
		}

		// per direction wiener denominator	
		Vec2d.Real fDenom =  wFilter.getIntermediateDenominator( angIdx, wienParam);	
		result.times( fDenom );
		    
		// output the wiener denominator
		if (visualFeedback>1) {
		    Vec2d.Real wd = fDenom.duplicate();
		    wd.reciproc();
		    wd.normalize();
		    Transforms.swapQuadrant( wd );
		    pwSt2.addImage( wd, String.format(
			"a%1d: OTF/Wiener all bands",angIdx ));
		}
		
		pwSt2.addImage( SimUtils.pwSpec( result ) ,String.format(
		    "a%1d: all bands",angIdx));
		spSt2.addImage( SimUtils.spatial( result, imgClipScale ) ,String.format(
		    "a%1d: all bands",angIdx));
	    
		// power spectra before shift
		if (visualFeedback>2) { 
		    for (int i=0; i<(par.nrBand()*2-1) ;i++)  
		    pwSt.addImage( SimUtils.pwSpec( separate[i] ), String.format(
			"a%1d, sep%1d, seperated band", angIdx, i));
		}
	   
	    }


	}   
	
	// -- done loop all pattern directions, 'fullResult' now holds the image --

	Vec2d.Real fullResultImage = null;

	// ------------------------------------------------------------------------
	// Wiener filter on outout
	// ------------------------------------------------------------------------

	if ( param.useWienerFilter() ) {

	    Tool.tell("Applying Wiener filter");

	    // multiply by wiener denominator
	    Vec2d.Real denom = wFilter.getDenominator( wienParam );
	    fullResult.times(denom);
	    
	    if (visualFeedback>0) {
		pwSt2.addImage(  SimUtils.pwSpec( fullResult), "full (w/o APO)");
		spSt2.addImage(  SimUtils.spatial(fullResult, imgClipScale), "full (w/o APO)");
	    }

	    // apply apotization filter
	    Vec2d.Cplx apo = Vec2d.createCplx(2*w,2*h);
	    otfPr.writeApoVector( apo, apoB, apoF);
	    
	    if (visualFeedback>1) {
		Vec2d.Real tmp = Vec2d.createReal(2*w,2*h);
		tmp.copy( apo );
		Transforms.swapQuadrant( tmp );
		pwSt2.addImage( tmp, "Apo filter");
	    }
	    
	    fullResult.times(apo);
	    
	    fullResultImage = SimUtils.spatial( fullResult, imgClipScale);

	    if (spSt2 != null) 
		spSt2.addImage( fullResultImage, "full result");


	    if (visualFeedback>0) {
		pwSt2.addImage( SimUtils.pwSpec( fullResult), "full result");
	    }


	    // Add wide-field for comparison
	    if (visualFeedback>=0 || widefieldResult != null || filteredWidefieldResult != null ) {
	    
		Tool.tell("Computing wide-field");
		
		// obtain the low freq result
		Vec2d.Cplx lowFreqResult = Vec2d.createCplx( param, 2);
		
		// have to do the separation again, result before had the OTF multiplied
		for (int angIdx = 0; angIdx < param.nrDir(); angIdx ++ ) {
		    
		    final SimParam.Dir par = param.dir(angIdx);
		    
		    Vec2d.Cplx [] separate  = Vec2d.createArrayCplx( par.nrComp(), w, h);
		    BandSeparation.separateBands( inFFT[angIdx] , separate , 
			par.getPhases(), par.nrBand(), par.getModulations());

		    Vec2d.Cplx tmp  = Vec2d.createCplx( param, 2 );
		    SimUtils.placeFreq( separate[0],  tmp);
		    lowFreqResult.add( tmp );
		}	
		
		// now, output the widefield
		if (visualFeedback>0)
		    pwSt2.addImage( SimUtils.pwSpec(lowFreqResult), "Widefield" );
		if (visualFeedback>=0)
		    spSt2.addImage( SimUtils.spatial(lowFreqResult, imgClipScale), "Widefield" );
		
		if (widefieldResult!=null) { 
		    widefieldResult.copy( SimUtils.spatial(lowFreqResult, imgClipScale));
		    Tool.trace("generating widefield output");
		}
	    
		// otf-multiply and wiener-filter the wide-field
		otfPr.otfToVector( lowFreqResult, 0, 0, 0, false, false ); 

		Vec2d.Real lfDenom = wFilter.getWidefieldDenominator( wienParam );
		lowFreqResult.times( lfDenom );
	       
		// mask out freq. that could not have passed 
		// (does not really change the image)
		otfPr.maskOtf( lowFreqResult, 0,0 );
		
		// could apodize the result, but that would just re-apply 
		// the OTF (more or less)
		// Vec2d.Cplx apoLowFreq = Vec2d.createCplx(2*w,2*h);
		// otfPr.writeApoVector( apoLowFreq, 0.4, 1.2);
		// lowFreqResult.times(apoLowFreq);
		
		if (visualFeedback>0)
		    pwSt2.addImage( SimUtils.pwSpec( lowFreqResult), "filtered Widefield" );
		if (visualFeedback>=0)
		    spSt2.addImage( SimUtils.spatial( lowFreqResult, imgClipScale), "filtered Widefield" );
		
		if (filteredWidefieldResult!=null) {
		    filteredWidefieldResult.copy( SimUtils.spatial(lowFreqResult,imgClipScale));
		    Tool.trace("generating filtered widefield output");
		}

	    }
	}	

	// ------------------------------------------------------------------------
	// Richardson-Lucy deconvolution on output
	// ------------------------------------------------------------------------

	if ( param.useRLonOutput() ) {

	    Tool.tell("Applying RL deconvolution filter");

	    
	    // generate the effective SIM OTF
	    Vec2d.Cplx otfSim = Vec2d.createCplx( param, 2);
	    
	    Vec2d.Cplx otfTmpPos = Vec2d.createCplx( param, 2);
	    Vec2d.Cplx otfTmpNeg = Vec2d.createCplx( param, 2);

	    for (int ang = 0; ang<param.nrDir(); ang++) {
	    
		final SimParam.Dir par = param.dir(ang);
		
		for (int b=0; b<param.nrBand(); b++) {
		    
		    otfPr.writeOtfVector( otfTmpPos, b,  par.px(b),  par.py(b));
		    otfPr.writeOtfVector( otfTmpNeg, b, -par.px(b), -par.py(b));
		    
		    if (b==0) {
			otfTmpPos.scal( 1.f/param.nrDir());
			otfSim.add( otfTmpPos );
		    } else {
			otfSim.add( otfTmpPos );
			otfSim.add( otfTmpNeg );
		    }

		}
	    } 

	    
	    // visualize the SIM OTF
	    if (visualFeedback>1) {
		Vec2d.Real otfVis = Vec2d.createReal( param,2);
		otfVis.copy( otfSim );
		Transforms.swapQuadrant( otfVis );
		pwSt2.addImage( otfVis, "SIM OTF added up");
	    }
	    
	    // output result before RL-filtering
	    if (visualFeedback>0) {
		pwSt2.addImage(  SimUtils.pwSpec( fullResult), "full unfiltered");
		spSt2.addImage(  SimUtils.spatial(fullResult), "full unfiltered");
	    }

	    // deconvolve the result
	    RLDeconvolution.deconvolve( fullResult, otfSim, 
		param.getRLiterations(), true);
	    
	    fullResultImage = SimUtils.spatial( fullResult, imgClipScale);

	    if (spSt2 != null) 
		spSt2.addImage( fullResultImage, "full result (RL)");


	    if (visualFeedback>0) {
		pwSt2.addImage( SimUtils.pwSpec( fullResult), "full result (RL)");
	    }


	    // Add wide-field for comparison
	    if (visualFeedback>=0 || widefieldResult != null || filteredWidefieldResult != null ) {
	    
		Tool.tell("Computing wide-field");
		
		// obtain the low freq result
		Vec2d.Cplx lowFreqResult = Vec2d.createCplx( param, 2);
		
		// have to do the separation again, result before had the OTF multiplied
		for (int angIdx = 0; angIdx < param.nrDir(); angIdx ++ ) {
		    
		    final SimParam.Dir par = param.dir(angIdx);
		    
		    Vec2d.Cplx [] separate  = Vec2d.createArrayCplx( par.nrComp(), w, h);
		    BandSeparation.separateBands( inFFT[angIdx] , separate , 
			par.getPhases(), par.nrBand(), par.getModulations());

		    Vec2d.Cplx tmp  = Vec2d.createCplx( param, 2 );
		    SimUtils.placeFreq( separate[0],  tmp);
		    lowFreqResult.add( tmp );
		}	
		
		// now, output the widefield
		if (visualFeedback>0)
		    pwSt2.addImage( SimUtils.pwSpec(lowFreqResult), "Widefield" );
		if (visualFeedback>=0)
		    spSt2.addImage( SimUtils.spatial(lowFreqResult, imgClipScale), "Widefield" );
	    
		if (widefieldResult!=null) { 
		    widefieldResult.copy( SimUtils.spatial(lowFreqResult, imgClipScale));
		    Tool.trace("generating widefield output");
		}

		// deconvolve the wide-field 
		Vec2d.Cplx zeroOrderOtf = Vec2d.createCplx(param,2);
		otfPr.writeOtfVector( zeroOrderOtf, 0,0,0);
		RLDeconvolution.deconvolve( lowFreqResult, zeroOrderOtf, 
		    param.getRLiterations(), true);
		
		if (visualFeedback>0)
		    pwSt2.addImage( SimUtils.pwSpec( lowFreqResult), "filtered Widefield" );
		if (visualFeedback>=0)
		    spSt2.addImage( SimUtils.spatial( lowFreqResult, imgClipScale), "filtered Widefield" );
		
		if (filteredWidefieldResult!=null) {
		    filteredWidefieldResult.copy( SimUtils.spatial(lowFreqResult,imgClipScale));
		    Tool.trace("generating filtered widefield output");
		}


	    }
	}	

	// -----------------------------------------------------------------------

	// stop timers
	if (tRec!=null) tRec.stop();	

	Tool.tell("done.");

	// output parameters
	Tool.trace(" ---- Reconstruction ---- ");
	Tool.trace( "\n"+param.prettyPrint(true));
	
	if (visualFeedback>=0) {
	    pwSt.display();
	    spSt.display();
	    pwSt2.display();
	}
	
	if (spSt2 != null)
	    spSt2.display();

	return fullResultImage;
    }


} public final class Vec2d {

    /** Vec2d is static only */
    private Vec2d() { 
	throw new AssertionError();	
    };
    
    /** Return a 2-dimensional vector sized w x h, real-valued */
    public static Real createReal( int w, int h) {
	return Vec.vf.createReal2D(w,h);
    }
 
    /** Return a 2-dimensional vector sized w x h, complex-valued  */
    public static Cplx createCplx( int w, int h) {
	return Vec.vf.createCplx2D(w,h);
    }
    
    /** Return a 2-dimensional vector, width and height of 'size'
     *  scaled by 'scale' */
    public static Real createReal( Size size, int scale ) {
	return createReal( scale*size.vectorWidth(), scale*size.vectorHeight() );
    }

    /** Return a 2-dimensional complex-valued vector, width and height of 'size'
     *  scaled by 'scale' */
    public static Cplx createCplx( Size size, int scale ) {
	return createCplx( scale*size.vectorWidth(), scale*size.vectorHeight() );
    }
    
    /** Return a 2-dimensional vector, width and height of 'size' */
    public static Real createReal( Size size ) {
	return createReal(size, 1);
    }

    /** Return a 2-dimensional complex-valued vector, width and height of 'size' */
    public static Cplx createCplx( Size size ) {
	return createCplx(size, 1);
    }
 
 
    /** Return a n-element array of 2-dimensional vector,
     *  sized w x h with real-valued elements */
    public static Real [] createArrayReal(int n, int w, int h) {
	Real [] ret = new Real[n];
	for (int i=0; i<n; i++)
	    ret[i] = createReal(w,h);
	return ret;
    }
    
    /** Return a n-element array of 2-dimensional vector,
     *  sized w x h with complex-valued elements */
    public static Cplx [] createArrayCplx(int n, int w, int h) {
	Cplx [] ret = new Cplx[n];
	for (int i=0; i<n; i++)
	    ret[i] = createCplx(w,h);
	return ret;
    }

    /** Throws a runtime exception if not all vectors are 
     *	not of same width and height. */
    public static void failSize( Size ... v ) {
	if (v.length<2) return;
	int w = v[0].vectorWidth(), h=v[0].vectorHeight();
	for (int i=1;i<v.length;i++)
	    if ((v[i].vectorWidth() != w) || (v[i].vectorHeight() != h))
		throw new RuntimeException("Vector len mismatch"
		    +i+" "+w+" "+h+" "+v[i].vectorWidth()+" "+v[i].vectorHeight());
    }
    
    /** Throws a runtime exception if not all vectors are 
     *  of the same, square size (N = width = height). Returns the size N.*/
    public static int checkSquare( Size ... v ) {
	if (v[0].vectorWidth()!=v[0].vectorHeight())
	    throw new RuntimeException("Vector len mismatch "
	        +v[0].vectorWidth()+" "+v[0].vectorHeight());
	failSize(v);
	return v[0].vectorWidth();
    }

    /** Object with width and height */
    public interface Size {
	/** Width */
	public int vectorWidth();
	/** Height */
	public int vectorHeight();
    }

    // =======================================================
    
    /**  Real-valued, 2-dimensional vector, float base-type */
    public interface Real extends Vec.Real, Size { 

	/** Get element x,y */
	public float get(int x, int y);
	/** Set element x,y */
	public void set(int x, int y, float v);
	
	/** Return a duplicate / clone of the input vector */
	public Real duplicate() ;
	
	/** Copy the vector x into a region of this vector.
	 *  @param ulx x-coordinate of upper left corner where to paste  
	 *  @param uly y-coordinate of upper left corner where to paste  
	 *  @param zero Zero the output vector before paste */
	public void paste( Vec2d.Real xIn, int ulx, int uly, boolean zero );

	/** width of the vector */
	public int vectorWidth();

	/** height of the vector */
	public int vectorHeight();
    
	/** project (sum along z, from start to end) a Vec3d into this vector */
	public void project( Vec3d.Real in, int start, int end );
	
	/** project (sum along z) a Vec3d into this vector */
	public void project( Vec3d.Real in );
	
	/** get slice nr n */
	public void slice( Vec3d.Real in, int n );
    
	/** Copy an input image into this vector.
	 *  The input image is assumed as 16bit unsigned(!) short values,
	 *  arranged short[ x + width * y ].
	 *  @param in The array to copy
	 * */
	public void setFrom16bitPixels( short [] in );
    
    }
    
    // =======================================================
    
    /** Complex-valued, 2-dimensional vector, float base-type */
    public interface Cplx extends Vec.Cplx, Size { 

	
	/** Get element x,y */
	public org.fairsim.linalg.Cplx.Float get(int x, int y);
	/** Set element x,y */
	public void set(int x, int y, org.fairsim.linalg.Cplx.Float v);
	/** Set element x,y */
	public void set(int x, int y, org.fairsim.linalg.Cplx.Double v);
	
	/** Return a duplicate / clone of the input vector */
	public abstract Cplx duplicate() ; 

	/** Compute the in-place fft */
	public void fft2d( boolean inverse );

	/** Copy the vector x into a region of this vector.
	 *  Short-hand for paste( xIn, yOut, 0, 0, xWidth(), xHeight(), ulx, uly, zero )
	 *  @param ulx x-coordinate of upper left corner where to paste  
	 *  @param uly y-coordinate of upper left corner where to paste  
	 *  @param zero Zero the output vector before paste */
	public void paste( Vec2d.Cplx xIn, int ulx, int uly, boolean zero );

	/** Copy a freq-space vector into this one.
	 * @param xIn The vector to copy into this one 
	 * @param xOff The offset in x-direction 
	 * @param yOff The offset in y-direction */
	public void pasteFreq( Vec2d.Cplx xIn , int xOff, int yOff );

	/** pasteFreq(xIn, 0, 0) */
	public void pasteFreq( Vec2d.Cplx xIn );



	/** Multiply a vector with Fourier shift theorem phases.
	 *  @param kx x-coordinate of shift
	 *  @param ky y-coordinate of shift */
	public void fourierShift( double kx, double ky);


	/** width of the vector */
	public int vectorWidth();

	/** height of the vector */
	public int vectorHeight();

	/** project (sum along z, from start to end) a Vec3d into this vector */
	public void project( Vec3d.Cplx in, int start, int end );
	
	/** project (sum along z) a Vec3d into this vector */
	public void project( Vec3d.Cplx in );
    
	/** get slice nr n */
	public void slice( Vec3d.Cplx in, int n );

	/** Copy an input image into this vector.
	 *  The input image is assumed as 16bit unsigned(!) short values,
	 *  arranged short[ x + width * y ]. The imaginary components
	 *  of the vector will be set to zero.
	 *  @param in The array to copy
	 * */
	public void setFrom16bitPixels( short [] in );

    }
   
    // =======================================================
    
    /** Copy a region of vector xIn into this vector, generic implementation.
     *  @param xIn Input vector
     *  @param yOut Vector to paste into
     *  @param ulxin x-coordinate upper left corner of region to be copied
     *  @param ulyin y-coordinate upper left corner of region to be copied
     *  @param w width of region to be copied
     *  @param h height of region to be copied
     *  @param ulxout x-coordinate upper left corner of output position
     *  @param ulyout y-coordinate upper left corner of output position 
     *  @param zeroOutput If true, vector gets zero'd before pasting data
     *  */
    public static void paste( Vec2d.Real xIn, Vec2d.Real yOut,
	final int ulxin, final int ulyin, final int w, final int h,
	final int ulxout, final int ulyout, boolean zeroOutput ) {

	float [] in  =  xIn.vectorData();
	float [] out = yOut.vectorData();

	// Checks:
	if (ulxin<0 || ulyin<0 || w<0 || h<0 || ulyout<0 || ulxout<0 )
	    throw new RuntimeException("Coordinates not positive");
	if ( ( ulxin + w > xIn.vectorWidth() ) ||
	     ( ulyin + h > xIn.vectorHeight() ) )
	    throw new RuntimeException("Input size mismatch");
	if ( ( ulxout + w  > yOut.vectorWidth() ) ||
	     ( ulyout + h > yOut.vectorHeight() ) )
	    throw new RuntimeException("Output size mismatch");
	
	// zero
	if ( zeroOutput)
	    yOut.zero();
	//Tool.trace(""+ulyout+" "+ulxout+" "+ulyin+" "+ulxin+" "+w+" "+h);
	
	final int outWidth = yOut.vectorWidth();
	final int inWidth  = xIn.vectorWidth();

	// copy
	for (int y=0;y<h;y++)
	for (int x=0;x<w;x++)
	    out[ (y + ulyout) * outWidth + ( x + ulxout ) ] =
		in[ (y + ulyin) * inWidth + ( x + ulxin ) ];

	yOut.syncBuffer();
    }
    
    /** Copy a region of vector xIn into this vector, generic implementation.
     *  @param xIn Input vector
     *  @param yOut Vector to paste into
     *  @param ulxin x-coordinate upper left corner of region to be copied
     *  @param ulyin y-coordinate upper left corner of region to be copied
     *  @param w width of region to be copied
     *  @param h height of region to be copied
     *  @param ulxout x-coordinate upper left corner of output position
     *  @param ulyout y-coordinate upper left corner of output position 
     *  @param zeroOutput If true, vector gets zero'd before pasting data
     *  */
    public static void paste( Vec2d.Cplx xIn, Vec2d.Cplx yOut,
	final int ulxin, final int ulyin, final int w, final int h,
	final int ulxout, final int ulyout, boolean zeroOutput ) {

	float [] in  =  xIn.vectorData();
	float [] out = yOut.vectorData();

	// Checks:
	if (ulxin<0 || ulyin<0 || w<0 || h<0 || ulyout<0 || ulxout<0 )
	    throw new RuntimeException("Coordinates not positive");
	if ( ( ulxin + w > xIn.vectorWidth() ) ||
	     ( ulyin + h > xIn.vectorHeight() ) )
	    throw new RuntimeException("Input size mismatch");
	if ( ( ulxout + w  > yOut.vectorWidth() ) ||
	     ( ulyout + h > yOut.vectorHeight() ) )
	    throw new RuntimeException("Output size mismatch");
	
	// zero
	if ( zeroOutput)
	    yOut.zero();
	//Tool.trace(""+ulyout+" "+ulxout+" "+ulyin+" "+ulxin+" "+w+" "+h);
	
	final int outWidth = yOut.vectorWidth();
	final int inWidth  = xIn.vectorWidth();

	// copy
	for (int y=0;y<h;y++)
	for (int x=0;x<w;x++){
	    out[ ((y + ulyout) * outWidth + ( x + ulxout ))*2+0 ] =
		in[ ((y + ulyin) * inWidth + ( x + ulxin ))*2+0 ];
	    out[ ((y + ulyout) * outWidth + ( x + ulxout ))*2+1 ] =
		in[ ((y + ulyin) * inWidth + ( x + ulxin ))*2+1 ];

	}

	yOut.syncBuffer();
    }
    



} public final class MTool {



    /** Computes sin(x), precision 5.4e-9, approx 3x faster than Math.sin. */
    public static double fsin(double x) {
	x  = x/(Math.PI/2);	    // 0..2pi -> 0..4
	x -= Math.floor(0.25*x)*4;  // x now in 0..4

	// check which quadrat to return
	if (x>3) return -aproxSin01(4-x);
	if (x>1) return  aproxSin01(2-x);
	return aproxSin01(x);
    }
    
    /** Computes cos(x), precision 5.4e-9, approx 3x faster than Math.sin */
    public static double fcos(double x) {
	return fsin(x+Math.PI/2);
    }
    
    // http://stackoverflow.com/questions/523531/fast-transcendent-trigonometric-functions-for-java
    private static double aproxSin01(double x) {
  	double x2 = x * x;
  	return ((((.00015148419 * x2
		- .00467376557) * x2
            	+ .07968967928) * x2
           	- .64596371106) * x2
          	+ 1.57079631847) * x;
    }	

    /** Computes sqrt(x^2+y^2), without the overflow handling of Math.hypot(x,y).
     *  Math.hypot ensures correct results even if x^2 would overflow a double.
     *  If neither x^2 nor y^2 overflows, this will be 30x faster */
    public static double fhypot(double x, double y) {
	return Math.sqrt(x*x+y*y);
    }
    

    /** Returns the average of an array of doubles */
    public static double average( double [] in ) {
	double ret = 0;
	for ( double i : in ) ret += i/in.length;
	return ret;
    }

    /** Returns the average of an array of float.
     *	This uses an double internally as accumulator for higher precision. */
    public static float average( float [] in ) {
	double ret = 0;
	for ( float i : in ) ret += i/in.length;
	return (float)ret;
    }

    /** Returns greatest common divisor of a and b.
     * Uses Euclid's algorithm. */
    public static int greatestCommonDivisor(int a, int b)
    {
	while (b > 0)
	{
	    int c = b;
	    b = a % b; 
	    a = c;
	}
	return a;
    }



    /** Tests against java.lang.Math, Timings */
    public static void main(String [] args ) {

	// check approx
	double sindiv = 0;
	for (double r=-6*Math.PI;r<6*Math.PI;r+=0.01) {
	    double ds = Math.abs( fsin(r) - Math.sin(r) );
	    double dc = Math.abs( fcos(r) - Math.cos(r) );
	    //System.out.println(fsin(r)+" "+Math.sin(r));
	    if (ds>sindiv) sindiv=ds;
	    if (dc>sindiv) sindiv=ds;
	}
	System.out.println(String.format("sin/cos: Largest abs. error: %5.3e",sindiv));

	// timing
	double [] v1, v2;
	float [] v2f ;	
	v1 = new double[ 1000000 ];
	v2 = new double[ 1000000 ];
	Tool.Timer t1 = Tool.getTimer();
	Tool.Timer t2 = Tool.getTimer();
  
	// sin, cos
	for (int loop=0; loop<4; loop++) {
	    t1.start();
	    for ( int i=0;i< v1.length; i++)
		v1[i] = Math.sin( (i/100.) );
	    t1.stop();
	    
	    t2.start();
	    for ( int i=0;i< v2.length; i++)
		v2[i] = fsin( (i/100.) );
	    t2.stop();

	    System.out.println("Timing,   Math.sin vs fsin:   "+t1+" "+t2);
	}   
	
	// hypot
	for (int loop=0; loop<4; loop++) {
	    t1.start();
	    for ( int i=0;i< v1.length; i++)
		v1[i] = Math.hypot( (i/100.), (i%100) );
	    t1.stop();
	    
	    t2.start();
	    for ( int i=0;i< v2.length; i++)
		v2[i] = fhypot( (i/100.), (i%100) );
	    t2.stop();

	    System.out.println("Timing, Math.hypot vs fhypot: "+t1+" "+t2);
	}   
	
	
    }



} public abstract class Transforms {


    /** Two-dimensional FFT of the input vector. */
    static public void fft2d( Vec2d.Cplx in, boolean inverse ) {
	in.fft2d(inverse);
    }
    
    /** Three-dimensional FFT of the input vector. */
    static public void fft3d( Vec3d.Cplx in, boolean inverse ) {
	in.fft3d(inverse);
    }





    static public void runfft( Vec2d.Cplx in, boolean inverse ) {
	// get parameters
	final int w = in.vectorWidth();
	final int h = in.vectorHeight();
	// see if we have an instance already, otherwise make one
	FftProvider.Instance ffti = getOrCreateInstance(new FFTkey(w,h));
	float [] dat = in.vectorData();
	ffti.fftTransform( dat , inverse );
	in.syncBuffer();
    }

    static public void runfft3d( Vec3d.Cplx in, boolean inverse ) {
	// get parameters
	final int w = in.vectorWidth();
	final int h = in.vectorHeight();
	final int d = in.vectorDepth();
	// see if we have an instance already, otherwise make one
	FftProvider.Instance ffti = getOrCreateInstance(new FFTkey(w,h,d));
	float [] dat = in.vectorData();
	ffti.fftTransform( dat , inverse );
	in.syncBuffer();
    }



    /** One-dimensional FFT of the complex input vector. */
    static public void fft1d( Vec.Cplx in, boolean inverse ) {
	final int len = in.vectorSize();
	FftProvider.Instance ffti = getOrCreateInstance( new FFTkey(len));
	float [] dat = in.vectorData();
	ffti.fftTransform( dat, inverse );
	in.syncBuffer();
    }

    /** One-dimensional FFT of a standard float array. This assumes
     *  standard packing ( a[i*2] = real[i], a[i*2+1] =imag[i] ) of
     *  the input array */
    static public void fft1d( float [] in, boolean inverse ) {
	final int len = in.length/2;
	FftProvider.Instance ffti = getOrCreateInstance( new FFTkey(len));
	ffti.fftTransform( in, inverse );
    }

    // -----------------------------------------------------------
    // Instance management


    /** key to store instances */
    private static class FFTkey implements  Comparable<FFTkey> { 
	final int d,x,y,z ; 
	FFTkey( int xi ) {
	    d=1; x=xi; y=-1; z=-1;
	}
	FFTkey( int xi, int yi ) { 
	    d=2; x=xi; y=yi; z=-1;
	}
	FFTkey( int xi, int yi, int zi ) { 
	    d=3; x=xi; y=yi; z=zi;
	}
	@Override
	public int compareTo(FFTkey t) {
	    if (d != t.d) return (t.d -d );
	    if (x != t.x) return (t.x -x );
	    if (y != t.y) return (t.y -y );
	    if (z != t.z) return (t.z -z );
	    return 0;
	}
    }
    
    /** FFT instances */
    static private Map<FFTkey, FftProvider.Instance> instances; 
    
    /** static initialization of the instances list. */
    static {
	if (instances==null) {
	    instances = new TreeMap<FFTkey, FftProvider.Instance>();
	}
    }



    /** returns an instance, creates one if none exists */
    static protected FftProvider.Instance getOrCreateInstance(final FFTkey k) {
	FftProvider.Instance ffti = instances.get(k);
	if (ffti!=null) return ffti;
	//Tool.trace("FFT: creating new instance");
	if (k.d==1)
	    ffti = FftProvider.get1Dfft(k.x);
	if (k.d==2)
	    ffti = FftProvider.get2Dfft(k.x,k.y);
	if (k.d==3)
	    ffti = FftProvider.get3Dfft(k.x,k.y,k.z);
	if (ffti==null) 
	    throw new RuntimeException("Unsupported dimensions");
	instances.put( k , ffti );
	return ffti;
    }


    // ---------------------------------------------------------
    //
    // abtract functions for the actual transformations
    //
    // ---------------------------------------------------------

    /** 
     *  1D Fourier tranformation complexFloat2complexFloat.
     *  This has to be implemented in a subclass. For convenience,
     *  all necessary parameters are passed.
     */
    protected abstract void fft_1d_trans_c2c(  float [] x, boolean inverse );
    /** 
     *  2D Fourier tranformation complexFloat2complexFloat.
     *  This has to be implemented in a subclass. For convenience,
     *  all necessary parameters are passed.
     */
    protected abstract void fft_2d_trans_c2c(  float [] x, boolean inverse );

    /** 
     *  3D Fourier tranformation complexFloat2complexFloat.
     *  This has to be implemented in a subclass. For convenience,
     *  all necessary parameters are passed.
     */
    protected abstract void fft_3d_trans_c2c(  float [] x, boolean inverse );



    /** checks if input is 2^n */
    static public boolean powerOf2(int l) {
	    int i=2;
	    while(i<l) i *= 2;
	    return (i==l);
    }
    
    // ---------------------------------------------------------
    //
    // Power spectrum calculations
    //
    // ---------------------------------------------------------

    /** Computes a power spectrum, assuming the input vector is an FFT spectrum.
     *	Note: The spectrum will be clipped and quadrant-swapped for display purposes. */
    static public void computePowerSpectrum( Vec2d.Cplx inV, Vec2d.Real outV) {
	    computePowerSpectrum( inV, outV, true );
    }
    
    /** Computes a power spectrum, assuming the input vector is an FFT spectrum.
     *  Note: The spectrum is clipped to log(max)-log(min)<30, and can be quadrand-swapped. */
    static public void computePowerSpectrum( Vec2d.Cplx inV, Vec2d.Real outV, 
	boolean swap_quad ) {


	// TODO: maybe speed this up to higher efficience
	final int w= inV.vectorWidth();
	final int h= inV.vectorHeight();
	float [] in  =  inV.vectorData();
	float [] out  = outV.vectorData();

	if (( w != outV.vectorWidth()  ) || ( h != outV.vectorHeight() ))
	    throw new RuntimeException("Vector size mismatch");



	// calculate min and max for scaling
	float min = Float.MAX_VALUE;
	float max = Float.MIN_VALUE;
	for ( float i : in ) {
	    if (i<min) min=i;
	    if (i>max) max=i;
	}
	
	// reduce the range 
	max = (float)Math.log(max);
	min = (float)Math.log(min);
	if (Float.isNaN(min) || max-min>30)
	    min = max-30;

	
	for (int y=0;y<h;y++)
	for (int x=0;x<w;x++) {
	    float r = (float)Math.sqrt(
		Math.pow( in[2*(y*w+x)+0] , 2)+Math.pow( in[2*(y*w+x)+1] , 2));
	    r = (float)(((Math.log(r) - min)/(max-min)));
	    if (Float.isNaN(r) || r<0) r=0f;

	    if (swap_quad) {
		int xo = (x<w/2)?(x+w/2):(x-w/2);
		int yo = (y<h/2)?(y+h/2):(y-h/2);
		out[    yo*w    + xo ] = r;
	    } else {
		out[    y*w    + x ] = r;
	    }
	}
    
	outV.syncBuffer( );

    }

    /** Swap quadrands. */
    static public void swapQuadrant(Vec2d.Cplx in) {
	final int w = in.vectorWidth();
	final int h = in.vectorHeight();
	for (int y=0;y<h/2;y++)
	for (int x=0;x<w/2;x++) {
	    // 1 <-> 3
	    Cplx.Float tmp1 = in.get(x,y);
	    Cplx.Float tmp3 = in.get(x+w/2,y+h/2);
	    in.set(x,y,tmp3);
	    in.set(x+w/2,y+h/2,tmp1);
	    // 2 <-> 3
	    Cplx.Float tmp2 = in.get(x,y+h/2);
	    Cplx.Float tmp4 = in.get(x+w/2,y);
	    in.set(x,y+h/2,tmp4);
	    in.set(x+w/2,y,tmp2);
	}
    }
    
    /** Swap quadrands. */
    static public void swapQuadrant(Vec2d.Real in) {
	final int w = in.vectorWidth();
	final int h = in.vectorHeight();
	for (int y=0;y<h/2;y++)
	for (int x=0;x<w/2;x++) {
	    // 1 <-> 3
	    float tmp1 = in.get(x,y);
	    float tmp3 = in.get(x+w/2,y+h/2);
	    in.set(x,y,tmp3);
	    in.set(x+w/2,y+h/2,tmp1);
	    // 2 <-> 3
	    float tmp2 = in.get(x,y+h/2);
	    float tmp4 = in.get(x+w/2,y);
	    in.set(x,y+h/2,tmp4);
	    in.set(x+w/2,y,tmp2);
	}
    }
    
    /** Swap quadrands. */
    static public void swapQuadrant(Vec3d.Cplx in) {
	final int w = in.vectorWidth();
	final int h = in.vectorHeight();
	final int d = in.vectorHeight();
	for (int z0=0;z0<d/2;z0++)
	for (int y0=0;y0<h/2;y0++)
	for (int x0=0;x0<w/2;x0++) {
	    int x1=x0+w/2, y1=y0+h/2, z1=z0+d/2;
	    Cplx.Float tmpA, tmpB;
	    // 000 <> 111
	    tmpA = in.get(x0,y0,z0);
	    tmpB = in.get(x1,y1,z1);
	    in.set(x0,y0,z0,tmpB);
	    in.set(x1,y1,z1,tmpA);
	    // 010 <> 101
	    tmpA = in.get(x0,y1,z0);
	    tmpB = in.get(x1,y0,z1);
	    in.set(x0,y1,z0,tmpB);
	    in.set(x1,y0,z1,tmpA);
	    // 100 <> 011
	    tmpA = in.get(x1,y0,z0);
	    tmpB = in.get(x0,y1,z1);
	    in.set(x1,y0,z0,tmpA);
	    in.set(x0,y1,z1,tmpB);
	    // 001 <> 110
	    tmpA = in.get(x0,y0,z1);
	    tmpB = in.get(x1,y1,z0);
	    in.set(x0,y0,z1,tmpB);
	    in.set(x1,y1,z0,tmpA);
	}
    }
    
    /** Swap quadrands. */
    static public void swapQuadrant(Vec3d.Real in) {
	final int w = in.vectorWidth();
	final int h = in.vectorHeight();
	final int d = in.vectorHeight();
	for (int z0=0;z0<d/2;z0++)
	for (int y0=0;y0<h/2;y0++)
	for (int x0=0;x0<w/2;x0++) {
	    int x1=x0+w/2, y1=y0+h/2, z1=z0+d/2;
	    float tmpA, tmpB;
	    // 000 <> 111
	    tmpA = in.get(x0,y0,z0);
	    tmpB = in.get(x1,y1,z1);
	    in.set(x0,y0,z0,tmpB);
	    in.set(x1,y1,z1,tmpA);
	    // 010 <> 101
	    tmpA = in.get(x0,y1,z0);
	    tmpB = in.get(x1,y0,z1);
	    in.set(x0,y1,z0,tmpB);
	    in.set(x1,y0,z1,tmpA);
	    // 100 <> 011
	    tmpA = in.get(x1,y0,z0);
	    tmpB = in.get(x0,y1,z1);
	    in.set(x1,y0,z0,tmpA);
	    in.set(x0,y1,z1,tmpB);
	    // 001 <> 110
	    tmpA = in.get(x0,y0,z1);
	    tmpB = in.get(x1,y1,z0);
	    in.set(x0,y0,z1,tmpB);
	    in.set(x1,y1,z0,tmpA);
	}
    }
    
    
    
    /* Return a vector containing phases for a Fourier shift theorems shift to kx,ky.
     *  @param N Width and heigt of vector
     *  @param kx x-coordinate of shift
     *  @param ky y-coordinate of shift
     *  @param fast Use faster, but less precise sin/cos (see {@link MTool#fsin}) */
    /*
    static public Vec2d.Cplx createShiftVector( 
	final int N, final double kx, final double ky, final boolean fast ) {
	Vec2d.Cplx shft = Vec2d.createCplx(N,N);
	final float [] val = shft.vectorData();
	
	// run outer loop in parallel
	new SimpleMT.PFor(0,N) {
	    public void at(int y) {
		for (int x=0; x<N; x++) {
		    float phaVal = (float)(2*Math.PI*(kx*x+ky*y)/N);
		    if (fast) {
			val[ (y*N+x)*2+0 ] = (float)MTool.fcos( phaVal );
			val[ (y*N+x)*2+1 ] = (float)MTool.fsin( phaVal );
		    } else {
			val[ (y*N+x)*2+0 ] = (float)Math.cos( phaVal );
			val[ (y*N+x)*2+1 ] = (float)Math.sin( phaVal );
		    }
		}
	    }
	};

	shft.syncBuffer();

	return shft;
    } */

    /* See {@link #createShiftVector}, with fast='false' */
    /*
    static public Vec2d.Cplx createShiftVector(int N, double kx, double ky ) {
	return createShiftVector(N, kx, ky, false);
    } */

    /** Multiply a vector with Fourier shift theorem phases.
     *  Vector has to be of square size (w==h).
     *  @param kx x-coordinate of shift
     *  @param ky y-coordinate of shift
     *  @param fast Use faster, but less precise sin/cos (see {@link MTool#fsin}) */
    @Deprecated
    static public void timesShiftVector( final Vec2d.Cplx vec,
	final double kx, final double ky, final boolean fast ) {
	vec.fourierShift( kx, ky );
    }


    static public void runTimesShiftVector( final Vec2d.Cplx vec,
	final double kx, final double ky, final boolean fast ) {
	final float [] val = vec.vectorData();
	final int N = Vec2d.checkSquare( vec );
	
	// run outer loop in parallel
	new SimpleMT.PFor(0,N) {
	    public void at(int y) {
		for (int x=0; x<N; x++) {
		    float phaVal = (float)(2*Math.PI*(kx*x+ky*y)/N);
		    float si,co;
		    if (fast) {
			co = (float)MTool.fcos( phaVal );
			si = (float)MTool.fsin( phaVal );
		    } else {
			co = (float)Math.cos( phaVal );
			si = (float)Math.sin( phaVal );
		    }
		    // get
		    float re = val[ (y*N+x)*2+0 ] ;
		    float im = val[ (y*N+x)*2+1 ] ;
		    // set
		    val[ (y*N+x)*2+0 ] = Cplx.multReal( re, im, co, si );
		    val[ (y*N+x)*2+1 ] = Cplx.multImag( re, im, co, si );
		}
	    }
	};

	vec.syncBuffer();
    }

    /** See {@link #timesShiftVector}, with fast='false' */
    static public void timesShiftVector( final Vec2d.Cplx vec, double kx,  double ky ) {
	timesShiftVector( vec, kx, ky, false );
    }



	
}
import java.util.List;
import java.util.ArrayList;
import javax.swing.JPanel;
import javax.swing.JDialog;
import javax.swing.JFrame;
import javax.swing.JButton;
import javax.swing.JLabel;
import javax.swing.JCheckBox;
import javax.swing.JOptionPane;
import javax.swing.SwingWorker;
import javax.swing.JSlider;
import javax.swing.event.ChangeListener;
import javax.swing.event.ChangeEvent;
import javax.swing.JProgressBar;
import javax.swing.JTabbedPane;
import javax.swing.Box;
import javax.swing.BoxLayout;
import javax.swing.BorderFactory;
import java.awt.Dimension;
import java.awt.ComponentOrientation;
import java.awt.Component;
import java.awt.Color;
import java.awt.event.WindowAdapter;
import java.awt.event.WindowEvent;
import java.awt.GridBagLayout;
import java.awt.GridBagConstraints;
import java.awt.event.ActionListener;
import java.awt.event.ActionEvent;
import javax.swing.event.PopupMenuListener;
import javax.swing.event.PopupMenuEvent;
import org.fairsim.utils.Tool;
import org.fairsim.utils.ImageSelector;
import org.fairsim.utils.ImageDisplay;
import org.fairsim.sim_algorithm.OtfProvider;
import org.fairsim.sim_algorithm.SimParam;
import org.fairsim.sim_algorithm.SimUtils;
import org.fairsim.sim_algorithm.SimAlgorithm;
import org.fairsim.linalg.Vec2d;
import org.fairsim.linalg.MTool;
import org.fairsim.linalg.Transforms;
/*
This file is part of Free Analysis and Interactive Reconstruction
for Structured Illumination Microscopy (fairSIM).

fairSIM 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, either version 2 of the License, or
(at your option) any later version.

fairSIM 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.

You should have received a copy of the GNU General Public License
along with fairSIM.  If not, see <http://www.gnu.org/licenses/>
*/

package org.fairsim.sim_gui;










/**
 * GUI elements to set up OTF
 **/
public class ImageControl { 
    
    private final JPanel ourContent = new JPanel();  
    private final JLabel ourState = new JLabel("no image selected");
    private volatile boolean running = false;

    private final ImageSelector imgSelect ;
    //private final ImageSelector.Callback imgSelectCB;
    
    private final JFrame baseframe;
    private final SimParam simParam;
    private final FairSimGUI fsGUI;
    private ImageDisplay.Factory idpFactory ;

    // global references to be accessed both by importImage and showSlices
    private JButton importImageButton	= null;
    private JCheckBox importTimelapseBox= null;
    private JButton batchButton		= null;
    private JButton showSlicesButton	= null;
    private ImageDisplay sliceSelector	= null;
    private Tiles.LComboBox<Integer> sliceBox;
	
    final Tiles.LComboBox<ImageSelector.ImageInfo> imgBox;
    
    private Tiles.LComboBox<String> bgrBox;
    private Tiles.LNSpinner bgrSpinner;

    // video position slider
    private final JSlider videoPosSlider;
    final JLabel  videoPosLabel = new JLabel( String.format("t % 5d", 0));
    final JLabel maxTimePointsLabel;
    final Tiles.LNSpinner zSliceVideoSpinner ; 
    final String [] updateModes = {"off", "widefield", "recon", "recon+par.est."};
    final Tiles.LComboBox<String> videoAutoUpdateMode =
     new Tiles.LComboBox<String>("auto-update", updateModes); 

    // prefactor correction spinners and auto selector
    Tiles.LNSpinner  []   prefactorAngSpinner; 
    Tiles.LNSpinner [][] prefactorPhaSpinner; 
    private JCheckBox prefactorAutoUpdateAngle;
    private JCheckBox prefactorAutoUpdatePhase;
    final Tiles.LComboBox<String> prefactorMethodBox
	= new Tiles.LComboBox<String>("estimate, by", "average", "median"); 


    // the images
    Vec2d.Real [][] theImages    =null;
    Vec2d.Cplx [][] theFFTImages =null;
    double pxlSize;
    boolean pxlSet;
    
    private ImageDisplay rawDataDisplay   = null;
    private ImageDisplay simPreviewDisplay = null;

    // video functions
    boolean isVideoStack = false;
    int videoStackPositionTime = 0;
    int videoStackPositionZ =0;
    boolean sucessfulImport = false;
    int totalTimePoints = 1;

    AutoUpdateThread autoUpdater = null;

    public JPanel getPanel() {
	return ourContent;
    }

    /** Contructor, initializes image list. */
    public ImageControl( final JFrame baseframe, ImageSelector is, 
	final ImageDisplay.Factory imgFactory, 	FairSimGUI fsg, final SimParam sp ) {

	// initialize variables
	this.baseframe    = baseframe;
	this.imgSelect    = is;
	this.simParam     = sp;
	this.idpFactory   = imgFactory;
	this.fsGUI        = fsg;

	JTabbedPane imgTabs = new JTabbedPane();

	// create our panel
	ourContent.setLayout(new BoxLayout(ourContent, BoxLayout.PAGE_AXIS));
	ourContent.setBorder(BorderFactory.createTitledBorder("1 - Image Selector") );

	// setup label
	ourState.setAlignmentX(Component.CENTER_ALIGNMENT);
	ourState.setForeground(Color.RED);

	// setup selector box
	ImageSelector.ImageInfo [] openImages = imgSelect.getOpenImages();
	imgBox = new Tiles.LComboBox<ImageSelector.ImageInfo>(
	    "Img", null, true, openImages);
	importImageButton  = new JButton("import");
	importTimelapseBox  = new JCheckBox("timelapse mode");
	
	zSliceVideoSpinner = 	
	    new Tiles.LNSpinner( "zSlices" , 1, 1, 200, 1 ); 

	imgBox.box.setPrototypeDisplayValue("some really long image filename"+
	    "as these will often be used");

	maxTimePointsLabel = new JLabel("t-max: n/a");
	batchButton = new JButton("batch");

	JPanel row1 = new JPanel();
	row1.setLayout(new BoxLayout(row1, BoxLayout.LINE_AXIS));
	row1.add(Box.createHorizontalGlue());
	row1.add( imgBox );
	row1.add(Box.createRigidArea(new Dimension(5,1)));
	row1.add( importImageButton );
	
	JPanel row2 = new JPanel();
	row2.add(Box.createHorizontalGlue());
	row2.add( importTimelapseBox );
	row2.add( zSliceVideoSpinner);
	row2.add(Box.createRigidArea(new Dimension(5,1)));
	row2.add( maxTimePointsLabel );
	row2.add(Box.createRigidArea(new Dimension(5,1)));
	row2.add( videoAutoUpdateMode );
	row2.add(Box.createHorizontalGlue());
	
	// add the video-select sliders
	JPanel row3 = new JPanel();
	
	videoPosSlider = new JSlider(JSlider.HORIZONTAL, 1,100,1);
	videoPosSlider.setEnabled(false);

	JPanel sliders = new JPanel(new GridBagLayout());
	GridBagConstraints c = new GridBagConstraints();	
	
	c.gridx=0; c.gridy=0; c.gridwidth=10; c.gridheight=1;
	sliders.add( videoPosSlider,c);
	c.gridx=11; c.gridwidth=1;
	sliders.add( videoPosLabel,c );
	c.gridx=12; c.gridwidth=1;
	sliders.add( Box.createRigidArea( new Dimension(5,1) ));
	c.gridx=13; c.gridwidth=1;
	sliders.add( batchButton );


	// preprocessing control
	JPanel preprocessPanel = new JPanel(new GridBagLayout());
	//preprocessPanel.setBorder(
	//    BorderFactory.createTitledBorder("Intensity correction pre-processing") );

	c.gridwidth=1; c.gridheight=1;

	prefactorAngSpinner = new Tiles.LNSpinner[ sp.nrDir() ];
	prefactorPhaSpinner = new Tiles.LNSpinner[ sp.nrDir() ][]; 

	for (int ang = 0 ; ang < sp.nrDir(); ang++) {
	    
	    prefactorAngSpinner[ang]  = new Tiles.LNSpinner(String.format("a%1d",ang),
		1, .01, 100, 0.01);
	  
	    final int angCopy = ang;
	    
	    prefactorAngSpinner[ang].addNumberListener( new Tiles.NumberListener() {
		@Override
		public void number( double n, Tiles.LNSpinner e) {
		    sp.dir(angCopy).setAngleIntensityFactor( n );	
		};
	    });

	    c.gridx=0; c.gridy=ang+1;

	    preprocessPanel.add( prefactorAngSpinner[ang], c); 

	    prefactorPhaSpinner[ang] = new Tiles.LNSpinner[ sp.dir(ang).nrPha()];
	    for (int pha = 0 ; pha < sp.dir(ang).nrPha(); pha++) {
		prefactorPhaSpinner[ang][pha] = new Tiles.LNSpinner(String.format("p%1d",pha),
		1, .01, 100, 0.01);
	    
		final int phaCopy = pha;
		
		prefactorPhaSpinner[ang][pha].addNumberListener( new Tiles.NumberListener() {
		    @Override
		    public void number( double n, Tiles.LNSpinner e) {
			sp.dir(angCopy).setPhaseIntensityFactor( phaCopy, n );	
		    };
	        });

	    
		c.gridx=pha+1;
		preprocessPanel.add( prefactorPhaSpinner[ang][pha], c); 
	    
	    }
	}

	c.gridx=0; c.gridy=0; c.gridwidth=3; 
	preprocessPanel.add( prefactorMethodBox, c );

	c.gridx=3; c.gridwidth=1;
	
	JButton estAngleButton = new JButton("angles");
	preprocessPanel.add( estAngleButton, c );
	
	estAngleButton.addActionListener( new ActionListener() {
	    public void actionPerformed(ActionEvent e) {
		runEstimateAngleVariation(prefactorMethodBox.getSelectedIndex());
	    }
	});

	c.gridx=4;
	JButton estPhaseButton = new JButton("phases");
	preprocessPanel.add( estPhaseButton, c );
	estPhaseButton.addActionListener( new ActionListener() {
	    public void actionPerformed(ActionEvent e) {
		runEstimatePhaseVariation(prefactorMethodBox.getSelectedIndex());
	    }
	});

	c.gridx=0; c.gridy=sp.nrDir()+1;

	JButton showCorrectedImagesButton = new JButton("show corrected stack");
	preprocessPanel.add( showCorrectedImagesButton, c);
	
	showCorrectedImagesButton.addActionListener( new ActionListener() {
	    @Override
	    public void actionPerformed(ActionEvent e) {
		if (!sucessfulImport || theImages==null) {
		    JOptionPane.showMessageDialog( baseframe,
		    "Please import an image first", "fairSIM error",
		     JOptionPane.ERROR_MESSAGE);
		    return;
		}
		
		ImageDisplay correctedDataDisplay = 
		    idpFactory.create(theImages[0][0].vectorWidth(), theImages[0][0].vectorHeight(), 
			"Intensity corrected input images");

		for (int ang=0; ang<sp.nrDir(); ang++) {
		    for (int pha=0; pha<sp.dir(ang).nrPha(); pha++) {

			Vec2d.Real corImg = theImages[ang][pha].duplicate();
			corImg.scal( (float)sp.dir(ang).getIntensityQuotient(pha));
			correctedDataDisplay.addImage( corImg, String.format("a%d p%d",ang,pha));
		    }
		
		}

		correctedDataDisplay.display();


	    }
	});

	JButton resetValuesButton = new JButton("reset");
	c.gridx=1;
	preprocessPanel.add( resetValuesButton,c );
	
	c.gridy++;
	prefactorAutoUpdateAngle = new JCheckBox("auto-update angle");
	prefactorAutoUpdatePhase = new JCheckBox("phase");
	c.gridx=0;
	preprocessPanel.add( prefactorAutoUpdateAngle,c );
	c.gridx=2;
	preprocessPanel.add( prefactorAutoUpdatePhase,c );

	resetValuesButton.addActionListener( new ActionListener() {
	    @Override
	    public void actionPerformed(ActionEvent e) {
		for (int ang=0; ang<sp.nrDir(); ang++) {
		   prefactorAngSpinner[ang].spr.setValue(1);
		   for (int pha=0; pha<sp.dir(ang).nrPha(); pha++) {
			prefactorPhaSpinner[ang][pha].spr.setValue(1);
		   }
		}
	    }
	});

	
	// add all content to the panel
	JPanel tab1 = new JPanel();
	JPanel tab2 = new JPanel();
	JPanel tab3 = new JPanel();

	tab1.setLayout(new BoxLayout(tab1, BoxLayout.PAGE_AXIS));
	tab2.setLayout(new BoxLayout(tab2, BoxLayout.PAGE_AXIS));
	tab3.setLayout(new BoxLayout(tab3, BoxLayout.PAGE_AXIS));

	tab1.add( Box.createVerticalGlue());
	tab1.add(ourState);
	tab1.add(Box.createRigidArea(new Dimension(1,5)));
	tab1.add(row1);
	tab1.add( Box.createVerticalGlue());
	
	tab2.add( Box.createVerticalGlue());
	tab2.add(row2);
	tab2.add(Box.createRigidArea(new Dimension(1,5)));
	tab2.add(sliders);
	tab2.add( Box.createVerticalGlue());
	//ourContent.add(Box.createRigidArea(new Dimension(1,5)));


	tab3.add( preprocessPanel );
	tab3.add(Box.createRigidArea(new Dimension(1,5)));

	imgTabs.addTab("Load", tab1);
	imgTabs.addTab("Timelapse", tab2);
	imgTabs.addTab("Corrections", tab3);
    
	ourContent.add( imgTabs);

	// ------ CALLBACKS and LOGIC ------
	// refresh when opened
	imgBox.box.addPopupMenuListener( new PopupMenuListener() {
	    public void popupMenuCanceled( PopupMenuEvent e ) {};
	    public void popupMenuWillBecomeInvisible( PopupMenuEvent e ) {};
	    public void popupMenuWillBecomeVisible( PopupMenuEvent e ) {
		//System.out.println("note");	
		imgBox.newElements( imgSelect.getOpenImages());
	    };
	});


	// import images / time-lapse
	importImageButton.addActionListener( new ActionListener() {
	    public void actionPerformed( ActionEvent e ) {
		
		if ( isVideoStack ) {
		    // import as time-laspe stack
		    importImageDialog( imgBox.getSelectedItem(), 
			(int)zSliceVideoSpinner.getVal(), videoStackPositionTime );
		    isVideoStack = true;
		} else {
		    // import as stadard stack
		    importImageDialog( 
			imgBox.getSelectedItem(), 
			imgBox.getSelectedItem().depth / simParam.getImgPerZ(), 0 );
		
		    isVideoStack = false;
		}
	    }
	});
	
	// change import mode
	importTimelapseBox.addActionListener( new ActionListener() {
	    public void actionPerformed( ActionEvent e ) {
		changeTimelapseMode( importTimelapseBox.isSelected());
	    }
	});

	// video position slider change position
	videoPosSlider.addChangeListener( new ChangeListener() {
	    public void stateChanged(ChangeEvent e) {
		videoStackPositionTime = videoPosSlider.getValue()-1;
		videoPosLabel.setText( 
		    String.format(" t:  % 5d", videoStackPositionTime+1));
		simParam.signalRuntimeChange();
	    }
	});

	// selected image changed
	imgBox.addSelectListener( 
	    new Tiles.SelectListener<ImageSelector.ImageInfo>() {
	    public void selected( ImageSelector.ImageInfo inf, int i ) {
		updateIfVideoCanBeUsed();	
	    }
	});

	// selected number of z-planes changed
	zSliceVideoSpinner.addNumberListener( new Tiles.NumberListener() {
	    public void number(double d, Tiles.LNSpinner e) {
		updateIfVideoCanBeUsed();
	    }
	});

	// change the auto-update mode
	videoAutoUpdateMode.addSelectListener( new Tiles.SelectListener<String>() {
	    public void selected( String e, int i) {
		
		if (!sucessfulImport) {
		    JOptionPane.showMessageDialog( baseframe,
		     "Please run a successful image import first", "fairSIM error",
		     JOptionPane.ERROR_MESSAGE);
		    videoAutoUpdateMode.setSelectedIndex(0);
		    return;
		}
		
		if (i>1 && simParam.otf() == null) {
		    JOptionPane.showMessageDialog( baseframe,
		     "Please run a successful single SIM reconstruction first", 
		     "fairSIM error",
		     JOptionPane.ERROR_MESSAGE);
		    videoAutoUpdateMode.setSelectedIndex(0);
		    return;
		}
		
		runAutoUpdate(i);
	    
	    }
	});

	// run the batch reconstruction dialog
	batchButton.addActionListener( new ActionListener() {
	    public void actionPerformed( ActionEvent e) {
		if (!sucessfulImport || simParam.otf() == null) {
		    JOptionPane.showMessageDialog( baseframe,
		     "Please run a successful single SIM reconstruction first", 
		     "fairSIM error",
		     JOptionPane.ERROR_MESSAGE);
		    return;
		}
		openBatchDialog();
	    }
	});


	// update the video function
	changeTimelapseMode(false);
    }

    /** check if the image is o.k. */
    boolean checkImage( ImageSelector.ImageInfo img, int zSlices ) {
	if ((img==null)||(img.depth==0)) {
	    JOptionPane.showMessageDialog( baseframe,
	     "No image selected", "fairSIM error",
	     JOptionPane.ERROR_MESSAGE);
	    return false;
	}

	int ipz = simParam.getImgPerZ(); 
	if ( img.depth % ipz != 0 || (img.depth/ipz) % zSlices != 0 ) {
	    String err=String.format(
		"<html>Number of images should be a multiple of<br>"+
		"ang x pha (%d) x z (%d), but is %d. </html>",
		simParam.getImgPerZ(), zSlices, img.depth );
	    JOptionPane.showMessageDialog( baseframe,
	     err, "fairSIM error",
	     JOptionPane.ERROR_MESSAGE);
	    return false;
	}
    
	return true;
    }


    /** Update to the user if the stack could be a video */
    void updateIfVideoCanBeUsed() {


	ImageSelector.ImageInfo inf = imgBox.getSelectedItem(); 
	if ( inf == null ) {
	    videoPosSlider.setEnabled(false);
	    importImageButton.setEnabled(false);
	    maxTimePointsLabel.setText("no image");	    
	    return;
	}
	
	int zplSel = (int)zSliceVideoSpinner.getVal();  

	// check if with the currently selected number of images,
	// the image could be imported as video
	
	if ( inf.depth % ( zplSel * simParam.getImgPerZ() ) != 0 ) {
	    videoPosSlider.setEnabled(false);
	    importImageButton.setEnabled(false);
	    videoAutoUpdateMode.setEnabled(false);
	    
	    videoPosLabel.setText("t: n/a");
	    maxTimePointsLabel.setText("t-max: n/a");
	    maxTimePointsLabel.setToolTipText("Total number of images divided by "+
		"number of z-slices is no integer, "+
		"please select matching number of z-slices");
	} else {
	    importImageButton.setEnabled(true);
	    videoPosSlider.setEnabled(true);
	    videoAutoUpdateMode.setEnabled(true);

	    totalTimePoints = inf.depth / zplSel / simParam.getImgPerZ(); 
	    if ( videoPosSlider.getValue() >= totalTimePoints ) {
		videoPosSlider.setValue(0);
	    } else {
		videoPosLabel.setText( 
		    String.format(" t:  % 5d", videoStackPositionTime+1));
	    }
	    videoPosSlider.setMaximum( totalTimePoints );
	    maxTimePointsLabel.setText("t-max: "+ totalTimePoints);
	    maxTimePointsLabel.setToolTipText("Number of available time points "+
		"at selected z");

	}

    }

    /** switch timelapse mode on and off */
    void changeTimelapseMode( boolean onoff ) {

	isVideoStack = onoff;
	
	sucessfulImport = false;
	runAutoUpdate(0);
	
	videoPosSlider.setEnabled(isVideoStack);
	zSliceVideoSpinner.setEnabled(isVideoStack);
    
	if (isVideoStack) {
	    updateIfVideoCanBeUsed();
	} else {
	    importImageButton.setEnabled(true);
	    videoAutoUpdateMode.setEnabled(true);
	}

    }


    /** start auto-update */
    void runAutoUpdate( int mode ) {

	if ( autoUpdater != null ) {
	   autoUpdater.updateMode = 0;
	   //autoUpdater.interrupt();
	   try {
		autoUpdater.join();
	   } catch (InterruptedException ex) {
		Tool.trace("Thread join interruption, should not happen");
	   }
	   autoUpdater = null;
	   videoAutoUpdateMode.setBackground( null );	
	   if (simPreviewDisplay != null)
		simPreviewDisplay.drop();
	}
	
	// mode 0 stops the thread and 
	if (mode==0) return;

	// all other modes are written to the thread
	// which is started if not existing
	autoUpdater = new AutoUpdateThread(mode);
	Tool.trace("Starting auto update, mode "+updateModes[mode]);
	autoUpdater.start();
    }


    /** this runs a background thread updating the images on the fly */
    class AutoUpdateThread extends Thread {
	
	int updateMode;
	AutoUpdateThread(int initialMode) {
	    updateMode = initialMode;
	}

	@Override
	public void run() {

	    long lastUpdate = 0;
		    
	    int simWidth  = theFFTImages[0][0].vectorWidth()*2;
	    int simHeight = theFFTImages[0][0].vectorHeight()*2;

	    simPreviewDisplay = 
		idpFactory.create(simWidth, simHeight, "SIM preview");
	    
	    simPreviewDisplay.addImage( Vec2d.createReal(simWidth, simHeight), 
		"SIM preview Image");

	    while(true) {

		// stop the thread if the import was not successful
		if (!sucessfulImport) {
		    videoAutoUpdateMode.setSelectedIndex(0);
		    Tool.trace("AUTO-UPDATE: Stopping auto updater, "+
		    "as last import was unsuccessful");
		    return;
		}

		// stop the thread if the update mode is 0
	        if (updateMode == 0) {
		    Tool.trace("AUTO-UPDATE: Stopping auto updater, as requested");
		    return;
		}


		// see if any updates have to be performed
		boolean doUpdate = simParam.compareRuntimeTimestamp( lastUpdate );

		// wait a few hundred ms if there is nothing to update
		if (!doUpdate) {
		    //Tool.trace("AUTO-UPDATE: no change, nothing to update. ");
		    try {
			sleep(400);
		    } catch (InterruptedException ex) {}
		    continue;
		}

		// run the automatic update for the preview..
		videoAutoUpdateMode.setBackground( Color.YELLOW );	
		lastUpdate = System.currentTimeMillis();
	
		// update the input images
		Tool.trace("AUTO-UPDATE: updating wide-field view ");

		importImages(imgBox.getSelectedItem(), videoStackPositionZ, 
		    videoStackPositionTime, true ); 

		// update the correction factors
		if (prefactorAutoUpdateAngle.isSelected()) {
		    runEstimateAngleVariation(prefactorMethodBox.getSelectedIndex());
		}
		if (prefactorAutoUpdatePhase.isSelected()) {
		    runEstimatePhaseVariation(prefactorMethodBox.getSelectedIndex());
		}

		// update the parameter estimation
		if (updateMode>2) {

		    SimAlgorithm.estimateParameters( 
			simParam, theFFTImages, 
			fsGUI.parc.getFitBand(), 
			fsGUI.parc.getFitExclude(), 
			null, 0, null);

		}

		// update the SIM reconstruction
		if (updateMode>1) {

		    if (simParam.otf() != null ) {
		
			Tool.trace("AUTO-UPDATE: updating wide-field view ");

			Vec2d.Real simRecon = SimAlgorithm.runReconstruction( 
			    simParam, theFFTImages, 
			    null,  0, false, SimParam.CLIPSCALE.BOTH, null);
		   
			if (simPreviewDisplay != null) {
			    simPreviewDisplay.setImage( simRecon, 0, "SIM image");
			    simPreviewDisplay.display();
			}
		    
		    } else {
			Tool.trace("AUTO-UPDATE: No OTF / parameters set for SIM");
		    }
		}
		
		videoAutoUpdateMode.setBackground( Color.GREEN );	
	    
	    }
	}
    }





    /** set the currently selected image */
    void importImageDialog( final ImageSelector.ImageInfo img, 
	final int zSlices, final int tPosition ) {

	Tool.trace("zSlice: "+zSlices);

	// TODO: Image size checks!
	if (!checkImage( img, zSlices ))
	    return;
	importImageButton.setEnabled(false);

	totalTimePoints = img.depth / zSlices / simParam.getImgPerZ();

	pxlSet  = ( img.micronsPerPxl > 0);
	pxlSize = ( pxlSet )?(img.micronsPerPxl*1000):(80);
	
	if (( pxlSize < 20 ) ||( pxlSize > 500 )) {
	    JOptionPane.showMessageDialog(baseframe,
		"<html>Pixel size in file not within 20nm..500nm<br >"+
		"Size in file is: "+String.format("%5.2f</html>", pxlSize),
		"Ignoring pixel size",
		JOptionPane.WARNING_MESSAGE		
		);
	    pxlSet  = false;
	    pxlSize = 80;
	}
	
	// slice selector - only shown if the image has enough slices
	JPanel p1 = new JPanel();
	p1.setLayout(new BoxLayout(p1, BoxLayout.LINE_AXIS));
	p1.setBorder(BorderFactory.createTitledBorder("Slice") );

	sliceBox = new Tiles.LComboBox<Integer>("#",1);
	sliceBox.box.setPrototypeDisplayValue("123");
	
	if (img.depth/totalTimePoints > simParam.getImgPerZ() ) {
	    // only enable the box if there is more than one slice
	    Integer [] n = new Integer[ 
		img.depth/simParam.getImgPerZ()/totalTimePoints ];
	    for (int i=0; i<n.length; i++) {
		n[i] = i+1;
	    }
	    sliceBox.newElements( n.length/2, n );
	    sliceBox.box.setEnabled(true );
	} else {
	    sliceBox.box.setEnabled(false );
	}
 
	// button to display slice projection
	showSlicesButton = new JButton("show slices");
	showSlicesButton.setToolTipText("Display wide-field projection of slices");
	showSlicesButton.addActionListener( new ActionListener() {
	    public void actionPerformed(ActionEvent e) {
		showSlices(img, zSlices);
	    }
	});

	p1.add( Box.createHorizontalGlue());
	p1.add( sliceBox );
	p1.add( Box.createRigidArea( new Dimension(5,0)));
	p1.add( showSlicesButton );
	p1.add( Box.createHorizontalGlue());

	// pixel size selector
	JPanel p2 = new JPanel();
	p2.setLayout(new BoxLayout(p2, BoxLayout.PAGE_AXIS));
	p2.setBorder(BorderFactory.createTitledBorder("Pixel size") );

	final JLabel pxlSizeLabel = new JLabel("from image file");
	pxlSizeLabel.setForeground(Color.BLUE);
	pxlSizeLabel.setAlignmentX(Component.CENTER_ALIGNMENT);

	Tiles.LNSpinner pxlSizeSpinner = 
	    new Tiles.LNSpinner( "nm" , pxlSize, 10, 500, 0.1 ); 

	if ( !pxlSet ) {
	    pxlSizeLabel.setText("not set");
	    pxlSizeLabel.setForeground(Color.RED);
	} 
	
	pxlSizeSpinner.addNumberListener( new Tiles.NumberListener() {
	    public void number( double pxl, Tiles.LNSpinner ign ) {
		pxlSize = pxl;
		pxlSet = true;
		pxlSizeLabel.setText("set manually to");
		pxlSizeLabel.setForeground(Color.GREEN.darker());
	    }
	});
	
	p2.add( pxlSizeLabel );
	p2.add( Box.createRigidArea( new Dimension( 5,5)));
	p2.add( pxlSizeSpinner );

	// background subtraction
	bgrSpinner = new Tiles.LNSpinner( "val" , 50, 0, 5000, 5 ); 
	bgrBox = new Tiles.LComboBox<String>("Use?", "no", "yes");
	bgrBox.box.setToolTipText("Subtract a constant background value?");

	JPanel bgrPanel = new JPanel();
	bgrPanel.setLayout( new BoxLayout( bgrPanel, BoxLayout.LINE_AXIS));
	bgrPanel.setBorder(BorderFactory.createTitledBorder("Background substraction") );

	bgrPanel.add( bgrBox );
	bgrPanel.add( bgrSpinner );

	// dialog	
	final JDialog imageDialog = new JDialog(baseframe,
	    "Image import settings", false);

	// bugfix: re-enable the select button if the dialog is closed
	imageDialog.addWindowListener( new WindowAdapter() {
	    @Override
	    public void windowClosed(WindowEvent e) {
			importImageButton.setEnabled(true);
	    }
	});


	JButton ok = new JButton("Import");
	JButton cl = new JButton("Cancel");
	JPanel p0 = new JPanel();
	p0.add( ok );
	p0.add( cl );

	// perform the image import
	ok.addActionListener( new ActionListener() {
	    public void actionPerformed(ActionEvent e) {
		if (!pxlSet) {	
		    JOptionPane.showMessageDialog(baseframe,
			"No pixel size set",
			"Pixel size",
			JOptionPane.ERROR_MESSAGE
		    );
		    return;
		}
		if (running) {	
		    JOptionPane.showMessageDialog(baseframe,
			"Import still running",
			"Import active",
			JOptionPane.ERROR_MESSAGE
		    );
		    return;
		}
		imageDialog.dispose();
		(new SwingWorker<Object, Object> () {
		    @Override
		    public Object doInBackground() {
			running = true;	
			videoStackPositionZ = sliceBox.getSelectedItem()-1; 
			importImages(img, videoStackPositionZ, tPosition,true);
			return null;
		    }
		    @Override
		    protected void done() {
			running = false;
			importImageButton.setEnabled(true);
			sucessfulImport = true;
			simParam.signalRuntimeChange();
		    }
		}).execute();
		
	    }
	});

	cl.addActionListener( new ActionListener() {
	    public void actionPerformed(ActionEvent e) {
		importImageButton.setEnabled(true);
		imageDialog.dispose();
	    }
	});

	// pack and display dialog
	JPanel dialog = new JPanel();
	dialog.setLayout( new BoxLayout( dialog, BoxLayout.PAGE_AXIS));
	dialog.add( p1 );
	dialog.add( p2 );
	dialog.add( bgrPanel );
	dialog.add( p0 );
	
	imageDialog.add( dialog );
	imageDialog.pack();
	imageDialog.setVisible(true);

    }

    



    /** import a set of images to be reconstructed */ 
    void importImages( final ImageSelector.ImageInfo img, 
	int zPos, int tPos, boolean showOutput ) {
    
	// images per z-plane, z-planes
	int ipz  = simParam.getImgPerZ(); 
	int zpl  = img.depth / ipz / totalTimePoints ;

	theImages    = new Vec2d.Real[ simParam.nrDir() ][];
	theFFTImages = new Vec2d.Cplx[ simParam.nrDir() ][];

	ourState.setText("Importing images...");
	ourState.setForeground(Color.BLUE);
	ourState.repaint(20);

	// square images
	final boolean doResize =
	    ((img.width != img.height) | ( img.width % 32 != 0));
	int imgSize = next32( Math.max(img.width, img.height) );
	
	final ImageDisplay curSlice = (showOutput)?(idpFactory.create( 
	    imgSize, imgSize, "Slice to reconstruct")):(null);
    
	// setup dislay for raw data
	if (showOutput) {
	    // reset and close for size mismatched
	    if ( rawDataDisplay != null && 
		(	rawDataDisplay.width() != imgSize 
		    ||	rawDataDisplay.height() != imgSize
		    ||	rawDataDisplay.getCount() != ipz+1 )) {
		rawDataDisplay.drop();
		rawDataDisplay = null;
	    }
	    // create with right size
	    if (rawDataDisplay == null) {
		rawDataDisplay = idpFactory.create( imgSize, imgSize, 
		"SIM: raw input data");
		for (int i=0; i<ipz+1; i++)
		    rawDataDisplay.addImage( Vec2d.createReal(imgSize, imgSize), 
			"Placeholder");

	    }
	} else {
	    if (rawDataDisplay != null)
		rawDataDisplay.drop();
	    rawDataDisplay=null;
	}

	Vec2d.Real widefield = Vec2d.createReal( imgSize, imgSize );
	int imageCount=0;

	for (int d=0; d<simParam.nrDir(); d++) {
	    
	    theImages[d]    = new Vec2d.Real[ simParam.dir(d).nrPha() ];
	    theFFTImages[d] = new Vec2d.Cplx[ simParam.dir(d).nrPha() ];

	    for (int p=0; p<simParam.dir(d).nrPha(); p++) {
	    
		int pos = simParam.getImgSeq().calcPosWithTime( d,p,zPos, tPos,
		    simParam.nrDir(), simParam.nrPha(), zpl );

		imageCount++;

		// import and store images  // TODO: move the FFT to somewhere else
		Vec2d.Real curImg   = imgSelect.getImage(img,pos).duplicate();

		// if we have to do background subtraction
		if (bgrBox.box.getSelectedIndex()==1) {
		    double r = SimUtils.subtractBackground( curImg,
			bgrSpinner.getVal());
		}
		
		SimUtils.fadeBorderCos( curImg , 10);

		// if we have to rescale the images...
		if (doResize) {
		    Tool.trace("Resizing input to square "+imgSize+"x"+imgSize);
		    Tool.tell("Input resize to "+imgSize+"x"+imgSize);
		    theImages[d][p] = Vec2d.createReal( imgSize, imgSize);
		    theImages[d][p].paste(curImg,0,0, false);
		} else {
		    theImages[d][p] = curImg;
		}
		
		
		theFFTImages[d][p] = Vec2d.createCplx( theImages[d][p] );
		theFFTImages[d][p].copy( theImages[d][p] );
		Transforms.fft2d( theFFTImages[d][p], false );
		
		if (rawDataDisplay != null) {
		    
		    rawDataDisplay.setImage( theImages[d][p], imageCount, 
		    String.format( 
		    "Image p: %2d/%2d a: %2d/%2d z: %2d/%2d t: %2d",
		    p+1, simParam.dir(d).nrPha(),
		    d+1, simParam.nrDir(), zPos+1, zpl, tPos+1 ));
		    
		    widefield.add( theImages[d][p]);
		}
	    }
	}
    
	if (rawDataDisplay!=null) {
	    widefield.scal( 1.f/simParam.getImgPerZ() );
	    rawDataDisplay.setImage( widefield, 0, "Image: proj. widefield");
	    rawDataDisplay.display();
	}

	ourState.setText(String.format("%s (slice z %d, t %d plx %3.0fnm)", img.name, zPos+1, tPos+1, pxlSize)); 
	ourState.setForeground(Color.GREEN.darker());

	simParam.setPxlSize( theImages[0][0].vectorWidth(), pxlSize/1000. );
    }

    /** calcuate angle-by-angle variation */
    void runEstimateAngleVariation(int method) {
    
	if (!sucessfulImport || theImages==null) {
	    JOptionPane.showMessageDialog( baseframe,
	     "Please import an image first", "fairSIM error",
	     JOptionPane.ERROR_MESSAGE);
	    return;
	}
		
	Tool.trace("Running angle-by-angle variation estimation:");

	double [] vals = new double[ simParam.nrDir() ];
	for (int ang=0; ang < simParam.nrDir(); ang++) {

	    // sum up all the phases
	    Vec2d.Real sum = Vec2d.createReal( theImages[ang][0] );
	    sum.add( theImages[ang] );
	    
	    if (method==0) {
		vals[ang] = sum.avr()/simParam.dir(ang).nrPha();
		Tool.trace("angle "+ang+": average "+String.format("%7.4f",vals[ang]));
	    }
	    if (method==1) {
		vals[ang] = sum.median()/simParam.dir(ang).nrPha();
		Tool.trace("angle "+ang+": median "+String.format("%7.4f",vals[ang]));
	    }
	}


	// calculate and set the correction factors
	double avr = MTool.average( vals ) ;