package imseProc.proc.demod;

import java.io.File;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.nio.DoubleBuffer;
import jafama.FastMath;
import oneLiners.OneLiners;

import org.eclipse.swt.widgets.Composite;
import binaryMatrixFile.BinaryMatrixFile;
import descriptors.gmds.GMDSSignalDesc;
import otherSupport.PhaseUnwrap2D;
import signals.Signal;
import signals.gmds.GMDSSignal;
import edu.emory.mathcs.jtransforms.fft.DoubleFFT_2D;
import imseProc.core.ByteBufferImage;
import imseProc.core.ComplexImage;
import imseProc.core.IMSEProc;
import imseProc.core.ImagePipeController;
import imseProc.core.Img;
import imseProc.core.ImgProcPipe;
import imseProc.core.ImgSource;

/** The main demodulation code, takes ComplexImage FFT input images and produces
 *  real images with a variety of output modes.
 *  
 * @author oliford
 */
public class DSHDemod extends ImgProcPipe {
	
	/** Components from incoming FFT */
	public static final int nComps = 5;
	public static final int COMPONENT_00 = 0;
	public static final int COMPONENT_p0 = 1;
	public static final int COMPONENT_pp = 2;
	public static final int COMPONENT_pm = 3;
	public static final int COMPONENT_0p = 4;	
	public static final String componentNames[] = new String[]{
		"(0, 0)", 
		"(+, 0)", 
		"(+, +)", 
		"(+, -)",
		"(0, +)",
	};
	/** Component selection boxes from incoming FFT image, set by controller */ 
	private int selections[][] = null;
	
	public static final String outputModeNames[] = new String[]{   
			"Amp (0, 0)", "Phase (0, 0)",
			"Amp (+, 0)", "Phase (+, 0)",
			"Amp (+, +)", "Phase (+, +)",
			"Amp (+, -)", "Phase (+, -)",
			"Amp (0, +)", "Phase (0, +)",
			"Pol 2|+,+| / |+,0|",
			"Pol 2|+,-| / |+,0|",
			"Pol 4(+,+)(+,-) / (+,0)²", 
			"Phase 4(+,+)(+,-) / (+,0)²", 
			"Contrast [Σ[|(+,+)|,|(+,-)|,|(+,0)|] / (|(0,0)| - |(0,0)[0]|)"
		};
		
	/** Output modes */ 
	public static final int OUTPUT_MODE_AMP_00 = 0;
	public static final int OUTPUT_MODE_PHS_00 = 1;
	public static final int OUTPUT_MODE_AMP_p0 = 2;
	public static final int OUTPUT_MODE_PHS_p0 = 3;
	public static final int OUTPUT_MODE_AMP_pp = 4;
	public static final int OUTPUT_MODE_PHS_pp = 5;
	public static final int OUTPUT_MODE_AMP_pm = 6;
	public static final int OUTPUT_MODE_PHS_pm = 7;
	public static final int OUTPUT_MODE_AMP_0p = 8;
	public static final int OUTPUT_MODE_PHS_0p = 9;
	public static final int OUTPUT_MODE_POL_pp = 10;
	public static final int OUTPUT_MODE_POL_pm = 11;
	public static final int OUTPUT_MODE_POL_px = 12;
	public static final int OUTPUT_MODE_PHS_px = 13;
	public static final int OUTPUT_MODE_CNT_pp = 14;
	
	/** Primary Output mode, set by controller */ 
	private int outputMode = OUTPUT_MODE_POL_px;
	
	/** Input interlacing mode, and what to do with it */
	public static final int INTERLACE_MODE_OUTPUT_ALL = 0;
	public static final int INTERLACE_MODE_OUTPUT_EVEN = 1;
	public static final int INTERLACE_MODE_OUTPUT_ODD = 2;
	public static final int INTERLACE_MODE_PAIR_DIFF = 3;
	public static final int INTERLACE_MODE_PAIR_SUM = 4;
	public static final int INTERLACE_MODE_ELLIP_COMP_A = 5;
	public static final int INTERLACE_MODE_ELLIP_COMP_B = 6;
	
	public static final String interlaceModeOpts[] = new String[]{
		"Show all",
		"Show even",
		"Show odd",
		"Phase Diff",
		"Phase Avg",
		"Ellip Comp. A", //ellipticity compensation polarisation
		"Ellip Comp. B", //experimental ellipticity compensation polarisation
	};
	
	private int interlaceMode = INTERLACE_MODE_OUTPUT_ALL;
	
	private boolean unwrapPhases = false;
	
	/** Global intrinsic contrast compensation */
	private double adjustP0 = 1.0, adjustPP = 1.0, adjustPM = 1.0;
	
	/** Calibration image config */
	private String calibImgPath;
	private boolean calibPathValid = false; // set true when calibImgPath was read successfully 
	private int calibImgIdx = -1; //input index of calib image
	private double calibImgAng = 22.5 * Math.PI/180;
	private ByteBufferImage calibImage = null; //image of calib factor µ
	
	private double cmptCmlpxData0[][];
		
	/** Only for keeping track of selection saving */
	private int currentPulse = -1;
	
	public DSHDemod() { super(ByteBufferImage.class); }
	
	public DSHDemod(ImgSource source, int selectedIndex) {
		super(ComplexImage.class, source, selectedIndex);	
	}
	
	@Override
 	protected void preCalc(boolean settingsHadChanged){
		super.preCalc(settingsHadChanged);
		
		String experiment = IMSEProc.getMetaExp(connectedSource);
		int pulse = IMSEProc.getMetaPulse(connectedSource);
		
		//if there's no selections
		if(selections == null){
			//first, if we have a pulse number, see if there's any in the DB
			try{
				loadSelectionsFromGMDS(experiment, pulse);
				currentPulse = pulse;				
			}catch(Exception e){ 
				System.err.println("DSHDemod: Error loaded selecions from DB: " + e.getMessage());
			}
			
			if(selections != null){
				System.out.println("DSHDemod: Loaded selections from GMDS DB");
			}else{
				//try auto selections
				//if possible, do this from the our 'selected' image
				
				Img selImg = getSelectedImageFromConnectedSource();
			
				selections = ComponentsFinder.autoInitSelections((ComplexImage)selImg);
				updateAllControllers();
			}
		}
		
		if(pulse > 0 && currentPulse != pulse){
			//if the pulse has changed and there are saved selection in the new pulse, load those
			try{
				loadSelectionsFromGMDS(experiment, pulse);
				currentPulse = pulse;
			}catch(Exception e){ 
				System.err.println("Pulse changed, but couldn't load selection from new pulse. Keeping old. Err: " + e.getMessage());
			}
			currentPulse = pulse;
		}
				
		//we need to process the calibration image first
		
		Img calibSrcImg = connectedSource.getImage(calibImgIdx);		
		if(calibImgPath != null && calibImgPath.length() > 0){ //preferably from an external path
			loadExternalCalibImage();
			
		}else if(calibSrcImg != null){ //or from an image in this set 
			if(calibImage == null || calibImage.getWidth() != outWidth || calibImage.getHeight() != outHeight){
				if(calibImage != null)
					calibImage.destroy();

				calibImage = new ByteBufferImage(this, -1, outWidth, outHeight, ByteBufferImage.DEPTH_DOUBLE, false);
			}
			attemptCalc(calibImage, new Img[]{ calibSrcImg }, settingsHadChanged);
			
		}else{ //otherwise, not
			calibImage = null;
		}
		
		//if there's a 'time' sequence, apply the interlace selection
		if(interlaceMode != INTERLACE_MODE_OUTPUT_ALL){
			double time[] = (double[])connectedSource.getSeriesMetaData("/time");
			if(time != null){
				double newTime[] = new double[time.length/2];
				for(int i=0; i < newTime.length; i++){
					newTime[i] = (time[i*2] + time[i*2+1]) / 2;
				}
				setSeriesMetaData("/time", newTime);
			}
		}
	}
 	
 	/** Image allocation */
	protected boolean checkOutputSet(int nImagesIn){
		outWidth = inWidth;
		outHeight = inHeight;
	    int nImagesOut = (interlaceMode == INTERLACE_MODE_OUTPUT_ALL) ? nImagesIn : (nImagesIn/2);
		
        //allocate or re-use
        ByteBufferImage newOutSet[] = ByteBufferImage.checkBulkAllocation(this, (ByteBufferImage[])imagesOut, outWidth, outHeight, ByteBufferImage.DEPTH_DOUBLE, nImagesOut, ByteOrder.LITTLE_ENDIAN);
        if(newOutSet != imagesOut){
        	imagesOut = newOutSet;
        	return true;
        }
        
        return false;
	}
	
	protected int[] sourceIndices(int outIdx){		
		switch(interlaceMode){
			case INTERLACE_MODE_OUTPUT_ALL:
					return new int[]{ outIdx };				
			case INTERLACE_MODE_OUTPUT_EVEN:
					return new int[]{ 2*outIdx };
			case INTERLACE_MODE_OUTPUT_ODD:
				return new int[]{ 2*outIdx+1 };
			default:
				return new int[]{ 2*outIdx, 2*outIdx+1 };				
		}
	}

	/** Reusable temporary image for ellipticity info */
	private ByteBufferImage ellipImage = null;
	
	@Override
	public boolean doCalc(Img imageOutGen, Img sourceSet[], boolean settingsHadChanged) throws InterruptedException {
		ByteBufferImage imageOut = (ByteBufferImage)imageOutGen;
		
		//check all source images are ComplexImage  
		for(int j=0; j < sourceSet.length; j++){
			if(!(sourceSet[j] instanceof ComplexImage))
				throw new IllegalArgumentException("sourceSet["+j+"] is "+sourceSet[j].toString()+" not a ComplexImage");				
		}

		// Deal with the interlaced modes that use both input images first
		if(interlaceMode == INTERLACE_MODE_PAIR_DIFF || interlaceMode == INTERLACE_MODE_PAIR_SUM
				|| interlaceMode == INTERLACE_MODE_ELLIP_COMP_A || interlaceMode == INTERLACE_MODE_ELLIP_COMP_B){

			double evenCmptCmlpxData[][] = calcNeededComplexData((ComplexImage)sourceSet[0]);				
			double oddCmptCmlpxData[][] = calcNeededComplexData((ComplexImage)sourceSet[1]);

			//ellipticity compensation modes require normal mode calcs of both images first
			if(interlaceMode == INTERLACE_MODE_ELLIP_COMP_B){ 
				 calcEllipCompensatedThetaB(imageOut, evenCmptCmlpxData, oddCmptCmlpxData);

			}else if(interlaceMode == INTERLACE_MODE_ELLIP_COMP_A){
				if(ellipImage == null || ellipImage.getWidth() != outWidth || ellipImage.getHeight() != outHeight){
					if(calibImage != null)
						calibImage.destroy();

					ellipImage = new ByteBufferImage(this, -1, outWidth, outHeight, ByteBufferImage.DEPTH_DOUBLE, false);
				}
				
				ellipImage.startWriting(); try{	
					switch(outputMode){
						case OUTPUT_MODE_POL_pp: 
							calcTheta(ellipImage, evenCmptCmlpxData, COMPONENT_pp, COMPONENT_p0, false);
							calcTheta(imageOut, oddCmptCmlpxData, COMPONENT_pp, COMPONENT_p0, false);
							break;
						case OUTPUT_MODE_POL_pm: 
							calcTheta(ellipImage, evenCmptCmlpxData, COMPONENT_pm, COMPONENT_p0, false);
							calcTheta(imageOut, oddCmptCmlpxData, COMPONENT_pm, COMPONENT_p0, false); 
							break;	
						case OUTPUT_MODE_POL_px: 
							calcBalanced(ellipImage, evenCmptCmlpxData, false, false); 
							calcBalanced(imageOut, oddCmptCmlpxData, false, false); 
							break;
					}
				}finally{ ellipImage.endWriting(); }
	
				calcEllipCompensatedTheta(imageOut, ellipImage);		

			}else{ //phase sum/diff images
				switch(outputMode){
					case OUTPUT_MODE_PHS_00: calcPhaseDiffImage(imageOut, evenCmptCmlpxData, oddCmptCmlpxData, COMPONENT_00); break;
					case OUTPUT_MODE_PHS_p0: calcPhaseDiffImage(imageOut, evenCmptCmlpxData, oddCmptCmlpxData, COMPONENT_p0); break;
					case OUTPUT_MODE_PHS_pp: calcPhaseDiffImage(imageOut, evenCmptCmlpxData, oddCmptCmlpxData, COMPONENT_pp); break;
					case OUTPUT_MODE_PHS_pm: calcPhaseDiffImage(imageOut, evenCmptCmlpxData, oddCmptCmlpxData, COMPONENT_pm); break;	
					case OUTPUT_MODE_PHS_0p: calcPhaseDiffImage(imageOut, evenCmptCmlpxData, oddCmptCmlpxData, COMPONENT_0p); break;						
				default:
					throw new IllegalArgumentException("Interlaced pair sum/diff only works for phase output modes");
				}
			}

		}else{ //a normal 1-1 image to image mode (although it might still be interlaced input)

			double cmptCmlpxData[][] = calcNeededComplexData((ComplexImage)sourceSet[0]);

			if(outputMode == OUTPUT_MODE_CNT_pp && imageOut == imagesOut[0]){
				cmptCmlpxData0 = cmptCmlpxData;
			}

			switch(outputMode){
				case OUTPUT_MODE_POL_pp: calcTheta(imageOut, cmptCmlpxData, COMPONENT_pp, COMPONENT_p0, true); break;
				case OUTPUT_MODE_POL_pm: calcTheta(imageOut, cmptCmlpxData, COMPONENT_pm, COMPONENT_p0, true); break;	
				case OUTPUT_MODE_POL_px: calcBalanced(imageOut, cmptCmlpxData, false, imageOut != calibImage); break; //don't finish if producing calib image
				case OUTPUT_MODE_PHS_px: calcBalanced(imageOut, cmptCmlpxData, true,  imageOut != calibImage); break; //don't finish if producing calib image
				case OUTPUT_MODE_AMP_00: calcAmpImage(imageOut, cmptCmlpxData, COMPONENT_00); break;
				case OUTPUT_MODE_PHS_00: calcPhaseImage(imageOut, cmptCmlpxData, COMPONENT_00, imageOut != calibImage); break;
				case OUTPUT_MODE_AMP_p0: calcAmpImage(imageOut, cmptCmlpxData, COMPONENT_p0); break;
				case OUTPUT_MODE_PHS_p0: calcPhaseImage(imageOut, cmptCmlpxData, COMPONENT_p0, imageOut != calibImage); break;
				case OUTPUT_MODE_AMP_pp: calcAmpImage(imageOut, cmptCmlpxData, COMPONENT_pp); break;
				case OUTPUT_MODE_PHS_pp: calcPhaseImage(imageOut, cmptCmlpxData, COMPONENT_pp, imageOut != calibImage); break;
				case OUTPUT_MODE_AMP_pm: calcAmpImage(imageOut, cmptCmlpxData, COMPONENT_pm); break;
				case OUTPUT_MODE_PHS_pm: calcPhaseImage(imageOut, cmptCmlpxData, COMPONENT_pm, imageOut != calibImage); break;	
				case OUTPUT_MODE_AMP_0p: calcAmpImage(imageOut, cmptCmlpxData, COMPONENT_0p); break;
				case OUTPUT_MODE_PHS_0p: calcPhaseImage(imageOut, cmptCmlpxData, COMPONENT_0p, imageOut != calibImage); break;
				case OUTPUT_MODE_CNT_pp: calcContrastImage(imageOut, cmptCmlpxData, cmptCmlpxData0, COMPONENT_pp); break;
			}

		}
		
		return true;
	}

	private double[][] calcNeededComplexData(ComplexImage imgT){
		double cmptCmlpxData[][] = new double[nComps][];	

		if(outputMode == OUTPUT_MODE_AMP_00 || outputMode == OUTPUT_MODE_PHS_00
				|| outputMode == OUTPUT_MODE_CNT_pp || outputMode == OUTPUT_MODE_CNT_pp){
			cmptCmlpxData[COMPONENT_00] = selectComponent(COMPONENT_00, imgT, 1.0);
		}
		
		if(outputMode == OUTPUT_MODE_AMP_p0 || outputMode == OUTPUT_MODE_PHS_p0
				|| outputMode == OUTPUT_MODE_POL_pp || outputMode == OUTPUT_MODE_POL_pm 
				|| outputMode == OUTPUT_MODE_POL_px  || outputMode == OUTPUT_MODE_PHS_px|| outputMode == OUTPUT_MODE_CNT_pp){
			cmptCmlpxData[COMPONENT_p0] = selectComponent(COMPONENT_p0, imgT, adjustP0);
		}
		
		if(outputMode == OUTPUT_MODE_AMP_pp || outputMode == OUTPUT_MODE_PHS_pp 
				|| outputMode == OUTPUT_MODE_POL_pp || outputMode == OUTPUT_MODE_POL_px  || outputMode == OUTPUT_MODE_PHS_px
				|| outputMode == OUTPUT_MODE_CNT_pp || outputMode == OUTPUT_MODE_CNT_pp){
			cmptCmlpxData[COMPONENT_pp] = selectComponent(COMPONENT_pp, imgT, adjustPP);
		}
		
		if(outputMode == OUTPUT_MODE_AMP_pm || outputMode == OUTPUT_MODE_PHS_pm
				|| outputMode == OUTPUT_MODE_POL_pm || outputMode == OUTPUT_MODE_POL_px 
				|| outputMode == OUTPUT_MODE_PHS_px|| outputMode == OUTPUT_MODE_CNT_pp){
			cmptCmlpxData[COMPONENT_pm] = selectComponent(COMPONENT_pm, imgT, adjustPM);
		}
		
		if(outputMode == OUTPUT_MODE_AMP_0p || outputMode == OUTPUT_MODE_PHS_0p
				|| outputMode == OUTPUT_MODE_POL_px || outputMode == OUTPUT_MODE_PHS_px){
			cmptCmlpxData[COMPONENT_0p] = selectComponent(COMPONENT_0p, imgT, 1);
		}
		
		return cmptCmlpxData;
	}
	

	/** Selects the given component from the FT and inverse FTs it, returning the complex data array */ 
	private double[] selectComponent(int sel,  ComplexImage imageIn, double adjust) {
		if(imageIn == null || selections == null){
			return null;
		}
		
		int sx0 = selections[sel][0], sy0 = selections[sel][1];
		int sw = selections[sel][2], sh = selections[sel][3];
		int sxC = selections[sel][4], syC = selections[sel][5];
		int sx1 = sx0 + sw - 1, sy1 = sy0 + sh - 1;
		int iw = imageIn.getWidth(), ih = imageIn.getHeight();
		
		if(sx1 >= iw || sy1 >= ih
				//sw <= 0 || sh <= 0 || sxC < sx0 || sxC > sx1 || syC < sy0 || syC > sy1 || 
				){
			System.err.println("DSHDemod.demod("+sel+"): Invalid selection / centre freq.");
			return null;
		}
	
		try{
			
			//pr00,pi00,pr12,pi12...
			//and then rows go [+ve cols, -ve cols, +ve cols ....]
			//and then the block of rows pairs are [ +++++, ------] in y		
			double cmplxData[] = new double[iw*ih*2];
			
			DoubleFFT_2D fft = new DoubleFFT_2D(ih, iw);
			
			imageIn.startReading();
			try{
				for(int iy = sy0; iy <= sy1; iy++){
					for(int ix = sx0; ix <= sx1; ix++){
						//we need the img coords (ix,iy) in frequency coords (x,y)
						int x = ix - (iw/2);
						int y = iy - (ih/2);
						if(x < 0) x += iw; //and remeber that the -ve x are on the RHS 
						if(y < 0) y += ih; //and -ve y at the bottom
						
						cmplxData[2*y*iw + x*2] = adjust * imageIn.getReal(ix, iy);
						cmplxData[2*y*iw + x*2+1] = adjust * imageIn.getImag(ix, iy);
					}		
				}
			}finally{
				imageIn.endReading();
			}
						
			fft.complexInverse(cmplxData, false);
			
			return cmplxData;
			
		}catch(InterruptedException e){
			System.err.println("DSHDemod.demod("+sel+"): Interrupted while waiting for image lock.");
			return null;
		}
		
	}
	
	private void calcAmpImage(ByteBufferImage imageOut, double cmplxData[][], int sel){
		if(cmplxData[sel] == null) {
			imageOut.invalidate();
			return;
		}
		
		int iw = imageOut.getWidth();
		int ih = imageOut.getHeight();
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				double real = cmplxData[sel][y*iw*2 + x*2];
				double imag = cmplxData[sel][y*iw*2 + x*2+1];
				
				imageOut.setPixelValue(x, y, FastMath.sqrt(real*real + imag*imag)); 
			}
		}
		
		imageOut.endWriting(); //force ranges invalid		
		imageOut.imageChanged(false);
		
	}
	
	private void calcContrastImage(ByteBufferImage imageOut, double cmplxData[][], double cmplxData0[][], int sel){
		if(cmplxData[sel] == null){
			imageOut.invalidate();
			return;
		}
				
		int iw = outWidth;
		int ih = outHeight;
		for(int y=0; y < ih; y++){			
			for(int x=0; x < iw; x++){
				double realPP = cmplxData[COMPONENT_pp][y*iw*2 + x*2];
				double imagPP = cmplxData[COMPONENT_pp][y*iw*2 + x*2+1];
				double absPP = FastMath.sqrt(realPP*realPP + imagPP*imagPP);

				double realPM = cmplxData[COMPONENT_pm][y*iw*2 + x*2];
				double imagPM = cmplxData[COMPONENT_pm][y*iw*2 + x*2+1];
				double absPM = FastMath.sqrt(realPM*realPM + imagPM*imagPM);

				double realP0 = cmplxData[COMPONENT_p0][y*iw*2 + x*2];
				double imagP0 = cmplxData[COMPONENT_p0][y*iw*2 + x*2+1];
				double absP0 = FastMath.sqrt(realP0*realP0 + imagP0*imagP0);

				double real00 = cmplxData[COMPONENT_00][y*iw*2 + x*2];
				double imag00 = cmplxData[COMPONENT_00][y*iw*2 + x*2+1];
				double abs00 = FastMath.sqrt(real00*real00 + imag00*imag00);
				
				double real000 = cmplxData0[COMPONENT_00][y*iw*2 + x*2];
				double imag000 = cmplxData0[COMPONENT_00][y*iw*2 + x*2+1];
				double abs000 = FastMath.sqrt(real000*real000 + imag000*imag000);
				
				if(x == 217 && y  == 189 )
					System.out.println("rar");
				 
				imageOut.setPixelValue(x, y, (abs00 <= abs000) ? Double.NaN : ((absPP + absPM + absP0) / (abs00 - abs000))); 
			}
		}
		
		imageOut.endWriting(); //force ranges invalid		
		imageOut.imageChanged(false);
		
	}
	
	
	
	private void calcPhaseImage(ByteBufferImage imageOut, double cmplxData[][], int sel, boolean subtractCalib){
		if(cmplxData[sel] == null || selections == null){
			imageOut.invalidate();
			return;
		}
		
		int sx0 = selections[sel][0], sy0 = selections[sel][1];
		int sw = selections[sel][2], sh = selections[sel][3];
		int sxC = selections[sel][4], syC = selections[sel][5];
		int sx1 = sx0 + sw - 1, sy1 = sy0 + sh - 1;
		
		int iw = imageOut.getWidth();
		int ih = imageOut.getHeight();
		int fx0 = sxC - (iw/2);
		int fy0 = syC - (ih/2);
		
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				double phase = (double)x * fx0 / iw + (double)y * fy0 / ih;
				double realOsc = FastMath.cos(-2 * Math.PI * phase);
				double imagOsc = FastMath.sin(-2 * Math.PI * phase);
				double realImg = cmplxData[sel][y*iw*2 + x*2];
				double imagImg = cmplxData[sel][y*iw*2 + x*2+1];
				
				double realD = (realOsc * realImg) - (imagOsc * imagImg);
				double imagD = (realOsc * imagImg) + (imagOsc * realImg);
				
				imageOut.setPixelValue(x, y, FastMath.atan2(imagD, realD));
			}
		}
		
		if(unwrapPhases){
			ByteBuffer phaseBBuff = imageOut.getWritableBuffer();
			DoubleBuffer phaseBuff = phaseBBuff.asDoubleBuffer();
			
			double phaseData[] = new double[iw*ih];
			phaseBuff.rewind();
			phaseBuff.limit(phaseBuff.capacity());
			phaseBuff.get(phaseData);
			PhaseUnwrap2D.unwrap(phaseData, iw, ih, false, false);
			phaseBuff.clear();
			phaseBuff.put(phaseData);
		}
		
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				imageOut.setPixelValue(x, y, imageOut.getPixelValue(x, y) * 180 / Math.PI);
			}
		}
		
		if(subtractCalib && calibImage != null && 
				calibImage.getWidth() == iw && calibImage.getHeight() == ih){ //phase from first image
			for(int y=0; y < ih; y++){
				for(int x=0; x < iw; x++){
					double v =  imageOut.getPixelValue(x, y) - calibImage.getPixelValue(x, y);
					if(v >= 180) v -= 360;
					if(v < -180) v += 360;
					imageOut.setPixelValue(x, y, v);
				}
			}
			
			
		}
	}
	
	private void calcPhaseDiffImage(ByteBufferImage imageOut, double evenCmplxData[][], double oddCmplxData[][], int sel){
		if(evenCmplxData[sel] == null || oddCmplxData[sel] == null || selections == null){
			imageOut.invalidate();
			return;
		}
		
		//int sx0 = selections[sel][0], sy0 = selections[sel][1];
		//int sw = selections[sel][2], sh = selections[sel][3];
		//int sxC = selections[sel][4], syC = selections[sel][5];
		//int sx1 = sx0 + sw - 1, sy1 = sy0 + sh - 1;
		
		
		int iw = imageOut.getWidth();
		int ih = imageOut.getHeight();
		//int fx0 = sxC - (iw/2);
		//int fy0 = syC - (ih/2);
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				//double phase = (double)x * fx0 / iw + (double)y * fy0 / ih;
				
				double realOsc = 1;//FastMath.cos(-2 * Math.PI * phase);
				double imagOsc = 0;//FastMath.sin(-2 * Math.PI * phase);
				double realImgE = evenCmplxData[sel][y*iw*2 + x*2];
				double imagImgE = evenCmplxData[sel][y*iw*2 + x*2+1];
				
				double realDE = (realOsc * realImgE) - (imagOsc * imagImgE);
				double imagDE = (realOsc * imagImgE) + (imagOsc * realImgE);
				
				//double phaseE = FastMath.atan2(imagDE, realDE);
				
				double realImgO = oddCmplxData[sel][y*iw*2 + x*2];
				double imagImgO = oddCmplxData[sel][y*iw*2 + x*2+1];
				
				double realDO = (realOsc * realImgO) - (imagOsc * imagImgO);
				double imagDO = (realOsc * imagImgO) + (imagOsc * realImgO);
				
				//double phaseO = FastMath.atan2(imagDO, realDO);
				
				if(interlaceMode == INTERLACE_MODE_PAIR_SUM){
					//imageOut.setPixelValue(x, y, ((phaseO + phaseE) % (Math.PI)));
					
					imageOut.setPixelValue(x, y, FastMath.atan2((imagDE*realDO + realDE*imagDO), (realDE*realDO - imagDE*imagDO)));
					
					
				}else{ //INTERLACE_MODE_PAIR_DIFF
					
					imageOut.setPixelValue(x, y, FastMath.atan2((imagDE*realDO - realDE*imagDO), (realDE*realDO + imagDE*imagDO)));
					
					//imageOut.setPixelValue(x, y, ((phaseO - phaseE) % (Math.PI)));
					
				}
				
			}
		}
		
		if(unwrapPhases){
			ByteBuffer phaseBBuff = imageOut.getWritableBuffer();
			DoubleBuffer phaseBuff = phaseBBuff.asDoubleBuffer();
			
			double phaseData[] = new double[iw*ih*8];
			phaseBuff.rewind();
			phaseBuff.limit(phaseBuff.capacity());
			phaseBuff.get(phaseData);
			PhaseUnwrap2D.unwrap(phaseData, iw, ih, false, false);
			phaseBuff.clear();
			phaseBuff.put(phaseData);
		}
		
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				imageOut.setPixelValue(x, y, imageOut.getPixelValue(x, y) * 180 / Math.PI);
			}
		}
		
	}
	
	private void calcTheta(ByteBufferImage imageOut, double cmplxData[][], int selA, int selB, boolean finish) {
		if(cmplxData[selA] == null || cmplxData[selB] == null ){
			imageOut.invalidate();
			return;
		}
				
		int iw = imageOut.getWidth();
		int ih = imageOut.getHeight();
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				double real, imag;
				
				real = cmplxData[selA][y*iw*2 + x*2];
				imag = cmplxData[selA][y*iw*2 + x*2+1];
				double ampA = FastMath.sqrt(real*real + imag*imag);

				real = cmplxData[selB][y*iw*2 + x*2];
				imag = cmplxData[selB][y*iw*2 + x*2+1];
				double ampB = FastMath.sqrt(real*real + imag*imag);
				
				if(finish){
					double ang = FastMath.atan2(2*ampA, ampB)/2 * 180 / Math.PI;				
					imageOut.setPixelValue(x, y, ang);
				}else{
					imageOut.setPixelValue(x, y, 2*ampA / ampB);
				}
			}
		}
	
		if(finish){
			imageOut.endWriting(); //force ranges invalid
			imageOut.imageChanged(false);
		}

	}
	
	/** Recover polarisation from (+,+)*(+,-)/(+,0)^2 which seems to balance the Γ(x') spectrum integral */
	private void calcBalanced(ByteBufferImage imageOut, double cmplxData[][], boolean phase, boolean finish) {
		if(cmplxData[COMPONENT_pp] == null ||cmplxData[COMPONENT_pm] == null ||
				cmplxData[COMPONENT_p0] == null){
			imageOut.invalidate();
			return;
		}
		
		double tanSq2CalAng = FastMath.pow2(FastMath.tan(2 * calibImgAng));
		
		int iw = imageOut.getWidth();
		int ih = imageOut.getHeight();
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				double ppRe = cmplxData[COMPONENT_pp][y*iw*2 + x*2];
				double ppIm = cmplxData[COMPONENT_pp][y*iw*2 + x*2+1];
				
				double pmRe = cmplxData[COMPONENT_pm][y*iw*2 + x*2];
				double pmIm = cmplxData[COMPONENT_pm][y*iw*2 + x*2+1];
				
				double p0Re = cmplxData[COMPONENT_p0][y*iw*2 + x*2];
				double p0Im = cmplxData[COMPONENT_p0][y*iw*2 + x*2+1];
				
				// nom = (+,+)*(+,-)
				double nomRe = ppRe * pmRe - ppIm * pmIm;
				double nomIm = ppRe * pmIm + ppIm * pmRe;

				// denom = (+,0)*(+,0)
				double denomRe = p0Re * p0Re - p0Im * p0Im;
				double denomIm = p0Re * p0Im + p0Im * p0Re;
				
				//res = nom / denom
				double c2d2 = denomRe*denomRe + denomIm*denomIm;				
				double resRe = (nomRe*denomRe + nomIm*denomIm) / c2d2;
				double resIm = (nomIm*denomRe - nomRe*denomIm) / c2d2;
				
				if(phase){

					imageOut.setPixelValue(x, y, FastMath.atan2(resIm, resRe)*180/Math.PI);
				}else{
						
					double tanSq2θ = -4*resRe; 
					
					if(x == imageOut.getWidth()/2 && y == imageOut.getHeight()/2){
						tanSq2θ = 1 * tanSq2θ;
					}
							
									
					if(finish){
						resRe *= 4;
						
						if(calibImage != null){
							
							// calibImage = tan(2*mang)² = µ² tan(2*calibAngImg)²
							// µ² = calibImage / tan(2*calibAngImg)²
							
							// tanSq2θ = tan(2*mang)² = µ² tan(2*cang)
							// cang = atan(tanSq2θ / µ²)/2
							// cang = atan(tanSq2θ * tan(2*calibAngImg)² / calibImage)/2						
							tanSq2θ *= tanSq2CalAng / calibImage.getPixelValue(x, y);						
							
						}
							
						double ang = FastMath.atan(FastMath.sqrt(tanSq2θ))/2;
						
						ang *= 180/Math.PI;
						
						//if(Double.isInfinite(ang) || Double.isNaN(ang))
						//	ang = -2;
						
						imageOut.setPixelValue(x, y, ang);
					}else{
						imageOut.setPixelValue(x, y, tanSq2θ);
					}
				}
			}
		}
		
		if(finish){
			imageOut.endWriting(); //force ranges invalid
			imageOut.imageChanged(false);
		}

	}
	
	/** Recover polarisation from (+,+)E*(+,-)E/(+,0)O^2
	 * E = even(FLC on), O = odd(FLC off)
	 * That should give theta - without contamination from chi */
	private void calcEllipCompensatedThetaB(ByteBufferImage imageOut, double cmplxDataEven[][], double cmplxDataOdd[][]) {
		if(cmplxDataEven[COMPONENT_pp] == null ||cmplxDataEven[COMPONENT_pm] == null ||
				cmplxDataEven[COMPONENT_p0] == null || cmplxDataOdd[COMPONENT_p0] == null ||
				cmplxDataOdd[COMPONENT_pp] == null || cmplxDataOdd[COMPONENT_pm] == null){
			imageOut.invalidate();
			return;
		}
		
		int iw = imageOut.getWidth();
		int ih = imageOut.getHeight();
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				double ppEvenRe = cmplxDataEven[COMPONENT_pp][y*iw*2 + x*2];
				double ppEvenIm = cmplxDataEven[COMPONENT_pp][y*iw*2 + x*2+1];
				
				double pmEvenRe = cmplxDataEven[COMPONENT_pm][y*iw*2 + x*2];
				double pmEvenIm = cmplxDataEven[COMPONENT_pm][y*iw*2 + x*2+1];
				
				double p0OddRe = cmplxDataOdd[COMPONENT_p0][y*iw*2 + x*2];
				double p0OddIm = cmplxDataOdd[COMPONENT_p0][y*iw*2 + x*2+1];
				
				// nom = (+,+)e*(+,-)e
				double nomRe = ppEvenRe * pmEvenRe - ppEvenIm * pmEvenIm;
				double nomIm = ppEvenRe * pmEvenIm + ppEvenIm * pmEvenRe;

				// denom = (+,0)o*(+,0)o
				double denomRe = p0OddRe * p0OddRe - p0OddIm * p0OddIm;
				double denomIm = p0OddRe * p0OddIm + p0OddIm * p0OddRe;
				
				//res = nom / denom
				double c2d2 = denomRe*denomRe + denomIm*denomIm;				
				double resRe = (nomRe*denomRe + nomIm*denomIm) / c2d2;
				double resIm = (nomIm*denomRe - nomRe*denomIm) / c2d2;
				
				double ang = FastMath.atan(Math.sqrt(FastMath.abs(-resRe * 4 - 1)))/2;
				
				//if(Double.isInfinite(ang) || Double.isNaN(ang))
				//	ang = -2;
				ang *= 180/ Math.PI;
				
				imageOut.setPixelValue(x, y, ang);			
			}
		}
		
		imageOut.endWriting(); //force ranges invalid
		imageOut.imageChanged(false);

	}
	
	private void calcEllipCompensatedTheta(ByteBufferImage thetaImg, ByteBufferImage ellipImg) {
		//for the ADSH system, the raw thetaImg isnt't actually tan²(2θ), but this:
		// X = -( cos²(2θ) + tan²(χ) ) / sin²(2θ)
		
		//so we recover θ from:
		// sin²(2θ) = (tan²(χ) + 1) / (1 - X) 
		
		if(thetaImg == null || ellipImg == null){			
			return;
		}
		
		int iw = thetaImg.getWidth();
		int ih = thetaImg.getHeight();
		for(int y=0; y < ih; y++){
			for(int x=0; x < iw; x++){
				double X = thetaImg.getPixelValue(x, y);
				double tanSqχ = ellipImg.getPixelValue(x, y);
				tanSqχ = -1 / tanSqχ;
				double a = FastMath.atan(FastMath.sqrt(tanSqχ));
				//a = 30 *Math.PI/180;
				tanSqχ = FastMath.pow2(FastMath.tan(a));
				double sinSq2θ = (tanSqχ + 1) / (1 - X);
				
				
				double ang;
				ang = FastMath.asin(FastMath.sqrt(sinSq2θ))/2;
				//ang =  FastMath.atan(FastMath.sqrt(tanSqχ));
				//ang = FastMath.atan(1 / FastMath.sqrt(-X))/2;
				thetaImg.setPixelValue(x, y, ang*180/Math.PI);
			}
		}
		
		thetaImg.endWriting(); //force ranges invalid
		thetaImg.imageChanged(false);

	}
	
	public void setWindow(int selection, int x0, int y0, int width, int height) {
		if(selections == null)
			selections = new int[nComps][6];
		selections[selection][0] = x0;
		selections[selection][1] = y0;
		selections[selection][2] = width;
		selections[selection][3] = height;
		selections[selection][4] = x0 + width/2;
		selections[selection][5] = y0 + height/2;
		settingsChanged = true;
		updateAllControllers();
		if(autoCalc)
			calc();
	}

	public void setCentralFreq(int selection, int x, int y) { 
		if(selections == null)
			selections = new int[nComps][6];
		selections[selection][4] = x;
		selections[selection][5] = y;
		settingsChanged = true;
		updateAllControllers();
		if(autoCalc)
			calc();
	}

	@Override
	public ImagePipeController createPipeController(Class interfacingClass, Object args[], boolean asSink) {
		if(interfacingClass == Composite.class){
			ImagePipeController controller = new DSHDemodSinkSWTController(this, (Composite)args[0], (Integer)args[1], asSink);
			controllers.add(controller);
			return controller;
		}
		
		return null;
	}
	
	@Override
	public DSHDemod clone() { return new DSHDemod(connectedSource, getSelectedSourceIndex()); }
		
	public int[][] getSelections() { return selections;	}
	
	public void setPhaseUnwrap(boolean unwrapPhases) { 
		if(this.unwrapPhases == unwrapPhases)
			return;
		
		this.unwrapPhases = unwrapPhases;
		this.settingsChanged = true;
		updateAllControllers();
		if(autoCalc)
			calc();
	}
		
	public void setOutputMode(int outputMode){
		if(this.outputMode == outputMode)
			return;
		
		this.outputMode = outputMode;
		this.settingsChanged = true;
		updateAllControllers();
		if(autoCalc)
			calc();
	}
	public int getOutputMode(){ return this.outputMode; }
	
	public void setComponentAdjust(double adjustP0, double adjustPP, double adjustPM) {
		if(this.adjustP0 == adjustP0 &&
				this.adjustPP == adjustPP &&
				this.adjustPM == adjustPM)
			return;
		
		this.adjustP0 = adjustP0;
		this.adjustPP = adjustPP;
		this.adjustPM = adjustPM;
		this.settingsChanged = true;
		updateAllControllers();
		if(autoCalc)
			calc();
	}
	
	public void setInterlaceMode(int interlaceMode) {
		this.interlaceMode = interlaceMode;		
		updateAllControllers();
		if(autoCalc)
			calc();
	}
	

	public void setCalibImage(String calibImgPath, int calibImg, double calibImgAng) {
		if( ((this.calibImgPath == null && calibImgPath == null) || (this.calibImgPath != null && this.calibImgPath.equals(calibImgPath))) &&
				this.calibImgIdx == calibImg && this.calibImgAng == calibImgAng)
			return;
		this.calibImgPath = calibImgPath;
		this.calibImgAng = calibImgAng;		
		this.calibImgIdx = calibImg;
		this.calibPathValid = false;
		this.settingsChanged = true;
		updateAllControllers();
		if(autoCalc)
			calc();
	}
	
	private void loadSelectionsFromGMDS(){
		if(connectedSource == null) return;

		String experiment = IMSEProc.getMetaExp(connectedSource);
		int pulse = IMSEProc.getMetaPulse(connectedSource);
		
		loadSelectionsFromGMDS(experiment, pulse);
		currentPulse = pulse;
	}
	
	public void loadSelectionsFromGMDS(String experiment, int pulse){
		GMDSSignalDesc sigDesc = new GMDSSignalDesc(pulse, experiment, "DSHDemod/Selections");
		
		GMDSSignal sig = (GMDSSignal)IMSEProc.globalGMDS().getSig(sigDesc);
		
		selections = (int[][])sig.get2DData();
		
		updateAllControllers();
		if(autoCalc)
			calc();
	}
	
	public void saveSelectionsToGMDS(int pulse){
		if(connectedSource == null || selections == null) return;
		
		String experiment = IMSEProc.getMetaExp(connectedSource);
		if(pulse < 0)
			pulse = IMSEProc.getMetaPulse(connectedSource);
		
		GMDSSignalDesc sigDesc = new GMDSSignalDesc(pulse, experiment, "DSHDemod/Selections");
		
		GMDSSignal sig = new GMDSSignal(sigDesc, selections.clone());

		IMSEProc.globalGMDS().writeToCache(sig);
	}
	
	private String autoSelectSyncObject = "autoSync";
	public void findAutoSelections(){
		IMSEProc.ensureFinalUpdate(autoSelectSyncObject, new Runnable() {
			@Override
			public void run() {
				Img selImg = getSelectedImageFromConnectedSource();
				settingsChanged = true;
				selections = ComponentsFinder.autoInitSelections((ComplexImage)selImg);
				calc();
				updateAllControllers();
			}
		});
	}
	
	private void loadExternalCalibImage(){		
		double calibData[][] = null;
		
		String path = calibImgPath;
		if(calibImgPath.matches(".*@[0-9.]*")){
			this.calibImgAng = OneLiners.mustParseDouble(calibImgPath.replaceFirst(".*@([0-9.]*)", "$1")) * Math.PI/180;
			path = calibImgPath.replaceFirst("(.*)@[0-9.]*", "$1");			
		}
		
		
		//see if it's a file first
		if(new File(path).canRead()){
			//as a Signal file?
			try {
				Signal sig = new Signal(path);
				calibData = (double[][]) sig.getDataAsType(double.class);
				
			}catch (Exception e) { }
			
			if(calibData == null){
				//binary matrix file?
				try{
					calibData = BinaryMatrixFile.mustLoad(path, false);
				}catch (Exception e) { }
			}
			
			if(calibData == null){
				System.err.println("Calib image path '"+calibImgPath+"' seems to be a file, but isn't a valid DataSignals or BinaryMatrix file");
				calibImage = null;
				return;
			}
		}else{
			try{
				//GMDS signal path?
				GMDSSignalDesc desc = new GMDSSignalDesc(path);
				GMDSSignal sig = (GMDSSignal)IMSEProc.globalGMDS().getSig(desc);
				calibData = (double[][]) sig.getDataAsType(double.class);
			}catch(Exception e){
				System.err.println("Calib image path '"+calibImgPath+"' not a file or a GMDS signal path.");
				return;
			}
		}
		
		calibImage = new ByteBufferImage(this, -1, calibData[0].length, calibData.length, ByteBufferImage.DEPTH_DOUBLE, false);
		try {
				calibImage.startWriting(); try{
			
				//doCalc() uses calibImage as if it were part of the DSH processing
				//i.e., if it was tan²(2θ)
				for(int iY=0; iY < calibData.length; iY++){
					for(int iX=0; iX < calibData[0].length; iX++){
						calibImage.setPixelValue(iX,  iY,  
								FastMath.pow2(FastMath.tan(2*calibData[iY][iX]*Math.PI/180)) );				
					}
				}
				
			}finally{ calibImage.endWriting(); }
				
		} catch (InterruptedException e) {
			e.printStackTrace();
		}
		calibImage.imageChanged(false);
		calibPathValid = true;
	}
	
	public String getCalibPath(){ return calibPathValid ? calibImgPath : null; }

}