FocusingMotors.java 176 KB
Newer Older
1
package com.elphel.imagej.calibration.hardware;
Andrey Filippov's avatar
Andrey Filippov committed
2

3 4
import java.io.BufferedWriter;
import java.io.File;
Andrey Filippov's avatar
Andrey Filippov committed
5 6 7
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
8
import java.io.FileWriter;
Andrey Filippov's avatar
Andrey Filippov committed
9 10 11 12 13 14 15 16 17 18 19
import java.io.IOException;
import java.io.InputStream;
import java.io.OutputStream;
import java.net.MalformedURLException;
import java.util.ArrayList;
import java.util.Properties;

import javax.xml.parsers.DocumentBuilder;
import javax.xml.parsers.DocumentBuilderFactory;
import javax.xml.parsers.ParserConfigurationException;

20 21
import org.apache.commons.configuration.ConfigurationException;
import org.apache.commons.configuration.XMLConfiguration;
Andrey Filippov's avatar
Andrey Filippov committed
22 23 24
import org.w3c.dom.Document;
import org.xml.sax.SAXException;

25 26 27 28 29 30 31
import com.elphel.imagej.calibration.CalibrationFileManagement;
import com.elphel.imagej.calibration.FocusingField;
import com.elphel.imagej.calibration.LensAdjustment;
//import com.elphel.imagej.calibration.hardware.CalibrationHardwareInterface.FocusingMotors.FocusingHistory;
//import com.elphel.imagej.calibration.hardware.CalibrationHardwareInterface.FocusingMotors.FocusingSharpness;
//import com.elphel.imagej.calibration.hardware.CalibrationHardwareInterface.FocusingMotors.FocusingHistory.FocusingState;
//import com.elphel.imagej.calibration.hardware.CalibrationHardwareInterface.FocusingMotors.FocusingSharpness.FocusingSharpnessItem;
32 33
import com.elphel.imagej.common.PolynomialApproximation;
import com.elphel.imagej.common.WindowTools;
Andrey Filippov's avatar
Andrey Filippov committed
34

Andrey Filippov's avatar
Andrey Filippov committed
35 36
import Jama.LUDecomposition;
import Jama.Matrix;
37 38 39 40
import ij.IJ;
import ij.Prefs;
import ij.gui.GenericDialog;
import ij.text.TextWindow;
Andrey Filippov's avatar
Andrey Filippov committed
41

42
public class FocusingMotors{ // should it be static?
Andrey Filippov's avatar
Andrey Filippov committed
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

    public int motor3rd_signed_hysteresis=100;
//    public double [][] motorsRange={{-32000,32000},{-32000,32000},{-32000,32000}}; // hardware limited to 16 bits signed
    public double [][] motorsRange={{-32000,25000},{-32000,25000},{-32000,25000}}; // need longer spring, bad hysteresis at the center, starting at +12000
	public String ipAddress="192.168.0.236";
	public int stepsPerRevolution=3584; //112*32
	public double threadPitch=0.4; // mm/revolution
	public double linearReductionRatio=4.0/38.0; // sensor movement per 3 screws together movement
	public double waitPerStep=0.005; //~ 2 sec/revolution
	public double waitExtra=1.0;
	public double calmMotors=1.0; // extra time to wait to calm
	public int    positionTolerance=2; // actual position should not differ than this from the required
	public int [] lastMovedTo=null;
    public int[] curpos=new int[3];
    public boolean [] lastDirs= {false,false,false}; // true if motor was going in positive direction
    public FocusingHistory focusingHistory=new FocusingHistory(1000);
    public FocusingSharpness focusingSharpness=new FocusingSharpness(1000);
    public int debugLevel=1;
    public String motorsStatePath="lensAdjustmentMotorPosition.xml"; // will be saved in ~/.imagej/
62
    public String prefsDirectory=null; // prefs.getPrefsDir()+Prefs.getFileSeparator(); // "~/.imagej/" - too early !!!
Andrey Filippov's avatar
Andrey Filippov committed
63 64 65 66 67 68
    private Properties motorProperties=new Properties();
    public boolean interactiveRestore=true; // ask for confirmation to restore motor position
    public boolean stateValid=false;
    public boolean motorHardwareInitialized=false; // check if motors need to be initialized, enabled and reset (cold start of the board)
    public boolean motorHardwareInitializationInProgress=false;
    public boolean motorsEnabled=false; // disable motors when UV was turned on
69

Andrey Filippov's avatar
Andrey Filippov committed
70 71 72 73 74 75
    public double [][] getRange(){return this.motorsRange;}
    public int getStepsPerRevolution(){return this.stepsPerRevolution;}
    public double getThreadPitch(){return this.threadPitch;}
    public double getLinearReductionRatio(){return this.linearReductionRatio;}
    public double getMicronsPerStep(){return  1000.0*this.threadPitch*this.linearReductionRatio/stepsPerRevolution;}
    public void   setLinearReductionRatio(double lrr){this.linearReductionRatio=lrr;}
76

77 78 79 80 81 82
    public String getPrefsDir(){
    	if (prefsDirectory==null) prefsDirectory=Prefs.getPrefsDir();
    	if (prefsDirectory!=null) prefsDirectory+=Prefs.getFileSeparator(); // "~/.imagej/"
    	return prefsDirectory;
    }

Andrey Filippov's avatar
Andrey Filippov committed
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
    public boolean isInitialized(){
    	return this.stateValid;
    }
    public boolean isEnabled(){
    	return this.motorsEnabled;
    }
    public void setEnable(boolean enable){
    	this.motorsEnabled=enable;
    }
    /**
     *  Run after the camera was turned off while program is still running
     */
    public void resetInitialization(){
    	this.stateValid=false;
        this.motorHardwareInitialized=false; // check if motors need to be initialized, enabled and reset (cold start of the board)
        this.motorHardwareInitializationInProgress=false;
        this.motorsEnabled=false; // disable motors when UV was turned on
    }

102

Andrey Filippov's avatar
Andrey Filippov committed
103 104 105
    public void setCalmMotors (double calmMotors){
    	this.calmMotors=calmMotors;
    }
106 107


Andrey Filippov's avatar
Andrey Filippov committed
108 109 110 111 112 113 114 115 116 117 118 119 120
    public void setHysteresis (int hysteresis){
    	this.motor3rd_signed_hysteresis=hysteresis;
    }
    public void setDebug (int debugLevel){
    	this.debugLevel=debugLevel;
    	System.out.println("Motor debug level is set to "+this.debugLevel);
    }

    public int[] moveElphel10364Motors(boolean wait, int [] positions, double sleep, boolean showStatus, boolean hysteresis)  {
    	 return moveElphel10364Motors(wait, positions, sleep, showStatus, "", hysteresis);
    }

    public void saveMotorState() throws IOException{
121
    	if (!this.stateValid) return; // do not save until allowed
Andrey Filippov's avatar
Andrey Filippov committed
122
    	for (int i=0;i<curpos.length;i++) this.motorProperties.setProperty("motor"+(i+1),this.curpos[i]+"");
123
    	String path=this.getPrefsDir()+this.motorsStatePath;
Andrey Filippov's avatar
Andrey Filippov committed
124 125 126 127 128 129 130 131 132 133 134 135 136
    	OutputStream os;
    	try {
    		os = new FileOutputStream(path);
    	} catch (FileNotFoundException e1) {
    		String msg="Failed to open configuration file: "+path;
    		IJ.showMessage("Error",msg);
    		throw new FileNotFoundException (msg);
    	}
/*    	if (os==null) {
    		String msg="Failed to open configuration file for writing: "+path;
    		IJ.showMessage("Error",msg);
    		throw new IllegalArgumentException (msg);
    	}
137
*/
Andrey Filippov's avatar
Andrey Filippov committed
138 139 140 141 142 143
    	try {
    		this.motorProperties.storeToXML(os,
    				"last updated " + new java.util.Date(), "UTF8");
    	} catch (IOException e) {
    		String msg="Failed to write XML motor state file: "+path;
    		IJ.showMessage("Error",msg);
144
    		extracted(msg);
Andrey Filippov's avatar
Andrey Filippov committed
145 146 147 148 149 150 151 152
    	}
    	try {
    		os.close();
    	} catch (IOException e) {
    		// TODO Auto-generated catch block
    		e.printStackTrace();
    	}
    }
153

Andrey Filippov's avatar
Andrey Filippov committed
154
    public int [] loadMotorState() throws IOException{
155
    	String path=this.getPrefsDir()+this.motorsStatePath;
Andrey Filippov's avatar
Andrey Filippov committed
156 157 158 159 160 161 162 163 164 165 166 167 168 169
    	InputStream is;
    	try {
    		is = new FileInputStream(path);
    	} catch (FileNotFoundException e) {
    		String msg="Failed to open configuration file: "+path;
    		System.out.println("Warning: "+msg);
    		return null; // file does not exist, use current motor position (if it is initialized)
    	}
    	try {
    		motorProperties.loadFromXML(is);

    	} catch (IOException e) {
    		String msg="Failed to read XML configuration file: "+path;
    		IJ.showMessage("Error",msg);
170
    		extracted(msg);
Andrey Filippov's avatar
Andrey Filippov committed
171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
    	}
    	try {
    		is.close();
    	} catch (IOException e) {
    		// TODO Auto-generated catch block
    		e.printStackTrace();
    	}
		int [] savedPosition=new int [this.curpos.length];
    	for (int i=0;i<curpos.length;i++){
    		String sValue=	this.motorProperties.getProperty("motor"+(i+1));
    		if (sValue!=null) {
    			savedPosition[i]=Integer.parseInt(sValue);
//    			if (this.debugLevel>2) System.out.println("Read savedPosition["+i+"]="+savedPosition[i]);
    			System.out.println("Read savedPosition["+i+"]="+savedPosition[i]);
    		} else {
    			String msg="motor"+(i+1)+" is undefined in "+path+". If the file is corrupted you may delete it.";
    			IJ.showMessage("Error",msg);
188
    			extracted(msg);
Andrey Filippov's avatar
Andrey Filippov committed
189 190 191 192
    		}
    	}
    	return savedPosition;
    }
193 194 195
	private void extracted(String msg) throws IOException {
		throw new IOException (msg);
	}
Andrey Filippov's avatar
Andrey Filippov committed
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
    public void checkEnabled(){
    	if (isEnabled()) return;
    	GenericDialog gd = new GenericDialog("Enable motors");
		gd.addMessage("Motors are currently disabled (i.e. SFE is being glued). Do you want to enable them?");
	    gd.showDialog();
	    if (gd.wasCanceled()) {
	    	throw new RuntimeException("Command canceled");
	    }
	    setEnable(true);
    }
    /**
     *  Check if the motor position is synchronized with disk, open dialog if needed (or throw) (use boolean flag - dialog or throw)
     */
    public void checkSynchronized(){
     	if (this.stateValid) return; // nothing to do
    	// first check if the hardware is initialized, if not - do that.
212
    	if (motorHardwareInitializationInProgress) return;
Andrey Filippov's avatar
Andrey Filippov committed
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
    	motorHardwareInitializationInProgress=true;
		boolean wasEnabled= true; // innocent until proven guilty
    	if (!motorHardwareInitialized){
    		readElphel10364Motors(); // will update this.curpos
    		int [] testPosition=this.curpos.clone();
    		testPosition[0]+=2*this.positionTolerance;
    		double savedCalmMotors=this.calmMotors;
    		this.calmMotors=2;
        	moveElphel10364Motors(
        			true, // will not wait forever - will time out
        			testPosition, //int [] positions,
        			0.0, // double sleep,
        			true, // boolean showStatus,
        			"Looking if motors are initialized - trying to move motor1", // String message,
        			false); // , boolean hysteresis)
    		wasEnabled= maxPosDiff(testPosition)<=this.positionTolerance;
        	if (!wasEnabled){ // seems not to be enabled, and supposedly not initialized
        		String msg="Initializing motor driver";
    			if (this.debugLevel>0) System.out.println(msg);
    			IJ.showStatus("Initializing motor driver");
            	commandElphel10364Motors("http://"+this.ipAddress+"/10364.php?init&enable");
            	moveElphel10364Motors(
            			true, //will timeout if noit initialized
            			testPosition, //int [] positions,
            			0.0, // double sleep,
            			true, // boolean showStatus,
            			"Trying to move motor1", // String message,
            			false); // , boolean hysteresis)
            	if (maxPosDiff(testPosition)>this.positionTolerance){//  failed to initialize
            		testPosition[0]-=2*this.positionTolerance; // restore back this.curpos may be wrong (probably not)
                	this.calmMotors=savedCalmMotors;
                	moveElphel10364Motors(
                			true, //boolean wait (Will be forever, if not initialized)
                			testPosition, //int [] positions,
                			0.0, // double sleep,
                			true, // boolean showStatus,
                			"Trying to move motor1", // String message,
                			false); // , boolean hysteresis)
                	motorHardwareInitializationInProgress=false;
            		String msg1="Failed to initialize motor driver board";
        			IJ.showMessage("Error",msg);
	    	    	throw new RuntimeException(msg1);
            	}
        	}
    		testPosition[0]-=2*this.positionTolerance; // restore back this.curpos may be wrong (probably not)
        	this.calmMotors=savedCalmMotors;
        	moveElphel10364Motors(
        			true, //boolean wait (Will be forever, if not initialized)
        			testPosition, //int [] positions,
        			0.0, // double sleep,
        			true, // boolean showStatus,
        			"Trying to move motor1", // String message,
        			false); // , boolean hysteresis)
        	motorHardwareInitializationInProgress=false;
    		motorHardwareInitialized=true;
    	}
269

Andrey Filippov's avatar
Andrey Filippov committed
270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313
    	int [] savedPosition=null;
   // 	if (this.stateValid) return; //
    	if (!this.stateValid) try {
    		savedPosition=loadMotorState();
    	} catch (IOException e) {
    		// IJ.showMessage is already run
    		e.printStackTrace();
    	}
    	// is it the same position
    	GenericDialog gd;
    	if (this.debugLevel>2) System.out.println("savedPosition="+((savedPosition==null)?"null":("{"+savedPosition[0]+","+savedPosition[1]+","+savedPosition[2]+"}"))+
    			" motorHardwareInitialized="+motorHardwareInitialized+ " motorHardwareInitializationInProgress="+motorHardwareInitializationInProgress);
    	if (savedPosition!=null){
    		if (maxPosDiff(savedPosition)>this.positionTolerance){
    			if (this.interactiveRestore){
    	    		gd = new GenericDialog("Select motor position restoration");
    	    		gd.addMessage("Motor position saved at the previous run of this program {"+savedPosition[0]+":"+savedPosition[1]+":"+savedPosition[2]+"]");
    	    		gd.addMessage("differs from the curent position reported by the hardware: {"+this.curpos[0]+":"+this.curpos[1]+":"+this.curpos[2]+"]");
    	    		gd.addMessage("");
    	    		gd.addMessage("You may restore the saved position, keep new coordinates (and overwrite saved data) or abort current command if you select \"Cancel\"");
    	    		gd.addMessage("");
    	    		gd.addMessage("If the camera was powered down between the program runs, as it "+(wasEnabled?"does not seem to be":"seems to be")+",");
    	    		gd.addMessage("\"Restore\" should be used, if it was not - powered down and the motors were just moved by another application,");
    	    		gd.addMessage("\"Use current\" is appropriate. Of course, it could be powered down and moved - in that case it would be difficult to restore");
    	    		gd.enableYesNoCancel("Restore", "Use current");
    	    	    gd.showDialog();
    	    	    if (gd.wasCanceled()) {
    	    	    	throw new RuntimeException("Command canceled");
    	    	    }
    	    	    if (gd.wasOKed()) {
        	        	motorHardwareInitializationInProgress=true;
    	    	    	restorePosition(savedPosition); // otherwise just keep current
    	    	    }
    			} else { // !interactiveRestore
    	        	motorHardwareInitializationInProgress=true;
    				restorePosition(savedPosition);
    			}
    		} else {
    			if (this.debugLevel>0) System.out.println("Motors did not move since last time this program ran");
    		}
    	} else { // no saved position
        	if (this.debugLevel>2) System.out.println("--savedPosition="+((savedPosition==null)?"null":("{"+savedPosition[0]+","+savedPosition[1]+","+savedPosition[2]+"}")));
			if (this.interactiveRestore){
	    		gd = new GenericDialog("Missing motor position file");
314
	    		gd.addMessage("No motor position file ("+this.getPrefsDir()+this.motorsStatePath+") found.");
Andrey Filippov's avatar
Andrey Filippov committed
315 316 317 318 319 320
	    		gd.addMessage("If you run the program for the first time it is normal, you may click \"OK\" and the new file will be created");
	    		gd.addMessage("If you select \"Cancel\" the current command will be aborted, you may check the file and restart the program.");
	    	    gd.showDialog();
	    	    if (gd.wasCanceled()) {
	    	    	throw new RuntimeException("Command canceled");
	    	    }
321
			}
Andrey Filippov's avatar
Andrey Filippov committed
322 323 324 325 326 327 328 329
    	}
    	motorHardwareInitializationInProgress=false;
    	this.stateValid=true;
    }
    private void restorePosition(int [] position){
//    	int [] inversePosition=new int [position.length];
//    	for (int i=0;i<inversePosition.length;i++) inversePosition[i]=-position[i];
    	if (this.debugLevel>2) System.out.println("restorePosition {"+position[0]+","+position[1]+","+position[2]+"}, motorHardwareInitializationInProgress="+motorHardwareInitializationInProgress);
330

Andrey Filippov's avatar
Andrey Filippov committed
331
    	commandElphel10364Motors("http://"+this.ipAddress+"/10364.php?m1="+position[0]+"&m2="+position[1]+"&m3="+position[2]+"&reset");
332 333

/*
Andrey Filippov's avatar
Andrey Filippov committed
334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354
 *      Code for old FPGA firmware
    	int [] wasCurPos=this.curpos.clone();
    	moveElphel10364Motors( // checks again and opens dialog
    			true, //boolean wait,
    			inversePosition, //int [] positions,
    			0.0, // double sleep,
    			true, // boolean showStatus,
    			"Restoring inversed saved position", // String message,
    			false); // , boolean hysteresis)
    	if (this.debugLevel>0) System.out.println("*** RESETING MOTORS at position {"+this.curpos[0]+","+this.curpos[1]+","+this.curpos[2]+"}");
    	commandElphel10364Motors("http://"+this.ipAddress+"/10364.php?reset"); // that would zero hardware motor registers
    	// Move back to the original location
    	moveElphel10364Motors(
    			true, //boolean wait,
    			wasCurPos, //int [] positions,
    			0.0, // double sleep,
    			true, // boolean showStatus,
    			"Restoring saved position", // String message,
    			true); // , boolean hysteresis)
    			*/
	}
355

Andrey Filippov's avatar
Andrey Filippov committed
356
    //if (maxPosDiff(positions)>this.positionTolerance){
357

Andrey Filippov's avatar
Andrey Filippov committed
358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
    public int[] moveElphel10364Motors(boolean wait, int [] positions, double sleep, boolean showStatus, String message, boolean hysteresis)  {
    	// global this.motor3rd_signed_hysteresis
    	int i;
    	if (this.debugLevel>2) System.out.println("moveElphel10364Motors(), before checks: motorHardwareInitialized="+motorHardwareInitialized+
    			" motorHardwareInitializationInProgress="+motorHardwareInitializationInProgress);

    	checkEnabled(); // will either enable or throw;
    	if (this.debugLevel>2) System.out.println("moveElphel10364Motors(), after checkEnabled(): motorHardwareInitialized="+motorHardwareInitialized+
    			" motorHardwareInitializationInProgress="+motorHardwareInitializationInProgress);
    	checkSynchronized();
    	if (this.debugLevel>2) System.out.println("moveElphel10364Motors(), after checkSynchronized(): motorHardwareInitialized="+motorHardwareInitialized+
    			" motorHardwareInitializationInProgress="+motorHardwareInitializationInProgress);
    	readElphel10364Motors(); // sets this.curpos
    	if (hysteresis) {
    		int [] hystPositions=positions.clone();
    		boolean hystRequired=false;
    		for (i=0;i<3;i++) if ((positions[i]-this.curpos[i])*this.motor3rd_signed_hysteresis<0){
    			hystPositions[i]-=this.motor3rd_signed_hysteresis;
    			hystRequired=true;
    		}
    		if (hystRequired) {
    			moveElphel10364Motors(true, hystPositions, 0, showStatus, message+((message.length()>0)?": c":"C")+"ompensating motors hysteresis -");
    		}
    	}
    	if (this.debugLevel>2) System.out.println("moveElphel10364Motors(), before actual move: motorHardwareInitialized="+motorHardwareInitialized+
    			" motorHardwareInitializationInProgress="+motorHardwareInitializationInProgress);
    	return moveElphel10364Motors(wait, positions, sleep, showStatus, message);
    }

    private int[] moveElphel10364Motors(boolean wait, int [] positions, double sleep, boolean showStatus, String message)  {
    	IJ.showStatus(message+"Moving m1="+positions[0]+" m2="+positions[1]+" m3="+positions[2]);
    	if (this.debugLevel>0) System.out.println(message+((message.length()>0)?"":"Moving motors")+" m1="+positions[0]+", m2="+positions[1]+", m3="+positions[2]+" wait="+wait);
    	return moveElphel10364Motors(wait, positions, sleep);
    }

    private int[] moveElphel10364Motors(boolean wait, int [] positions, double sleep)  {
    	moveElphel10364Motors(wait, positions);
    	return readElphel10364Motors(sleep);
    }

    private int[] moveElphel10364Motors(boolean wait, int [] positions)  {
    	int extraCalm=4; // sec
    	if (OK2moveElphel10364Motors(positions) || motorHardwareInitializationInProgress) { // OK to override ~OK
    		this.lastMovedTo=positions.clone();
    		String uri="http://"+this.ipAddress+"/10364.php?m1="+positions[0]+"&m2="+positions[1]+"&m3="+positions[2];
    		if (wait) {
    			double travel=0.0;

    			for (int i=0;i<positions.length;i++){
    				if (positions[i]!=this.curpos[i]) this.lastDirs[i]=(positions[i]>this.curpos[i]);
    				double t=Math.abs(positions[i]-this.curpos[i]);
    				if (travel<t) travel=t;
    			}
    			int waitSeconds=(int) Math.ceil(travel*this.waitPerStep+this.waitExtra);
    			uri+="&wait="+waitSeconds;
    			commandElphel10364Motors(uri+"&calm="+this.calmMotors);
    			if (maxPosDiff(positions)>this.positionTolerance){
    				System.out.println("Motors ["+this.curpos[0]+","+this.curpos[1]+","+this.curpos[2]+"] did not reach destination "+
    						"["+positions[0]+","+positions[1]+","+positions[2]+"], travel was="+travel+", waitSeconds="+waitSeconds+
    						". Waiting "+extraCalm+" seconds more...");
    				commandElphel10364Motors(uri+"&calm="+extraCalm);
    				if (maxPosDiff(positions)>this.positionTolerance){
    					System.out.println("Motors ["+this.curpos[0]+","+this.curpos[1]+","+this.curpos[2]+"] did not reach destination "+
    							"["+positions[0]+","+positions[1]+","+positions[2]+"] after extra waiting");
    				}
    			}
    		} else {
    			commandElphel10364Motors(uri);
    			if (this.debugLevel>0) System.out.println("Warning: send move command with no wait ["+positions[0]+","+positions[1]+","+positions[2]+"]");

    		}
429
// save position after each movement
Andrey Filippov's avatar
Andrey Filippov committed
430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465
    		try {
				saveMotorState();
			} catch (IOException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}
    		return this.curpos;
    	} else {
			if (this.debugLevel>0) System.out.println("Warning: OK2moveElphel10364Motors("+positions[0]+","+positions[1]+","+positions[2]+") returned FALSE!");
    	}
    	return readElphel10364Motors(); // just read
    }
    private int maxPosDiff(int [] positions){
    	int max=0;
    	for (int i=0;i< positions.length;i++){
    		if ((positions[i]-this.curpos[i])>max) max= positions[i]-this.curpos[i];
    		if ((this.curpos[i]-positions[i])>max) max= this.curpos[i]-positions[i];
    	}
    	return max;
    }

    private boolean OK2moveElphel10364Motors( int [] positions)  {
    	// public double [][]    this.motorsRange{{-5000,5000},{-5000,5000},{-5000,5000}};
    	return (positions[0]>=this.motorsRange[0][0]) && (positions[0]<= this.motorsRange[0][1]) &&
    	(positions[1]>=this.motorsRange[1][0]) && (positions[1]<= this.motorsRange[1][1]) &&
    	(positions[2]>=this.motorsRange[2][0]) && (positions[2]<= this.motorsRange[2][1]);
    }


    public int[] readElphel10364Motors()  {
    	return readElphel10364Motors(0.0);
    }

    public int[] readElphel10364Motors(double sleep)  {
//    	if (isEnabled)
//   	if (isInitialized)
466
//    	checkSynchronized(); // Do not check synchronization here - just read current position (OK to do even if motorst are disabled)
Andrey Filippov's avatar
Andrey Filippov committed
467 468 469 470 471 472 473 474 475 476 477
    	return commandElphel10364Motors("http://"+this.ipAddress+"/10364.php?sleep="+sleep);
    }

    private int[] commandElphel10364Motors(String url)  {
    	Document dom=null;
    	if (this.debugLevel>2)	System.out.println("commandElphel10364Motors("+url+")");
    	try {
    		DocumentBuilderFactory dbf = DocumentBuilderFactory.newInstance();
    		DocumentBuilder db = dbf.newDocumentBuilder();
    		dom = db.parse(url);
    		if (!dom.getDocumentElement().getNodeName().equals("motors")) {
478
    			IJ.showMessage("Root element: expected 'motors', got'" + dom.getDocumentElement().getNodeName()+"'");
Andrey Filippov's avatar
Andrey Filippov committed
479 480
    			return null;
    		}
481 482 483
    		this.curpos[0]=Integer.parseInt((dom.getDocumentElement().getElementsByTagName("motor1").item(0).getChildNodes().item(0)).getNodeValue());
    		this.curpos[1]=Integer.parseInt((dom.getDocumentElement().getElementsByTagName("motor2").item(0).getChildNodes().item(0)).getNodeValue());
    		this.curpos[2]=Integer.parseInt((dom.getDocumentElement().getElementsByTagName("motor3").item(0).getChildNodes().item(0)).getNodeValue());
Andrey Filippov's avatar
Andrey Filippov committed
484 485 486 487 488 489 490
    	} catch(MalformedURLException e){
    		System.out.println("Please check the URL:" + e.toString() );
    		return null;
    	} catch(IOException  e1){
    		IJ.showStatus("");
    		String error = e1.getMessage();
    		if (error==null || error.equals(""))  error = ""+e1;
491 492 493
    		System.out.println("commandElphel10364Motors ERRROR: "+error);
    		// *********** Temporary removed message box (usually "HOST UNREACHABLE") *************
//    		IJ.showMessage("commandElphel10364Motors ERRROR", ""+error);
Andrey Filippov's avatar
Andrey Filippov committed
494 495 496 497 498
    		return null;
    	}catch(ParserConfigurationException pce) {
    		pce.printStackTrace();
    		return null;
    	}catch(SAXException se) {
499
    		se.printStackTrace();
Andrey Filippov's avatar
Andrey Filippov committed
500 501
    		return null;
    	}
502
    	// prevent accumulating position errors
Andrey Filippov's avatar
Andrey Filippov committed
503 504
    	if ((this.lastMovedTo!=null) && (Math.abs(this.lastMovedTo[0]-this.curpos[0])<=this.positionTolerance) &&
    			(Math.abs(this.lastMovedTo[1]-this.curpos[1])<=this.positionTolerance) &&
505
    			(Math.abs(this.lastMovedTo[2]-this.curpos[2])<=this.positionTolerance)) this.curpos=this.lastMovedTo.clone();
Andrey Filippov's avatar
Andrey Filippov committed
506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525
    	//    	    	  this.lastMovedTo=positions.clone();
    	//	public int    positionTolerance=2; // actual position should not differ than this from the required

    	return this.curpos;
    }

    public void clearHistory(){
    	this.focusingHistory.clear();
    }

    public void addToHistory( String sTimestamp,double temperature, double[][] psfMetrics, double [][][][] fullResults){ // null OK
    	this.focusingHistory.add(
    			sTimestamp,
    			temperature,
    			this.curpos,
    			this.lastDirs,
    			psfMetrics,
    			fullResults
    	);
    }
526

527 528 529 530 531 532
	public FocusingField.FocusingFieldMeasurement getThisFFMeasurement(
			FocusingField focusingField,
			String sTimestamp,
			double temperature,
			double[][] psfMetrics,
			double [][][][] fullResults){
533

534 535 536 537 538
		return focusingField.getFocusingFieldMeasurement(
				sTimestamp,   //focusingState.getTimestamp(),
				temperature, //focusingState.getTemperature(),
				this.curpos,      //focusingState.motorsPos,
				fullResults); //focusingState.getSamples());
539 540 541 542
	}



Andrey Filippov's avatar
Andrey Filippov committed
543
/*
544
 * FocusingHistory.FocusingState focusingState
Andrey Filippov's avatar
Andrey Filippov committed
545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
    public void addToHistory( String sTimestamp,double temperature, double[][] psfMetrics){ // null OK
    	this.focusingHistory.add(
    			sTimestamp,
    			temperature,
    			this.curpos,
    			this.lastDirs,
    			psfMetrics,
    			null
    	);
    }
*/
    public void addToHistory(String sTimestamp,double temperature){
    	this.focusingHistory.add(
    			sTimestamp,
    			temperature,
    			this.curpos,
    			this.lastDirs,
    			null,
    			null);
    }

    public void setLastMetrics(double[][] psfMetrics){
    	this.focusingHistory.setLastMetrics(psfMetrics);
    }

    public void listHistory(
571
    		boolean useLMA,
Andrey Filippov's avatar
Andrey Filippov committed
572 573 574 575 576 577 578
    		String path,
    		String lensSerial,
    		String comment,
    		boolean showIndividualComponents,
    		boolean showIndividualSamples,
    		double weightRatioRedToGreen,
    		double weightRatioBlueToGreen,
579
			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
580 581 582
			double weightY // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
			){
    	listHistory(
583
    			useLMA,
Andrey Filippov's avatar
Andrey Filippov committed
584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599
    			path,
    			"",
    			lensSerial,
    			comment,
    			showIndividualComponents,
    			showIndividualSamples,
    			weightRatioRedToGreen,
    			weightRatioBlueToGreen,
    			weightK,
    			weightY,
    			false,
    			Double.NaN,
    			Double.NaN,
    			Double.NaN,
    			Double.NaN);
    }
600

Andrey Filippov's avatar
Andrey Filippov committed
601
    public void listHistory(
602
    		boolean useLMA,
Andrey Filippov's avatar
Andrey Filippov committed
603 604 605 606 607 608 609 610
    		String path,
    		String serialNumber,
    		String lensSerial,
    		String comment,
    		boolean showIndividualComponents,
    		boolean showIndividualSamples,
    		double weightRatioRedToGreen,
    		double weightRatioBlueToGreen,
611
			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
612 613 614 615 616 617 618 619 620
			double weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
			boolean justSummary,
			double pX0,
			double pY0,
			double result_lastKT,
			double result_allHistoryKT
			){

    	this.focusingHistory.list(
621
        		useLMA,
Andrey Filippov's avatar
Andrey Filippov committed
622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637
    			path,
        		serialNumber,
    			lensSerial,
    			comment,
    			showIndividualComponents,
    			showIndividualSamples,
    			weightRatioRedToGreen,
    			weightRatioBlueToGreen,
    			weightK,
    			weightY,
    			justSummary,
    			pX0,
    			pY0,
    			result_lastKT,
    			result_allHistoryKT);
    }
638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654
	public void saveHistoryXML(
			String path,
    		String serialNumber,
			String lensSerial, // if null - do not add average
    		String comment,
			double pX0,
			double pY0,
			double [][][] sampleCoord){ // x,y,r
		this.focusingHistory.saveXML(
    			path,
        		serialNumber,
    			lensSerial, // if null - do not add average
        		comment,
    			pX0,
    			pY0,
    			sampleCoord);
	}
655

656 657
	public FocusingField.FocusingFieldMeasurement getThisFFMeasurement(FocusingField focusingField){
		return getThisFFMeasurement(focusingField, -1);
658
	}
659 660 661 662 663 664 665 666 667 668 669 670
	public FocusingField.FocusingFieldMeasurement getThisFFMeasurement(FocusingField focusingField, int index){
		return focusingField. getFocusingFieldMeasurement(
				historyGetTimestamp(index),   //focusingState.getTimestamp(),
				historyGetTemperature(index), //focusingState.getTemperature(),
				historyGetMotors(index),      //focusingState.motorsPos,
				historyGetSamples(index)); //focusingState.getSamples());
	}
	private FocusingHistory.FocusingState getFocusingState(int index){
		if (index<0) index+=historySize();
		if ((index>=0) && (index<historySize())) return this.focusingHistory.history.get(index);
		return null;
	}
671

672 673 674 675 676
	public void addCurrentHistoryToFocusingField(
			FocusingField focusingField,
			int start,
			int end){
		for (int i=start;i<end;i++){
677
			addCurrentHistoryToFocusingField(focusingField,i);
678 679
		}
	}
680 681 682 683 684 685 686

	public void addCurrentHistoryToFocusingField(FocusingField focusingField){
		addCurrentHistoryToFocusingField(
				focusingField,
				0,
				historySize()
				);
687
	}
688 689 690 691 692 693 694 695 696 697

	public void addCurrentHistoryToFocusingField(
			FocusingField focusingField,
			int index){ // -1 - last (negative - from length)
		focusingField.addSample(
				historyGetTimestamp(index),   //focusingState.getTimestamp(),
				historyGetTemperature(index), //focusingState.getTemperature(),
				historyGetMotors(index),      //focusingState.motorsPos,
				historyGetSamples(index)); //focusingState.getSamples());
	}
698 699


700 701 702 703 704 705 706 707 708
	public String historyGetTimestamp(int index){
		return getFocusingState(index).getTimestamp();
	}

	public double historyGetTemperature(int index){
		return getFocusingState(index).getTemperature();
	}
	public int [] historyGetMotors(int index){
		return getFocusingState(index).getMotors();
709
	}
710

711 712 713 714 715
	public double [][][][] historyGetSamples(int index){
		return getFocusingState(index).getSamples();
	}

	public int historySize(){
Andrey Filippov's avatar
Andrey Filippov committed
716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783
    	return this.focusingHistory.history.size();
    }
    public void setLastProbed(){
    	this.focusingHistory.setLastProbed();
    }
    public double distFromProbed(){ // may return NaN
    	return this.focusingHistory.distFromProbed();
    }
    public double distFromProbed(int [] position){ // may return NaN
    	return this.focusingHistory.distFromProbed(position);
    }
    /**
     * Asks for the motor positions, moves motors, waits for the movement complete
     * @param newPosition array of 3 integers or null (will prompt current position
     * @param useHisteresis 1-element array (will be modified) - apply or not anti-hysteresis moves
     * @param showAuto if true, will show button "Auto", if it will be pressed - no movement
     * @return 0 - cancel, 1 - OK (moved), 2 - pressed "Auto"
     */
    public int dialogMove(
    		int [] newPosition, // null OK
    		boolean [] useHysteresis,
    		boolean showAuto
    ){
    	readElphel10364Motors();
    	GenericDialog gd = new GenericDialog("Moving Focus Motors");
    	int [] promptPosition=(newPosition==null)?this.curpos.clone():newPosition.clone();
    	for (int i=0;i<promptPosition.length;i++){
    		gd.addNumericField("Motor "+(i+1)+" "+(this.lastDirs[i]?"(+)":"(-)")+
    				((newPosition==null)?"":(" ("+this.curpos[i]+")"))+":", promptPosition[i], 0,5,"steps");
    	}
    	if (useHysteresis!=null) gd.addCheckbox("Use anti-hysteresis moves", useHysteresis[0]);
    	if (showAuto) gd.enableYesNoCancel("OK", "Auto");
    	gd.showDialog();
    	if (gd.wasCanceled()) return 0;
    	for (int i=0;i<promptPosition.length;i++){
    		promptPosition[i]= (int) gd.getNextNumber();
    	}
    	if (useHysteresis!=null)  useHysteresis[0]=gd.getNextBoolean();
    	else {
    		useHysteresis=new boolean[1];
    		useHysteresis[0]=false;
    	}
    	if (!gd.wasOKed()) return 2;
    	moveElphel10364Motors(
    			true, //boolean wait,
    			promptPosition,
    			0.0, //double sleep,
    			true, //boolean showStatus,
    			"",   //String message,
    			useHysteresis[0]); //boolean hysteresis)
    	// Ask something else
    	return 1;
    }

    public int dialogFocusing(
    		boolean isAdjusted,
    		int [] newPosition, // null OK
    		int dialogMode,  // 1 - manual, 2 auto first, 3 - auto normal, 4 fine focus (first), 5 fine focus (auto)
    		LensAdjustment.FocusMeasurementParameters focusMeasurementParameters
    ){
    	readElphel10364Motors();
    	GenericDialog gd = new GenericDialog("Moving Focus Motors  - mode="+dialogMode);
    	gd.addMessage("History size - "+historySize());
    	if (historySize()>0){
    		FocusingHistory.FocusingState focusingState=this.focusingHistory.history.get(historySize()-1);
			double fDistance=this.focusingHistory.getLensDistance(
					focusingState.getCenterResolutions(), // {R-sharpness,G-sharpness,B-sharpness}
	    			true, // boolean absolute, // return absolutely calibrated data
784
	    			focusMeasurementParameters.lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811
	    			focusMeasurementParameters.lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
	    			this.debugLevel); // motors debug level
        	if (!Double.isNaN(fDistance)){
        		gd.addMessage("Relative focal distance - "+IJ.d2s(fDistance,3)+" um");
        	}
    		if (historySize()>1){
    			FocusingHistory.FocusingState prevFocusingState=this.focusingHistory.history.get(historySize()-2);
    			double [] xyz={
    					focusingState.motorsPos[0]-prevFocusingState.motorsPos[0],
    					focusingState.motorsPos[1]-prevFocusingState.motorsPos[1],
    					focusingState.motorsPos[2]-prevFocusingState.motorsPos[2]};
    			gd.addMessage("Last travel "+IJ.d2s(Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]),0)+" steps");
    		}
    		double [][] metrics=focusingState.getMetrics(focusMeasurementParameters.weightRatioRedToGreen,focusMeasurementParameters.weightRatioBlueToGreen);
    		for (int i=0;i<metrics.length;i++) {
    			gd.addMessage(focusingState.getName(i)+
    					" Far/Near="+IJ.d2s(metrics[i][0],3)+
    					" TiltX="+IJ.d2s(metrics[i][1],3)+
    					" TiltY="+IJ.d2s(metrics[i][2],3)+
    					" R50%="+IJ.d2s(metrics[i][3],3)+
    					" A50%="+IJ.d2s(metrics[i][4],3)+
    					" B50%="+IJ.d2s(metrics[i][5],3)+
    					" R50%Center="+IJ.d2s(metrics[i][6],3)+
    					((i==3)?(" Weight="+IJ.d2s(100*metrics[i][7],3)+"%"):""));
    		}
    	}
//    	int [] promptPosition=(newPosition==null)?this.curpos.clone():newPosition.clone();
812

Andrey Filippov's avatar
Andrey Filippov committed
813 814 815 816 817 818 819 820 821 822 823 824 825
    	for (int i=0;i<newPosition.length;i++){
    		gd.addNumericField("Motor "+(i+1)+" "+(this.lastDirs[i]?"(+)":"(-)")+
    				((newPosition==null)?"":(" ("+this.curpos[i]+")"))+":", newPosition[i], 0,5,"steps");
    	}
    	gd.addNumericField("Weight ratio  red/green (use when combining defocusing data)",                    focusMeasurementParameters.weightRatioRedToGreen,    2,6,"");
    	gd.addNumericField("Weight ratio blue/green (use when combining defocusing data)",                    focusMeasurementParameters.weightRatioBlueToGreen,    2,6,"");
    	gd.addNumericField("Far/Near Focusing target (logariphm of average tangential-to-radial resolution)", focusMeasurementParameters.targetFarNear,    2,6,"");
		gd.addCheckbox    ("Use targetFarNear (radial/tangential resolution)  as a proxy for the distance",   focusMeasurementParameters.useRadialTangential); // true; // ignore lateral chromatic aberration (center OTF to 0,0)
		gd.addNumericField("Needed lens center distance (away from the \"best focus\"",                       focusMeasurementParameters.targetMicrons,    2,6,"um");
		gd.addNumericField("Tolerance of the center focal distance",                                          focusMeasurementParameters.toleranceMicrons,    2,6,"um");
		gd.addNumericField("Tolerance of the tilt ",                                                          focusMeasurementParameters.toleranceTilt,    2,6,"");
		gd.addNumericField("When under this scaled tolerance, correction step is scaled down twice",          focusMeasurementParameters.toleranceThreshold,    2,6,"");
		gd.addNumericField("Adjust 3 motors by parallel movement if focal distance error in the center exceeds this value", focusMeasurementParameters.parallelAdjustThreshold,    2,6,"um");
826

Andrey Filippov's avatar
Andrey Filippov committed
827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850
    	gd.addNumericField("When fitting far/near and tilts, samples 'decay' with this sigma",                focusMeasurementParameters.motorsSigma,      1,7,"motors steps");
		gd.addNumericField("Same for the 3 motors moving together (center focusing)",                         focusMeasurementParameters.motorsSigma3,      1,7,"motors steps");
    	gd.addNumericField("Low limit for motor sigma when fitting is in the final stage",                    focusMeasurementParameters.motorsMinSigma,      1,7,"motors steps");
    	gd.addNumericField("When walk is getting smaller, sigma be reduced proportionally",                   focusMeasurementParameters.motorsVarSigmaToTravel,      1,7,"motors steps");
    	gd.addNumericField("After each step new sigma will have this part of the calculated from the travel", focusMeasurementParameters.motorsFadeSigma,      1,7,"motors steps");
    	gd.addNumericField("For quadratic maximum the correction will be increased by (1+motorsOverShootToBalance) if there are less samples on the other side", focusMeasurementParameters.motorsOverShootToBalance,      1,7,"motors steps");
    	gd.addNumericField("Maximal allowed single-step focusing adjustment",                                   focusMeasurementParameters.maxStep,      1,7,"motors steps");
    	gd.addNumericField("How far to go to probe around the current point to measure derivatives",            focusMeasurementParameters.probeStep,      1,7,"motors steps");
    	gd.addCheckbox    ("Probe 6 measurements, if unchecked - only 4 (tetrahedron)",                      focusMeasurementParameters.probeSymmetrical);
    	gd.addNumericField("Re-run probing around in orthogonal directions, if the current position moved farther from the last probing one", focusMeasurementParameters.reProbeDistance, 1,7,"motors steps");
    	gd.addNumericField("Bias towards the last measurement: 0% use best fit for all, 100% make planes throgh the last point",100*focusMeasurementParameters.believeLast,      1,5,"%");
    	gd.addCheckbox    ("Move motors in the same direction to compensate for the hysteresis",             focusMeasurementParameters.compensateHysteresis);
		gd.addNumericField("Motor anti-hysteresis travel",                                                   focusMeasurementParameters.motorHysteresis, 0,7,"motors steps");

		gd.addNumericField("Sensor travel to motors travel (all 3 together), by design it is 4/38~=0.105",   focusMeasurementParameters.linearReductionRatio, 5,7,"");

		gd.addNumericField("Motors debug (1 - show moves, 2 show moves+hysteresis mioves)",                  focusMeasurementParameters.motorDebug,        0);

		gd.addMessage("Parameters to calculate lens distance measurements from 3-color PSF measurements");
		gd.addNumericField("Number of points to tabulate center focus parameters vs. focal distance",        focusMeasurementParameters.lensDistanceNumPoints,        0);
		gd.addNumericField("polynomial degree to approximate center focus parameters vs. focal distance ",   focusMeasurementParameters.lensDistancePolynomialDegree,        0);
		gd.addNumericField("Normalize overall sharpness (that depends on the lens quality and/or PSF parameters) to differential ones",            focusMeasurementParameters.lensDistanceWeightY,        3,5,"");
		gd.addNumericField("Use derivartives-dependent weights for components (1.0), equal weight - 0.0",    focusMeasurementParameters.lensDistanceWeightK,        3,5,"");
		gd.addMessage("");
851

Andrey Filippov's avatar
Andrey Filippov committed
852 853 854 855 856 857 858
    	gd.addNumericField("Minimal correction movement to initiate final series of corrections",            focusMeasurementParameters.minCorr,      1,5,"motors steps");
    	gd.addNumericField("Finish if this number of last corrections where below  minimum (previous input)",focusMeasurementParameters.numFinalCorr,      0);
    	gd.addNumericField("Unconditionally exit focusing adjustment after these number of itereations",     focusMeasurementParameters.maxAutoIterations,        0);
    	gd.addNumericField("Maximal allowed total motors travel for automatic correction",                   focusMeasurementParameters.maxAutoDistance, 1,7,"motors steps");
    	gd.addNumericField("If there are insufficient measurements to fit parabola - make this step",        focusMeasurementParameters.maxLinearStep, 1,7,"motors steps");
    	gd.addCheckbox    ("Run fine-focus adjustment for the center samples",                               (dialogMode==4)||(dialogMode==4));

859
    	// TODO: move to other dialog
Andrey Filippov's avatar
Andrey Filippov committed
860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890
    	gd.addCheckbox    ("Scan motors for the focus in the center samples",                                false);
    	gd.addNumericField("Motor single movement (all 3 motors) in scan focus mode (signed value)",         focusMeasurementParameters.scanStep, 0,7,"motors steps");
    	gd.addNumericField("Number of scan steps during (center) focus scanning",                            focusMeasurementParameters.scanNumber,        0);
		gd.addCheckbox    ("Scan focus in 2 directions, after the calibration estimate hysteresis (play)",   focusMeasurementParameters.scanHysteresis);
		gd.addNumericField("Number of scan steps during hysteresis (play) measurement",                      focusMeasurementParameters.scanHysteresisNumber, 0);

    	if      (dialogMode==1) gd.enableYesNoCancel("OK", "Auto Suggest");
    	else if (dialogMode==2) gd.enableYesNoCancel("Apply and Continue", "Retry");
    	else if (dialogMode==3) {
    		if (isAdjusted) gd.enableYesNoCancel("OK/Done", "Single Step"); // here OK - same as cancel?
    		else            gd.enableYesNoCancel("Continue", "Single Step");
    	}
    	else if (dialogMode==4) gd.enableYesNoCancel("Apply and Continue", "Retry");
    	else if (dialogMode==5) gd.enableYesNoCancel("Continue", "Single Step");
    	else {
    		String msg="Mode should be 1,2,3,4 or 5. Got "+dialogMode;
    		IJ.showMessage("Error",msg);
    		throw new IllegalArgumentException (msg);
    	}
		WindowTools.addScrollBars(gd);
    	gd.showDialog();
    	if (gd.wasCanceled()) return 0;
    	int [] wasPosition=newPosition.clone();
    	boolean wasModified=false;
    	for (int i=0;i<newPosition.length;i++){
    		newPosition[i]= (int) gd.getNextNumber();
    		wasModified|=(newPosition[i]!=wasPosition[i]);
    	}
    	focusMeasurementParameters.weightRatioRedToGreen=      gd.getNextNumber();
    	focusMeasurementParameters.weightRatioBlueToGreen=     gd.getNextNumber();
    	focusMeasurementParameters.targetFarNear=              gd.getNextNumber();
891
    	focusMeasurementParameters.useRadialTangential=        gd.getNextBoolean();
Andrey Filippov's avatar
Andrey Filippov committed
892 893 894 895 896 897 898 899
    	focusMeasurementParameters.targetMicrons=              gd.getNextNumber();
    	focusMeasurementParameters.toleranceMicrons=           gd.getNextNumber();
    	focusMeasurementParameters.toleranceTilt=              gd.getNextNumber();
    	focusMeasurementParameters.toleranceThreshold=         gd.getNextNumber();
    	focusMeasurementParameters.parallelAdjustThreshold=    gd.getNextNumber();

    	focusMeasurementParameters.motorsSigma=                gd.getNextNumber();
    	focusMeasurementParameters.motorsSigma3=               gd.getNextNumber();
900
    	focusMeasurementParameters.motorsMinSigma=             gd.getNextNumber();
Andrey Filippov's avatar
Andrey Filippov committed
901
    	focusMeasurementParameters.motorsVarSigmaToTravel=     gd.getNextNumber();
902
    	focusMeasurementParameters.motorsFadeSigma=            gd.getNextNumber();
Andrey Filippov's avatar
Andrey Filippov committed
903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937
    	focusMeasurementParameters.motorsOverShootToBalance=   gd.getNextNumber();

    	focusMeasurementParameters.maxStep=                    gd.getNextNumber();
    	focusMeasurementParameters.probeStep=                  gd.getNextNumber();
    	focusMeasurementParameters.probeSymmetrical=           gd.getNextBoolean();
    	focusMeasurementParameters.reProbeDistance=            gd.getNextNumber();
    	focusMeasurementParameters.believeLast=           0.01*gd.getNextNumber(); if (focusMeasurementParameters.believeLast>1.0) focusMeasurementParameters.believeLast=1.0; else if (focusMeasurementParameters.believeLast<0) focusMeasurementParameters.believeLast=0.0;
    	focusMeasurementParameters.compensateHysteresis=       gd.getNextBoolean();
    	focusMeasurementParameters.motorHysteresis=      (int) gd.getNextNumber();
    	focusMeasurementParameters.linearReductionRatio=       gd.getNextNumber();

    	focusMeasurementParameters.motorDebug=           (int) gd.getNextNumber();
    	focusMeasurementParameters.lensDistanceNumPoints=(int) gd.getNextNumber();
    	focusMeasurementParameters.lensDistancePolynomialDegree=(int) gd.getNextNumber();
    	focusMeasurementParameters.lensDistanceWeightY=        gd.getNextNumber();
    	focusMeasurementParameters.lensDistanceWeightK=        gd.getNextNumber();

    	focusMeasurementParameters.minCorr=                    gd.getNextNumber();
    	focusMeasurementParameters.numFinalCorr=         (int) gd.getNextNumber();
    	focusMeasurementParameters.maxAutoIterations=    (int) gd.getNextNumber();
    	focusMeasurementParameters.maxAutoDistance=            gd.getNextNumber();
    	focusMeasurementParameters.maxLinearStep=              gd.getNextNumber();
    	boolean fineFocus=                                     gd.getNextBoolean();

    	boolean scanMode=                                      gd.getNextBoolean();
    	focusMeasurementParameters.scanStep=             (int) gd.getNextNumber();
    	focusMeasurementParameters.scanNumber=           (int) gd.getNextNumber();
    	focusMeasurementParameters.scanHysteresis=            gd.getNextBoolean();
    	focusMeasurementParameters.scanHysteresisNumber= (int) gd.getNextNumber();

    	setHysteresis(focusMeasurementParameters.motorHysteresis);
    	setDebug(focusMeasurementParameters.motorDebug);

    	if (fineFocus){
    		if (gd.wasOKed()) return 7; // run fine focus
938
    		else return 6;              // run fine focus once, suggest result
Andrey Filippov's avatar
Andrey Filippov committed
939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955
    	}
    	if (scanMode) return 5; // regardless of the non-cancel buttons
    	if (gd.wasOKed()) {
    		if      (dialogMode==1) return 1;
    		else if (dialogMode==2) return 3;
    		else if (dialogMode==3) {
    			if (!isAdjusted) return 4;
    			else if (wasModified) return 1; // manual mode
    			else return 0; // all done
    		}
    		else if (dialogMode==4) return 4;
    		else if (dialogMode==5) return 4;
    	} else {
    		if      (dialogMode==1) return 2;
    		else if (dialogMode==2) return 2;
    		else if (dialogMode==3) return 2;
    	}
956
    	return 0; // should never get here
Andrey Filippov's avatar
Andrey Filippov committed
957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976
    }
    public int dialogFocusingSharpness
    (
    		boolean isAdjusted,
    		int [] newPosition, // null is not OK (may be modified)
    		int dialogMode, // 1 - manual, 2 auto first, 3 - auto normal
    		LensAdjustment.FocusMeasurementParameters focusMeasurementParameters
    ){
    	readElphel10364Motors();
    	GenericDialog gd = new GenericDialog("Moving Focus Motors - mode="+dialogMode);
    	if (isAdjusted ){ // should be true: (focusingSharpness.size()>0)
    		int [] currentPos=focusingSharpness.getPosition(); // last
    		double offCenter= 0.001*getMicronsPerStep()*(currentPos[0]+currentPos[1]+currentPos[2])/3.0;
    		String msg1="Lens focus is approximately "+IJ.d2s(Math.abs(offCenter),3)+" mm "+
			((offCenter>0)?"closer to":"farther from")+" the lens than the motor zero point.";
    		String msg2="You may screw the lens "+((offCenter>0)?"OUT":"IN")+ " by "+IJ.d2s(Math.abs(2.0*offCenter),3)+" turns ("+
			IJ.d2s(Math.abs(360*2.0*offCenter),0)+" degrees)";
    		if (this.debugLevel>0)System.out.println(msg1+"\n"+msg2);
    		gd.addMessage(msg1);
    		gd.addMessage(msg2);
977

Andrey Filippov's avatar
Andrey Filippov committed
978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026
    	}
    	gd.addMessage("History size - "+focusingSharpness.size());
    	if (focusingSharpness.size()>0){
    		double sharpness=focusingSharpness.getSharpness(); // last
    		int [] currentPos=focusingSharpness.getPosition(); // last
    		if (focusingSharpness.size()>1){
    			int [] prevPos=focusingSharpness.getPosition(focusingSharpness.size()-2); // pre-last
    			double [] xyz={
    					currentPos[0]-prevPos[0],
    					currentPos[1]-prevPos[1],
    					currentPos[2]-prevPos[2]};
    			//    				  gd.addMessage("Last travel "+IJ.d2s(Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]),0)+" steps");
    			gd.addMessage("Last center travel "+IJ.d2s((xyz[0]+xyz[1]+xyz[2])/3.0,1)+" steps");
    		}
    		gd.addMessage("Focus sharpness="+sharpness);
    	}
    	for (int i=0;i<newPosition.length;i++){
    		gd.addNumericField("Motor "+(i+1)+" "+(this.lastDirs[i]?"(+)":"(-)")+
    				" ("+this.curpos[i]+")"+":", newPosition[i], 0,5,"steps");
    	}
    	gd.addCheckbox    ("Move motors in the same direction to compensate for the hysteresis",             focusMeasurementParameters.compensateHysteresis);
		gd.addNumericField("Motor anti-hysteresis travel",                                                   focusMeasurementParameters.motorHysteresis, 0,7,"motors steps");
		gd.addNumericField("Motors debug (1 - show moves, 2 show moves+hysteresis mioves)",                  focusMeasurementParameters.motorDebug,        0);
    	gd.addNumericField("Minimal correction movement to initiate final series of corrections",            focusMeasurementParameters.minCorr,      1,5,"motors steps");
    	gd.addNumericField("Finish if this number of last corrections where below  minimum (previous input)",focusMeasurementParameters.numFinalCorr,      0);
    	gd.addNumericField("Unconditionally exit focusing adjustment after these number of itereations",     focusMeasurementParameters.maxAutoIterations,        0);
    	gd.addNumericField("Maximal allowed total motors travel for automatic correction",                   focusMeasurementParameters.maxAutoDistance, 1,7,"motors steps");
    	gd.addNumericField("when fitting parabola for focusing sharpness in the center, far measurements decay with this sigma", focusMeasurementParameters.motorsPreSigma, 1,7,"motors steps");
    	gd.addNumericField("If there are insufficient measurements to fit parabola - make this step",        focusMeasurementParameters.maxLinearStep, 1,7,"motors steps");

    	if      (dialogMode==1) gd.enableYesNoCancel("OK", "Auto Suggest");
    	else if (dialogMode==2) gd.enableYesNoCancel("Apply and Continue", "Retry");
    	else if (dialogMode==3){
    		if (isAdjusted) gd.enableYesNoCancel("OK/Done", "Single Step"); // here OK - same as cancel?
    		else            gd.enableYesNoCancel("Continue", "Single Step");
    	}
    	else {
    		String msg="Mode shou8ld be 1,2 or 3. Got "+dialogMode;
    		IJ.showMessage("Error",msg);
    		throw new IllegalArgumentException (msg);
    	}
    	gd.showDialog();
    	if (gd.wasCanceled()) return -1;
    	int [] wasPosition=newPosition.clone();
    	boolean wasModified=false;
    	for (int i=0;i<newPosition.length;i++){
    		newPosition[i]= (int) gd.getNextNumber();
    		wasModified|=(newPosition[i]!=wasPosition[i]);
    	}
1027

Andrey Filippov's avatar
Andrey Filippov committed
1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052
    	focusMeasurementParameters.compensateHysteresis=       gd.getNextBoolean();
    	focusMeasurementParameters.motorHysteresis=      (int) gd.getNextNumber();
    	focusMeasurementParameters.motorDebug=           (int) gd.getNextNumber();
    	focusMeasurementParameters.minCorr=                    gd.getNextNumber();
    	focusMeasurementParameters.numFinalCorr=         (int) gd.getNextNumber();
    	focusMeasurementParameters.maxAutoIterations=    (int) gd.getNextNumber();
    	focusMeasurementParameters.maxAutoDistance=            gd.getNextNumber();
    	focusMeasurementParameters.motorsPreSigma=             gd.getNextNumber();
    	focusMeasurementParameters.maxLinearStep=              gd.getNextNumber();
    	setHysteresis(focusMeasurementParameters.motorHysteresis);
    	setDebug(focusMeasurementParameters.motorDebug);

    	if (gd.wasOKed()) {
    		if      (dialogMode==1) return 1;
    		else if (dialogMode==2) return 3;
    		else if (dialogMode==3) {
    			if (!isAdjusted) return 4;
    			else if (wasModified) return 1; // manual mode
    			else return 0; // all done
    		}
    	} else {
    		if      (dialogMode==1) return 2;
    		else if (dialogMode==2) return 2;
    		else if (dialogMode==3) return 2;
    	}
1053
    	return 0; // should never get here
Andrey Filippov's avatar
Andrey Filippov committed
1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108
    }

    public void clearPreFocus(){
    	this.focusingSharpness.clear();
    }
    public void listPreFocus(String path, String comment){ // null - on screen
    	this.focusingSharpness.list(path,comment);
    }
    //   				  ,

    public void addPreFocus(
    		double sharpness){
    	addPreFocus(this.curpos, sharpness);
    }

    public void addPreFocus(
    		int [] motorsPos,
    		double sharpness){
    	this.focusingSharpness.add(motorsPos,sharpness);
    }
    public int sizePreFocus(){
    	return this.focusingSharpness.size();
    }
    public int sizeFocus(){
    	return this.focusingHistory.size();
    }
    public int [] solvePreFocus(
    		double motorsSigma,
    		double maxStep,        // Maximal motor move for each step
    		double linearStep,     // when no quadratic maximum
			double motorsOverShootToBalance,
    		int debugLevel
    ){
    	return focusingSharpness.solve(
    			getRange(),
    			motorsSigma,
    			maxStep,        // Maximal motor move for each step
    			linearStep,     // when no quadratic maximum
    			motorsOverShootToBalance,
    			debugLevel
    	);

    }

    public int [] solveSinglePoly(
    		boolean first,
    		double wR,  // weight of red in focus sharpness
    		double wG,
    		double wB,
    		boolean sameTiltOnly, // us only history where the difference between motors was the same as in the last sample
    		boolean centerOnly, // only use center samples
    		int polyDegree,      // polynomial order (just 2 always)
    		double motorsSigma,
    		double maxStep,        // Maximal motor move for each step
    		double linearStep,     // when no quadratic maximum
1109
    		double motorsMinSigma,
Andrey Filippov's avatar
Andrey Filippov committed
1110
    		double motorsVarSigmaToTravel,
1111
    		double motorsFadeSigma,
Andrey Filippov's avatar
Andrey Filippov committed
1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126
    		double motorsOverShootToBalance,
    		int debugLevel
    ){
    	return focusingHistory.solveSinglePoly(
    			first,
    			wR,  // weight of red in focus sharpness
    			wG,
    			wB,
    			sameTiltOnly, // us only history where the difference between motors was the same as in the last sample
    			centerOnly, // only use center samples
    			polyDegree,      // polynomial order (just 2 always)
    			getRange(),
    			motorsSigma,
    			maxStep,        // Maximal motor move for each step
    			linearStep,     // when no quadratic maximum
1127
    			motorsMinSigma,
Andrey Filippov's avatar
Andrey Filippov committed
1128
    			motorsVarSigmaToTravel,
1129
    			motorsFadeSigma,
Andrey Filippov's avatar
Andrey Filippov committed
1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172
    			motorsOverShootToBalance,
    			debugLevel
    	);

    }


    public class FocusingSharpness{
    	public ArrayList<FocusingSharpnessItem> history=null;
    	public FocusingSharpness(int size){
    		history=new ArrayList <FocusingSharpnessItem>(size);
    	}
    	public void clear(){
    		history.clear();
    	}
    	public int [] getPosition(){
    		if (this.history.size()<1) return null;
    		return this.history.get(this.history.size()-1).motorsPos;
    	}
    	public int [] getPosition(int index){
    		if ((index<0) || (index>=this.history.size())) return null;
    		return this.history.get(index).motorsPos;
    	}
    	public int size(){
    		return this.history.size();
    	}
    	public double getSharpness(){
    		if (this.history.size()<1) return Double.NaN;
    		return this.history.get(this.history.size()-1).sharpness;
    	}
    	public double getSharpness(int index){
    		if ((index<0) || (index>=this.history.size())) return Double.NaN;
    		return this.history.get(index).sharpness;
    	}
    	public void add(
    			int [] motorsPos,
    			double sharpness
    	){
    		this.history.add(new FocusingSharpnessItem(motorsPos,sharpness));
    	}

    	// moves motors all together, using average as a single parameter.
    	// tries a second-degree polynomial, if not enough data - linear (maxSteop)
1173
    	// TODO: modify to use PolunomialApproximation class
Andrey Filippov's avatar
Andrey Filippov committed
1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297
    	public int [] solve(
    			double [][] motorRange,
    			double motorsSigma,
    			double maxStep,        // Maximal motor move for each step
    			double linearStep,     // when no quadratic maximum
    			double motorsOverShootToBalance,
    			int debugLevel
    	){
    		double S0=0.0,SX=0.0,SX2=0.0,SX3=0.0,SX4=0.0,SF=0.0,SFX=0.0,SFX2=0.0;
    		double rangeCenter=(
    				motorRange[0][0]+motorRange[0][1]+
    				motorRange[1][0]+motorRange[1][1]+
    				motorRange[2][0]+motorRange[2][1])/6.0;
    		int [] origin=getPosition();
    		if (
    				(origin[0]<=motorRange[0][0]) ||
    				(origin[0]>=motorRange[0][1]) ||
    				(origin[1]<=motorRange[1][0]) ||
    				(origin[1]>=motorRange[1][1]) ||
    				(origin[2]<=motorRange[2][0]) ||
    				(origin[2]>=motorRange[2][1])) {
    			String message="Motors are off range: ["+origin[0]+", "+origin[1]+", "+origin[2]+"], range is ["+
    			motorRange[0][0]+".."+motorRange[0][1]+", "+
    			motorRange[1][0]+".."+motorRange[1][1]+", "+
    			motorRange[2][0]+".."+motorRange[2][1]+"]";
    			System.out.println ("*** ERROR *** "+message);
    			IJ.showMessage("ERROR",message);
    			return null;
    		}
    		double x0=(origin[0]+origin[1]+origin[2])/3.0;
    		double move=0;
    		if (size()==1){
    			move=(x0>rangeCenter)?-linearStep:linearStep;
    		} else {

    			double kExp=0.5/motorsSigma/motorsSigma;
    			for (int i=0;i<size();i++){
    				int [] position=getPosition(i);
    				double f=getSharpness(i);
    				double x=((position[0]+position[1]+position[2])/3.0-x0);
    				double w=Math.exp(-kExp*(x*x));
    				double x2=x*x;
    				double x3=x2*x;
    				S0+=  w;
    				SX+=  w*x;
    				SX2+= w*x2;
    				SX3+= w*x3;
    				SX4+= w*x3*x;
    				SF+=  w*f;
    				SFX+= w*f*x;
    				SFX2+=w*f*x2;
    			}
    			Matrix M,B;
    			if (size()<3){
    				double [][] aM={
    						{SX2,SX },
    						{SX, S0 }};
    				double [][] aB={
    						{SFX},
    						{SF}};
    				M=new Matrix(aM);
    				B=new Matrix(aB);
    			} else {
    				double [][] aM={
    						{SX4,SX3,SX2},
    						{SX3,SX2,SX },
    						{SX2,SX, S0 }};
    				double [][] aB={
    						{SFX2},
    						{SFX},
    						{SF}};
    				M=new Matrix(aM);
    				B=new Matrix(aB);

    				if (!(new LUDecomposition(M)).isNonsingular()){
    					if (debugLevel>0){
    						System.out.println("Singular matrix - using linear approximation");
    					}
    					double [][] aM2={
    							{SX2,SX },
    							{SX, S0 }};
    					double [][] aB2={
    							{SFX},
    							{SF}};
    					M=new Matrix(aM2);
    					B=new Matrix(aB2);
    				}

    			}
    			if (debugLevel>1){
    				System.out.println("M:");
    				M.print(10, 5);
    				System.out.println("B:");
    				B.print(10, 5);
    			}
    			double [] abc= M.solve(B).getColumnPackedCopy();
    			if (abc.length<3){
    				double [] tmp={0.0,abc[0],abc[1]};
    				abc=tmp;
    			}
    			if (debugLevel>1){
    				System.out.println("A:");
    				M.solve(B).print(10, 5);
    			}
    			move=0.0;
    			if (abc[0]>=0.0){ // no maximum at all
    				double pos=abc[0]*linearStep*linearStep+abc[1]*linearStep+abc[2];
    				double neg=abc[0]*linearStep*linearStep-abc[1]*linearStep+abc[2];
    				move=((pos>neg)?1:-1)*linearStep;
    			} else {
    				move=-abc[1]/(2*abc[0]);
//    				System.out.println("move="+move+" SX="+SX);
					if (move*SX<0){ // only balance in one direction and if center is "behind"
						double balancedMove=move-(SX-move)*S0;
						if (debugLevel>0) {
							System.out.println("Balancing (first estimation): move="+move+" balancedMove="+balancedMove);
						}
						if ((balancedMove/move)>motorsOverShootToBalance) balancedMove=move*(1.0+motorsOverShootToBalance);
						move= (int) Math.round(balancedMove);
						if (debugLevel>0) {
							System.out.println("Balancing limited move="+move);
						}
					}
    				if (Math.abs(move)>maxStep) move*=maxStep/Math.abs(move);
1298

Andrey Filippov's avatar
Andrey Filippov committed
1299 1300
    			}
    		}
1301
    		// limit to the range
Andrey Filippov's avatar
Andrey Filippov committed
1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320
    		for (int m=0; m<motorRange.length;m++) {
    			if (move>0){
    				if ((origin[m]+move)>motorRange[m][1]) move=motorRange[m][1]-origin[m];
    			} else {
    				if ((origin[m]+move)<motorRange[m][0]) move=motorRange[m][0]-origin[m];
    			}
    		}
    		int [] result = {
    				origin[0]+(int) move,
    				origin[1]+(int) move,
    				origin[2]+(int) move};
    		// see if it is an old result - then double the step?
    		boolean newPoint=true;
    		for (int i=0;i<size();i++){
    			int [] position=getPosition(i);
    			if ((position[0]==result[0])&&(position[1]==result[1])&&(position[2]==result[2])){
    				newPoint=false;
    				break;
    			}
1321
    		}
Andrey Filippov's avatar
Andrey Filippov committed
1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390
    		if (!newPoint) {
    			if (debugLevel>0) {
    				System.out.println("Got to the same point: ["+result[0]+", "+result[1]+", "+result[2]+"] again, doubling the step");
    			}
    			for (int m=0; m<motorRange.length;m++) {
    				if (move>0){
    					if ((result[m]+move)>motorRange[m][1]) move=motorRange[m][1]-result[m];
    				} else {
    					if ((result[m]+move)<motorRange[m][0]) move=motorRange[m][0]-result[m];
    				}
    			}
    			for (int m=0; m<motorRange.length;m++) result[m]+=move;
    		}
    		return result;
    	}


    	public void list(String path, String comment){
    		String header="#\tcenter travel\tCenter\tMotor1\tMotor2\tMotor3\tSharpness";
    		StringBuffer sb = new StringBuffer();
    		int [] prevPos=null;
    		for (int i=0;i<size();i++){
    			int [] motors=getPosition(i);
    			double sharpness=getSharpness(i);
    			sb.append((i+1)+"\t");
    			if (prevPos!=null) {
    				double [] xyz={motors[0]-prevPos[0],motors[1]-prevPos[1],motors[2]-prevPos[2]};
    				//	    	    		sb.append(""+IJ.d2s(Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]),0));
    				sb.append(""+IJ.d2s((xyz[0]+xyz[1]+xyz[2])/3.0,0)); // signed
    			}
    			sb.append("\t"+IJ.d2s((motors[0]+motors[1]+motors[2])/3.0,1)+"\t"+motors[0]+"\t"+motors[1]+"\t"+motors[2]+"\t"+IJ.d2s(sharpness,3)+"\n");
    			prevPos=motors;
    		}
    		if (path!=null) {
    			String footer="Comment: "+comment;
    			CalibrationFileManagement.saveStringToFile (
						path,
						header+"\n"+sb.toString()+"\n"+footer);
    		}else{
    			new TextWindow("Sharpness History", header, sb.toString(), 500,900);
    		}
    	}



    	public class FocusingSharpnessItem{
    		public int [] motorsPos=new int[3];
    		public double sharpness;
    		FocusingSharpnessItem(
    				int [] motorsPos,
    				double sharpness
    		){
    			this.motorsPos=motorsPos.clone();
    			this.sharpness=sharpness;
    		}
    	}
    }

    public class FocusingHistory{
    	private ArrayList<FocusingState> history=new ArrayList <FocusingState>(100);
    	private FocusingState focusingState = new FocusingState(); // just to access constant arrays
    	private double currentSigma=Double.NaN; // saves floating motor sigma between calls
    	private double [][]lensDistanceCalibration=null; // array to convert 3 color sharpnesses in the lens center to the distance
    	private double [][]polyCoeff=null; // polynomial approximations of the center color fractions and intensity
    	private double corrXMin;  // range of average motor coordinates to consider (should have the same tilt as last in history)
    	private double corrXMax;
    	private double micronsPerStep=Double.NaN;
    	private int    calibrationNumPoints; // number of calibration points in the table approximation;

1391
    	private int [] calibrationHistoryFromTo=null; // history indices, to-from (inclusive) used to calculate calibration
Andrey Filippov's avatar
Andrey Filippov committed
1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411
    	private int [] hysteresisHistoryFromTo=null; // history indices, to-from (inclusive) used to calculate hysteresis
    	private double measuredHysteresis=Double.NaN;
	    public int indexR=0; // index for Y data
	    public int indexB=1;
	    public int indexY=2;
	    public int indexRW=3;
	    public int indexBW=4;
	    public int indexYW=5;
	    public int indexSteps=6;
	    public int indexMicrons=7;

    	public FocusingHistory(int len){
    		this.history=new ArrayList <FocusingState>(len);
    	}
    	public double getMeasuredHysteresis(){ return this.measuredHysteresis;}
    	public void setCalibrationHistoryFromTo(int from, int to){
    		this.calibrationHistoryFromTo=new int [2];
    		this.calibrationHistoryFromTo[0]=from;
    		this.calibrationHistoryFromTo[1]=to;
    	}
1412

Andrey Filippov's avatar
Andrey Filippov committed
1413 1414 1415 1416 1417
    	public void setHysteresisHistoryFromTo(int from, int to){
    		this.hysteresisHistoryFromTo=new int [2];
    		this.hysteresisHistoryFromTo[0]=from;
    		this.hysteresisHistoryFromTo[1]=to;
    	}
1418

Andrey Filippov's avatar
Andrey Filippov committed
1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434
    	public int [] getPosition(){
    		if (this.history.size()<1) return null;
    		return this.history.get(this.history.size()-1).getMotors();
    	}
    	public int [] getPosition(int index){
    		if ((index<0) || (index>=this.history.size())) return null;
    		return this.history.get(index).getMotors();
    	}
    	public String getTimestamp(){
    		if (this.history.size()<1) return null;
    		return this.history.get(this.history.size()-1).getTimestamp();
    	}
    	public String getTimestamp(int index){
    		if ((index<0) || (index>=this.history.size())) return null;
    		return this.history.get(index).getTimestamp();
    	}
1435 1436


Andrey Filippov's avatar
Andrey Filippov committed
1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466
    	public double [][]  getMetrics(double wR, double wG, double wB){
    		if (this.history.size()<1) return null;
    		return getMetrics(this.history.size()-1,wR,wG,wB);
    	}
    	public double [][]  getMetrics(double wRtoG, double wBtoG){
    		if (this.history.size()<1) return null;
    		return getMetrics(this.history.size()-1,wRtoG,wBtoG);
    	}

    	public double [][]  getMetrics(int index, double wR, double wG, double wB){
    		if ((index<0) || (index>=this.history.size())) return null;
    		return this.history.get(index).getMetrics(wR,wG,wB);
    	}
    	public double [][]  getMetrics(int index, double wRtoG, double wBtoG){
    		if ((index<0) || (index>=this.history.size())) return null;
    		return this.history.get(index).getMetrics(wRtoG,wBtoG);
    	}

    	public double getPositionAverage(){
    		if (this.history.size()<1) return Double.NaN;
    		return this.history.get(this.history.size()-1).getMotorsAverage();
    	}
    	public double  getPositionAverage(int index){
    		if ((index<0) || (index>=this.history.size())) return Double.NaN;
    		return this.history.get(index).getMotorsAverage();
    	}
    	public double []  getCenterResolutions(){
    		if (this.history.size()<1) return null;
    		return getCenterResolutions(this.history.size()-1);
    	}
1467
    	public double [] getCenterResolutions(int index){
Andrey Filippov's avatar
Andrey Filippov committed
1468 1469
    		if ((index<0) || (index>=this.history.size())) return null;
    		return this.history.get(index).getCenterResolutions();
1470 1471 1472 1473 1474 1475 1476 1477 1478
    	}
    	public double [] getzTxTy(int index){
    		if ((index<0) || (index>=this.history.size())) return null;
    		return this.history.get(index).getzTxTy();
    	}
    	public double [] getzTxTy(){
    		if (this.history.size()<1) return null;
    		return getzTxTy(this.history.size()-1);
    	}
Andrey Filippov's avatar
Andrey Filippov committed
1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493
    	public int size(){
    		return this.history.size();
    	}
    	public int [] solveSinglePoly(
    			boolean first,
    			double wR,  // weight of red in focus sharpness
    			double wG,
    			double wB,
    			boolean sameTiltOnly, // us only history where the difference between motors was the same as in the last sample
    			boolean centerOnly, // only use center samples
    			int polyDegree,      // polynomial order (just 2 always)
    			double [][] motorRange,
    			double motorsSigma,
    			double maxStep,        // Maximal motor move for each step
    			double linearStep,     // when no quadratic maximum
1494
    			double motorsMinSigma,
Andrey Filippov's avatar
Andrey Filippov committed
1495
    			double motorsVarSigmaToTravel,
1496
    			double motorsFadeSigma,
Andrey Filippov's avatar
Andrey Filippov committed
1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581
    			double motorsOverShootToBalance,
    			int debugLevel
    	){
    		if (first || Double.isNaN(this.currentSigma)) this.currentSigma=motorsSigma;
    		if (polyDegree>2)polyDegree=2; // for now limit
    		int indexR50=centerOnly?6:3; // 6 is average R50 in the center samples,3 - in all samples
    		double rangeCenter=(
    				motorRange[0][0]+motorRange[0][1]+
    				motorRange[1][0]+motorRange[1][1]+
    				motorRange[2][0]+motorRange[2][1])/6.0;
    		int [] origin=getPosition();
    		if (
    				(origin[0]<=motorRange[0][0]) ||
    				(origin[0]>=motorRange[0][1]) ||
    				(origin[1]<=motorRange[1][0]) ||
    				(origin[1]>=motorRange[1][1]) ||
    				(origin[2]<=motorRange[2][0]) ||
    				(origin[2]>=motorRange[2][1])) {
    			String message="Motors are off range: ["+origin[0]+", "+origin[1]+", "+origin[2]+"], range is ["+
    			motorRange[0][0]+".."+motorRange[0][1]+", "+
    			motorRange[1][0]+".."+motorRange[1][1]+", "+
    			motorRange[2][0]+".."+motorRange[2][1]+"]";
    			System.out.println ("*** ERROR *** "+message);
    			IJ.showMessage("ERROR",message);
    			return null;
    		}
    		boolean []sameTiltmask=new boolean[size()];
    		int sameTiltSize=0;
    		for (int i=0;i<size();i++){
    			if (!sameTiltOnly) {
    				sameTiltmask[i]=true;
    				sameTiltSize++;
    			} else {
    				int [] position=getPosition(i);
    				if (((position[0]-position[1])==(origin[0]-origin[1])) && ((position[2]-position[1])==(origin[2]-origin[1]))){
    					sameTiltmask[i]=true;
    					sameTiltSize++;
    				} else sameTiltmask[i]=false;
    			}
    		}
    		if (debugLevel>1) System.out.println("solveSinglePoly(): sameTiltSize="+sameTiltSize+", size()="+size());
    		double x0=(origin[0]+origin[1]+origin[2])/3.0;
    		int move=0;
    		double w=wR+wG+wB;
    		wR/=w;
    		wG/=w;
    		wB/=w;
    		int iLinearStep=(int) linearStep;
    		if (sameTiltSize==1){
    			move= (x0>rangeCenter)?-iLinearStep:iLinearStep;
    		} else {
    			// create a data set for the polynomial approximation
    			double [][] data=new double [size()][3]; //x,f,w
    			double kExp=0.5/this.currentSigma/this.currentSigma;
    			//sameDiffOnly
    			double SX=0,S0=0;
    			for (int i=0;i<size();i++) if (sameTiltmask[i]){
    				data[i][0]=getPositionAverage(i)- x0;
    				double [][] metrics= this.history.get(i).getMetrics(wR,wG,wB); // here we need to average 1/value, so have to calculate
    				data[i][1]=wR/metrics[0][indexR50]+ wG/metrics[1][indexR50]+wB/metrics[2][indexR50];// "sharpness"
    				data[i][2]=Math.exp(-kExp*(data [i][0]*data [i][0]));
    				SX+=data[i][0]*data[i][2];
    				S0+=data[i][2];
    			} else {
    				data[i][2]=0.0;
    			}
    			if (S0>0) SX/=S0; // center of gravity of the weighted samples used for parabola fitting
    			// the following line did not cause PolynomialApproximation class to recompile after it was modified

    			//        				double [] polyCoeff=(new PolynomialApproximation(debugLevel)).polynomialApproximation1d(data, polyDegree);
    			PolynomialApproximation pa= new PolynomialApproximation(debugLevel);
    			double [] polyCoeff=pa.polynomialApproximation1d(data, polyDegree);

    			move=0;
    			// history can consist of single point repeated several times
    			if ((polyCoeff[2]==0.0) && (polyCoeff[1]==0.0)){
    				move=(x0>rangeCenter)?-iLinearStep:iLinearStep;
    			} else {
    				if (polyCoeff[2]>=0.0){ // no maximum at all

    					double pos=polyCoeff[2]*linearStep*linearStep+polyCoeff[1]*linearStep+polyCoeff[0];
    					double neg=polyCoeff[2]*linearStep*linearStep-polyCoeff[1]*linearStep+polyCoeff[0];
    					move=((pos>neg)?1:-1)*iLinearStep;
    				} else {
    					move=(int) (-polyCoeff[1]/(2*polyCoeff[2]));
1582
    					// modify move to balance around maximum when looking for the quadratic maximum
Andrey Filippov's avatar
Andrey Filippov committed
1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596
    					if (move*SX<0){ // only balance in one direction and if center is "behind"
    						double balancedMove=move-(SX-move)*S0;
    						if (debugLevel>0) {
    							System.out.println("Balancing (first estimation): move="+move+" balancedMove="+balancedMove);
    						}
    						if ((balancedMove/move)>motorsOverShootToBalance) balancedMove=move*(1.0+motorsOverShootToBalance);
    						move= (int) Math.round(balancedMove);
    						if (debugLevel>0) {
    							System.out.println("Balancing limited move="+move);
    						}
    					}
    					if (Math.abs(move)>maxStep) move*=maxStep/Math.abs(move);
    				}
    			}
1597 1598
    		}
    		// limit to the range
Andrey Filippov's avatar
Andrey Filippov committed
1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617
    		for (int m=0; m<motorRange.length;m++) {
    			if (move>0){
    				if ((origin[m]+move)>motorRange[m][1]) move=(int)motorRange[m][1]-origin[m];
    			} else {
    				if ((origin[m]+move)<motorRange[m][0]) move=(int)motorRange[m][0]-origin[m];
    			}
    		}
    		int [] result = {
    				origin[0]+ move,
    				origin[1]+ move,
    				origin[2]+ move};
    		// see if it is an old result - then double the step?
    		boolean newPoint=true;
    		for (int i=0;i<size();i++){
    			int [] position=getPosition(i);
    			if ((position[0]==result[0])&&(position[1]==result[1])&&(position[2]==result[2])){
    				newPoint=false;
    				break;
    			}
1618
    		}
Andrey Filippov's avatar
Andrey Filippov committed
1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677
    		if (!newPoint) {
    			int move0=move;
    			if (debugLevel>0) {
    				System.out.println("Got to the same point: ["+result[0]+", "+result[1]+", "+result[2]+"] again, doubling the step");
    			}
    			for (int m=0; m<motorRange.length;m++) {
    				if (move>0){
    					if ((result[m]+move)>motorRange[m][1]) move=(int) motorRange[m][1]-result[m];
    				} else {
    					if ((result[m]+move)<motorRange[m][0]) move=(int) motorRange[m][0]-result[m];
    				}
    			}
    			for (int m=0; m<motorRange.length;m++) result[m]+=move;
    			move+=move0; // full movement;
    		}
    		// modify sigma
    		this.currentSigma=(1.0-motorsFadeSigma)*this.currentSigma+motorsFadeSigma*(motorsVarSigmaToTravel*Math.abs(move));
    		if      (this.currentSigma>motorsSigma)    this.currentSigma=motorsSigma;
    		else if (this.currentSigma<motorsMinSigma) this.currentSigma=motorsMinSigma;
    		if (debugLevel>0) {
    			System.out.println("Current motor sigma is "+IJ.d2s(this.currentSigma,2));
    		}
    		return result;
    	}

    	public void setProperties(String prefix,Properties properties){
    		if (this.polyCoeff!=null){
    			properties.setProperty(prefix+"corrXMin",this.corrXMin+"");
    			properties.setProperty(prefix+"corrXMax",this.corrXMax+"");
    			properties.setProperty(prefix+"polyNumChannels",(this.polyCoeff.length)+"");
    			properties.setProperty(prefix+"polyDegree",(this.polyCoeff[0].length-1)+"");
    			for (int i=0;i<this.polyCoeff.length;i++) for (int j=0;j< (this.polyCoeff[0].length-1);j++)
        			properties.setProperty(prefix+"polyCoeff_"+i+"_"+j,this.polyCoeff[i][j]+"");
    		}
    	}

		public void getProperties(String prefix,Properties properties){
			if (properties.getProperty(prefix+"corrXMin")!=null)
				this.corrXMin=Double.parseDouble(properties.getProperty(prefix+"corrXMin"));
			if (properties.getProperty(prefix+"corrXMax")!=null)
				this.corrXMax=Double.parseDouble(properties.getProperty(prefix+"corrXMax"));
			if ((properties.getProperty(prefix+"polyNumChannels")!=null) && (properties.getProperty(prefix+"polyDegree")!=null)){
				this.polyCoeff=new double[Integer.parseInt(properties.getProperty(prefix+"polyNumChannels"))]
				                          [Integer.parseInt(properties.getProperty(prefix+"polyDegree"))+1];
    			for (int i=0;i<this.polyCoeff.length;i++) for (int j=0;j< (this.polyCoeff[0].length-1);j++){
    				if (properties.getProperty(prefix+"polyCoeff_"+i+"_"+j)==null){
    					String msg="Missing parameter \"polyCoeff_"+i+"_"+j+" in configuration file";
    					System.out.println("Error: "+msg);
    					IJ.showMessage("Error", msg);
    				} else {
    					this.polyCoeff[i][j]=Double.parseDouble(properties.getProperty(prefix+"polyCoeff_"+i+"_"+j));
    				}
    			}
			}

			// note: needs absolute calibration to convert steps to microns: optimalMotorPosition(
			// LensAdjustment.FocusMeasurementParameters focusMeasurementParameters,
			// double micronsPerStep,
			//  int debugLevel);
1678
	}
Andrey Filippov's avatar
Andrey Filippov committed
1679 1680 1681 1682
		public int [] optimalMotorPosition(
				LensAdjustment.FocusMeasurementParameters focusMeasurementParameters,
				double micronsPerStep,
				int debugLevel
1683
		) {
Andrey Filippov's avatar
Andrey Filippov committed
1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695
			if (this.polyCoeff==null) return null; // not calibrated at all
    		this.calibrationNumPoints=focusMeasurementParameters.lensDistanceNumPoints;
           	makeCalibrationTable(
        			this.polyCoeff,
        			this.calibrationNumPoints,
    				this.corrXMin,  // range of average motor coordinates to consider (should have the same tilt as last in history)
    				this.corrXMax,
        			debugLevel
        			);
			return absoluteCalibrate(
					micronsPerStep,
					focusMeasurementParameters.targetMicrons,
1696
					focusMeasurementParameters.lensDistanceWeightY,
Andrey Filippov's avatar
Andrey Filippov committed
1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715
					focusMeasurementParameters.lensDistanceWeightK,
					debugLevel);
		}

    	public int [] calibrateLensDistance(
    			String path, // or null
        		LensAdjustment.FocusMeasurementParameters focusMeasurementParameters,
    			boolean centerOnly,
    			double micronsPerStep,
    			  int debugLevel
    			){
    		if (this.calibrationHistoryFromTo==null){
    			return null; // not calibrated
    		}
    		if (focusMeasurementParameters.lensDistanceInteractive){
    			GenericDialog gd = new GenericDialog("Lens distance calculation");
    			gd.addNumericField("Number of division points",     focusMeasurementParameters.lensDistanceNumPoints, 0); // 1000
    			gd.addNumericField("Polynomial degree",                focusMeasurementParameters.lensDistancePolynomialDegree, 0); // 8
    			gd.addCheckbox    ("Center samples only",              centerOnly); // true
1716
    			gd.addNumericField("Weight Y",                           focusMeasurementParameters.lensDistanceWeightY, 3,5,""); // 0.5
Andrey Filippov's avatar
Andrey Filippov committed
1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763
    			gd.addNumericField("Weight K",                           focusMeasurementParameters.lensDistanceWeightK, 3,5,""); // 0.5
	    		gd.addCheckbox    ("Show results window from focal distance calibration", focusMeasurementParameters.lensDistanceShowResults);
	    		gd.addNumericField("Needed lens center distance (away from the \"best focus\"", focusMeasurementParameters.targetMicrons,    2,6,"um");

    			gd.addNumericField("Debug level",                       debugLevel,0);
    			gd.showDialog();
    			if (gd.wasCanceled()) return null;
    			focusMeasurementParameters.lensDistanceNumPoints=         (int) gd.getNextNumber();
    			focusMeasurementParameters.lensDistancePolynomialDegree=  (int) gd.getNextNumber();
    			centerOnly=                                                     gd.getNextBoolean();
    			focusMeasurementParameters.lensDistanceWeightY=                 gd.getNextNumber();
    			focusMeasurementParameters.lensDistanceWeightK=                 gd.getNextNumber();
    			focusMeasurementParameters.lensDistanceShowResults=             gd.getNextBoolean();
    			focusMeasurementParameters.targetMicrons=                       gd.getNextNumber();
    			debugLevel=                                               (int) gd.getNextNumber();
    		}
 //   		int sameTiltSize=0;
    		if (this.calibrationHistoryFromTo[0]<0)this.calibrationHistoryFromTo[0]=0;
    		if (this.calibrationHistoryFromTo[1]>=size())this.calibrationHistoryFromTo[1]=size()-1;
    		int [] origin=getPosition(this.calibrationHistoryFromTo[1]);
    		boolean []sameTiltmask=new boolean[size()];
    		for (int i=0;i<sameTiltmask.length;i++) sameTiltmask[i]=false;
    		for (int i=this.calibrationHistoryFromTo[0];i<=this.calibrationHistoryFromTo[1];i++){
    			int [] position=getPosition(i);
    			if (((position[0]-position[1])==(origin[0]-origin[1])) && ((position[2]-position[1])==(origin[2]-origin[1]))){
    				sameTiltmask[i]=true;
 //   				sameTiltSize++;
    			}
    		}
			// need to make sure the limits do not extend beyond actual samples (polynomials may go "crazy" there
    		this.corrXMin=getPositionAverage(this.calibrationHistoryFromTo[0]);
    		this.corrXMax=this.corrXMin;
    		for (int i=this.calibrationHistoryFromTo[0];i<=this.calibrationHistoryFromTo[1];i++) if (sameTiltmask[i]){
				double x=getPositionAverage(i);
				if (x<corrXMin) this.corrXMin=x;
				if (x>corrXMax) this.corrXMax=x;
			}
    		this.corrXMin=Math.floor(this.corrXMin)-1.0;
    		this.corrXMax=Math.ceil(this.corrXMax)+1.0;
    		this.calibrationNumPoints=focusMeasurementParameters.lensDistanceNumPoints;
    		if (debugLevel>0){
    			System.out.println( " Range is "+this.corrXMin+".."+this.corrXMax);
    		}
//    		xMin=corrXMin;
//    		xMax=corrXMax;

    		this.polyCoeff= calibrateFractions(
1764
    				this.calibrationHistoryFromTo[0], // first history data to consider
Andrey Filippov's avatar
Andrey Filippov committed
1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782 1783
        			this.calibrationHistoryFromTo[1],   // last history data to consider
    				focusMeasurementParameters.lensDistancePolynomialDegree,  // polynomial degree to fit
					centerOnly, // only use center samples
					this.corrXMin,  // range of average motor coordinates to consider (should have the same tilt as last in history)
					this.corrXMax,
					debugLevel);
        	makeCalibrationTable(
        			this.polyCoeff,
        			this.calibrationNumPoints,
    				this.corrXMin,  // range of average motor coordinates to consider (should have the same tilt as last in history)
    				this.corrXMax,
        			debugLevel
        			);

    		String header="#\ttimestamp\tcenter\tResol-R\tResol-G\tResol-B\tfrac-R\tfrac-B\tY\tindex\tT-frac-R\tT-frac-B\tT-Y\tT-X";
//    		int indexR50=centerOnly?6:3; // 6 is average R50 in the center samples,3 - in all samples
    		StringBuffer sb = new StringBuffer();
    		double sumErrors=0.0;
    		double numSummedErrors=0;
1784 1785 1786



Andrey Filippov's avatar
Andrey Filippov committed
1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816
			for (int i=this.calibrationHistoryFromTo[0];i<=this.calibrationHistoryFromTo[1];i++)	if (sameTiltmask[i]){
				double x=getPositionAverage(i);
				String timestamp=getTimestamp(i);
				if (timestamp==null) timestamp="n/a";
				sb.append(i+"\t"+timestamp+"\t"+IJ.d2s(x,0));
				double [] sharp=this.history.get(i).getCenterResolutions();

				double l2=sharp[0]*sharp[0] + sharp[1]*sharp[1] + sharp[2]*sharp[2];
				sb.append("\t"+IJ.d2s(sharp[0],3)+"\t"+IJ.d2s(sharp[1],3)+"\t"+IJ.d2s(sharp[2],3)); // r,g,b here
				sb.append("\t"+IJ.d2s(sharp[0]*sharp[0]/l2,3));
				sb.append("\t"+IJ.d2s(sharp[2]*sharp[2]/l2,3));
				sb.append("\t"+IJ.d2s(l2,3));
				int j0=0;
				for (int j=0;j<this.lensDistanceCalibration.length;j++){
					if (this.lensDistanceCalibration[j][6]>x){
						j0=j;
						break;
					}
				}
				 sb.append("\t"+j0);
				if (j0==0) sb.append("\t---\t---\t---\t---\t---\t---\t---");
				else {
					double k= (x-this.lensDistanceCalibration[j0-1][6])/(this.lensDistanceCalibration[j0][6]-this.lensDistanceCalibration[j0-1][6]);
					for (int j=0;j<3;j++){
						double d=k*this.lensDistanceCalibration[j0][j]+(1.0-k)*this.lensDistanceCalibration[j0-1][j];
						sb.append("\t"+IJ.d2s(d,3));
					}
					double d=getLensDistance(
			    			sharp, // {R-sharpness,G-sharpness,B-sharpness}
			    			false, // boolean absolute, // return absolutely calibrated data
1817
			    			focusMeasurementParameters.lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860
			    			focusMeasurementParameters.lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
			    			debugLevel
			    			);
					sb.append("\t"+IJ.d2s(d,0));
					numSummedErrors++;
					sumErrors+=(d-x)*(d-x);

				}
				sb.append("\n");
			}
			// Add hysteresis if defined
			if (this.hysteresisHistoryFromTo!=null) {
				this.measuredHysteresis=0.0;
				for (int i=this.hysteresisHistoryFromTo[0];i<=this.hysteresisHistoryFromTo[1];i++)	{
					double x=getPositionAverage(i);
					String timestamp=getTimestamp(i);
					if (timestamp==null) timestamp="n/a";
					sb.append(i+"\t"+timestamp+"\t"+IJ.d2s(x,0));
					double [] sharp=this.history.get(i).getCenterResolutions();

					double l2=sharp[0]*sharp[0] + sharp[1]*sharp[1] + sharp[2]*sharp[2];
					sb.append("\t"+IJ.d2s(sharp[0],3)+"\t"+IJ.d2s(sharp[1],3)+"\t"+IJ.d2s(sharp[2],3)); // r,g,b here
					sb.append("\t"+IJ.d2s(sharp[0]*sharp[0]/l2,3));
					sb.append("\t"+IJ.d2s(sharp[2]*sharp[2]/l2,3));
					sb.append("\t"+IJ.d2s(l2,3));
					int j0=0;
					for (int j=0;j<this.lensDistanceCalibration.length;j++){
						if (this.lensDistanceCalibration[j][6]>x){
							j0=j;
							break;
						}
					}
					sb.append("\t"+j0);
					if (j0==0) sb.append("\t---\t---\t---\t---\t---\t---\t---");
					else {
						double k= (x-this.lensDistanceCalibration[j0-1][6])/(this.lensDistanceCalibration[j0][6]-this.lensDistanceCalibration[j0-1][6]);
						for (int j=0;j<3;j++){
							double d=k*this.lensDistanceCalibration[j0][j]+(1.0-k)*this.lensDistanceCalibration[j0-1][j];
							sb.append("\t"+IJ.d2s(d,3));
						}
						double d=getLensDistance(
								sharp, // {R-sharpness,G-sharpness,B-sharpness}
								false, // boolean absolute, // return absolutely calibrated data
1861
								focusMeasurementParameters.lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
1862 1863 1864 1865 1866 1867 1868 1869 1870 1871 1872 1873 1874 1875 1876 1877 1878 1879
								focusMeasurementParameters.lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
								debugLevel
						);
						sb.append("\t"+IJ.d2s(d,0));
						if (this.measuredHysteresis<Math.abs((d-x))) this.measuredHysteresis=Math.abs((d-x));
					}
					sb.append("\n");
				}
			}
			double rms=Double.NaN;
			if (numSummedErrors>0) rms=Math.sqrt(sumErrors/numSummedErrors);
			String title= "Calibrating Lens Distance: RMS="+IJ.d2s(rms,3)+((this.hysteresisHistoryFromTo!=null)?(" hysteresis="+this.measuredHysteresis):"");

    		if (focusMeasurementParameters.showResults && focusMeasurementParameters.lensDistanceShowResults) new TextWindow(title, header, sb.toString(),1000,500);
    		//micronsPerStep
    		int [] optimalMotorPosition=absoluteCalibrate(
        			micronsPerStep,
        			focusMeasurementParameters.targetMicrons,
1880
        			focusMeasurementParameters.lensDistanceWeightY,
Andrey Filippov's avatar
Andrey Filippov committed
1881 1882 1883 1884 1885 1886
        			focusMeasurementParameters.lensDistanceWeightK,
        			debugLevel);
			double [] limits=absoluteCalibrationLimits() ; // here we know that calibration was made
			String footer=null;
    		if (optimalMotorPosition!=null) {
    			footer="Lens is calibrated from "+IJ.d2s(limits[0],3)+" microns to "+IJ.d2s(limits[1],3)+" microns from the \"best\" position\nBest offsets per color:";
1887
// Add positions of color best focus
Andrey Filippov's avatar
Andrey Filippov committed
1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913
    			double [] bestColors=colorFocusDistances();
    			String [] colorNames={"red","green","blue"};
    			for (int c=0;c< bestColors.length;c++) footer+=" "+colorNames[c]+"="+(Double.isNaN(bestColors[c])?"out of range":(IJ.d2s(bestColors[c],2)+"um"));
    		} else {
    			footer="Scanned interval does not include best focus position.";
    			if (limits[0]==0.0) footer+=" Focus seems to be closer (smaller motor step values).";
    			if (limits[1]==0.0) footer+=" Focus seems to be farther (higher motor step values).";
    		}
    		if (this.hysteresisHistoryFromTo!=null){
    			footer+="Measured (maximal) motor hysteresis is "+IJ.d2s(this.measuredHysteresis,0)+" motor steps ("+IJ.d2s(micronsPerStep*this.measuredHysteresis,1)+" microns)";
    		}
    		footer+="\nComment: "+focusMeasurementParameters.comment;
			if (debugLevel>0) System.out.println(footer);
			if (focusMeasurementParameters.lensDistanceShowResults) IJ.showMessage("Success",footer);
    		// colorFocusDistances()
			if (focusMeasurementParameters.saveResults && (path!=null)){
    			CalibrationFileManagement.saveStringToFile (
						path,
						header+"\n"+sb.toString()+"\n"+footer);

			}
    		return optimalMotorPosition;
    	}

    	public double getLensDistance(
    			boolean absolute, // return absolutely calibrated data
1914
    			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
1915 1916 1917 1918 1919 1920
    			double weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			int debugLevel
    			){
    		return getLensDistance(
    				getCenterResolutions(), // {R-sharpness,G-sharpness,B-sharpness}
        			absolute, // return absolutely calibrated data
1921
        			weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
1922 1923 1924
        			weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
        			debugLevel
        			);
1925
    	}
Andrey Filippov's avatar
Andrey Filippov committed
1926 1927 1928 1929

    	public double getLensDistance(
    			int index,
    			boolean absolute, // return absolutely calibrated data
1930
    			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
1931 1932 1933 1934 1935 1936
    			double weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			int debugLevel
    			){
    		return getLensDistance(
    				getCenterResolutions(index), // {R-sharpness,G-sharpness,B-sharpness}
        			absolute, // return absolutely calibrated data
1937
        			weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
1938 1939 1940
        			weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
        			debugLevel
        			);
1941 1942
    	}

Andrey Filippov's avatar
Andrey Filippov committed
1943 1944 1945
    	public double getLensDistance(
    			double [] sharpVector, // {R-sharpness,G-sharpness,B-sharpness}
    			boolean absolute, // return absolutely calibrated data
1946
    			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
    			double weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			int debugLevel
    			){
 //   		System.out.println("getLensDistance(): debugLevel="+debugLevel);
 //   		System.out.println("getLensDistance(): weightK="+weightK);
 //   		System.out.println("getLensDistance(): weightY="+weightY);

    		int debugThreshold=2;
    		boolean useLogY=(weightY<0);
    		if (useLogY) weightY=-weightY;
    		// Caclulate yMmax
    		if (this.lensDistanceCalibration==null) return Double.NaN;
    		if (sharpVector==null) return Double.NaN;
    		double yMax=0.0;
    		for (int i=0;i<this.lensDistanceCalibration.length;i++) if (this.lensDistanceCalibration[i][2]>yMax)yMax=this.lensDistanceCalibration[i][2];
    		if (useLogY) yMax=Math.log(yMax);
    		double wY=weightY/yMax;
    		double best=Double.NaN;
    		int iBest=0;
    		double [] weights=new double[3];
    		double [] sample=new double[3];
			double s=0.0;
			for (int j=0;j<3;j++){
				sample[j]=sharpVector[j]*sharpVector[j];
				s+= sample[j];
			}
			sample[0]/=s;
			sample[1]=sample[2]/s;
			sample[2]=useLogY?Math.log(s):s;
if (debugLevel>=debugThreshold) System.out.println("sharpVector= {"+sharpVector[0]+" "+sharpVector[1]+" "+sharpVector[2]+"}, weightK="+weightK+" weightY="+weightY);
if (debugLevel>=debugThreshold) System.out.println("Sample="+sample[0]+" "+sample[1]+" "+sample[2]+" yMax="+yMax);
    		double [] diff=new double[3];
    		for (int i=0;i<this.lensDistanceCalibration.length;i++){
    			weights[0]=this.lensDistanceCalibration[i][3];
    			weights[1]=this.lensDistanceCalibration[i][4];
    			weights[2]=this.lensDistanceCalibration[i][5]*wY;
    			double w=weights[0]+weights[1]+weights[2];
    			for (int j=0;j<3;j++) weights[j]=weightK*weights[j]/w+(1.0-weightK)/3.0;
    			diff[0]=sample[0]-this.lensDistanceCalibration[i][0];
    			diff[1]=sample[1]-this.lensDistanceCalibration[i][1]; // sample[2], not 1 here (blue)
    			diff[2]=wY*(sample[2]-(useLogY?Math.log(this.lensDistanceCalibration[i][2]):this.lensDistanceCalibration[i][2]));

    			double err=0.0;
    			for (int j=0;j<3;j++) err+=weights[j]*diff[j]*diff[j];
    			if ((i==0) || (err<best)){
    				best=err;
    				iBest=i;
    			}
if (debugLevel>=debugThreshold) System.out.println(i+" "+diff[0]+" "+diff[1]+" "+diff[2]+" "+weights[0]+" "+weights[1]+" "+weights[2]+" "+err);
    		}
    		return this.lensDistanceCalibration[iBest][absolute?7:6];
    	}
    	/**
    	 * Applies absolute (from the "best") calibration to the distance between the lens and the center of the sensor
    	 * after the calibration ion motor steps for the range that includes "best focus" was performed.
    	 * "Best" here is the maximal sum of reversed areas of the focal spots (50% intensity) for each of the 3 colors,
    	 * measured for the "center" samples. After the calibration is performed, array lensDistanceCalibration can be used
    	 * to find the distance from the lens to the sensor center knowing the focal spot sizes for each of 3 colors.
2005
    	 *
Andrey Filippov's avatar
Andrey Filippov committed
2006 2007 2008 2009 2010 2011 2012 2013
    	 * @param micronsPerStep lens movement per step (all 3 motors together)
    	 * @param debugLevel - if >1, will produce debug output
    	 * @return motor position corresponding to zero shift ("best focus") or null if it is not inside
    	 * the scanned interval. May use absoluteCalibrationLimits() to find out which of the limits is wrong.
    	 */
    	public int[] absoluteCalibrate(
    			double micronsPerStep,
    			double targetMicrons,
2014
    			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2015 2016 2017
    			double weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			int debugLevel){
// Temporary!
2018
    		debugLevel=debugLevel+0;
Andrey Filippov's avatar
Andrey Filippov committed
2019 2020 2021 2022 2023 2024
    		this.micronsPerStep=micronsPerStep;
    		if (this.lensDistanceCalibration==null) {
    			String msg="Lens travel is not yet measured/calibrated, impossible to apply absolute calibration";
    			IJ.showMessage("Error",msg);
    			throw new IllegalArgumentException (msg);
    		}
2025

Andrey Filippov's avatar
Andrey Filippov committed
2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037
    		double bestY=0.0;
    		int i0=0;
    		for (int i=0;i<this.lensDistanceCalibration.length;i++){
    			if (this.lensDistanceCalibration[i][this.indexY]>bestY){
    				bestY=this.lensDistanceCalibration[i][this.indexY];
    				i0=i;
    			}
    		}
    		// use local maximum if the curve has N- shape
    		if ((i0==0) || (i0>=(this.lensDistanceCalibration.length-1))){
    			int i1=0;
    			int i2=this.lensDistanceCalibration.length-1;
2038
    			while ((i1<(i2-1)) && (this.lensDistanceCalibration[i1+1][this.indexY]<this.lensDistanceCalibration[i1][this.indexY])) i1++;
Andrey Filippov's avatar
Andrey Filippov committed
2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078
    			while ((i2>(i1+1)) && (this.lensDistanceCalibration[i2-1][this.indexY]<this.lensDistanceCalibration[i2][this.indexY])) i2--;
    			if (i1<(i2-1)){
    				bestY=0.0;
    				for (int i=i1;i<i2;i++){
    					if (this.lensDistanceCalibration[i][this.indexY]>bestY){
    						bestY=this.lensDistanceCalibration[i][this.indexY];
    						i0=i;
    					}
    				}
    			}
  				System.out.println("absoluteCalibrate(): ***** Warning: Absolute maximum of approximated sharpness function is not reached, using local maximum at "+i0+"/"+this.lensDistanceCalibration.length);
    		}
    		double x0=i0;
    		if ((i0>0) && (i0<(this.lensDistanceCalibration.length-1))){
    			double [][] data= {
    					{x0-1,this.lensDistanceCalibration[i0-1][this.indexY]},
    					{x0,  this.lensDistanceCalibration[i0  ][this.indexY]},
    					{x0+1,this.lensDistanceCalibration[i0+1][this.indexY]}
    			};
    			double [] coeff=(new PolynomialApproximation()).polynomialApproximation1d(data, 2);
    			if ((coeff!=null) && (coeff[2]<0.0)) {
    				x0=-coeff[1]/(2.0*coeff[2]);
    			}
    		}
			if (debugLevel>1) System.out.println ("absoluteCalibrate(): i0="+i0+", interpolated(?)  x0="+IJ.d2s(x0,3)+" motor average="+this.lensDistanceCalibration[i0][this.indexSteps]);
			double micronsPerElement=
				(this.lensDistanceCalibration[this.lensDistanceCalibration.length-1][indexSteps]-this.lensDistanceCalibration[0][indexSteps])/
				(this.lensDistanceCalibration.length-1)*this.micronsPerStep;
    		for (int i=0;i<this.lensDistanceCalibration.length;i++) {
    			this.lensDistanceCalibration[i][this.indexMicrons]= (i-x0)*micronsPerElement;
    			if (debugLevel>1) System.out.println ("absoluteCalibrate(): this.lensDistanceCalibration["+i+"]["+this.indexMicrons+"]="+this.lensDistanceCalibration[i][this.indexMicrons]);
    		}
//    		int [] position= parallelPositionMove(targetMicrons, weightK, weightY, debugLevel); // "optimal" motor position (or null)
    		int [] position= parallelPositionMove(0.0, weightK, weightY, debugLevel); // "optimal" motor position (or null) // here - disregard targetMicrons
    		if (debugLevel>0) if (position!=null) System.out.println("Optimal focus in the center at position ["+position[0]+","+position[1]+","+position[2]+"]");
    		else System.out.println("Optimal focus is undefined");
    		return position;
    	}

    	/**
2079
    	 *
Andrey Filippov's avatar
Andrey Filippov committed
2080 2081 2082 2083 2084 2085
    	 * @param targetMicrons requested distance from the lens to the center, 0 is the "best focus" (between green and red), in microns
    	 * @param weightK scale between derivative-driven (1.0) and equal (0.0) weights in fitting measured parameters to find distance
    	 * @param weightY relative weight of overall (not color-differential) "sharpness" in fitting
    	 * @param debugLevel
    	 * @return new motors position
    	 */
2086

Andrey Filippov's avatar
Andrey Filippov committed
2087 2088
    	public int [] parallelPositionMove(
    			double targetMicrons,
2089
    			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2090 2091 2092 2093
    			double weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			int debugLevel
    			){
    		if (this.lensDistanceCalibration==null){
2094
    			if (debugLevel>0) System.out.println("parallelPositionMove(): this.lensDistanceCalibration==null");
Andrey Filippov's avatar
Andrey Filippov committed
2095 2096 2097 2098 2099 2100
    			return null;// not yet calibrated
    		}
    		int last=this.lensDistanceCalibration.length-1;

    		double [] limits=absoluteCalibrationLimits();
    		if ((limits==null) || ((limits[0]*limits[1])>0) || ((limits[0]==0) && (limits[1]==0))){
2101
    			if ((debugLevel>0) && (limits==null)) System.out.println("parallelPositionMove(): limits==null");
Andrey Filippov's avatar
Andrey Filippov committed
2102
    			if ((debugLevel>0) && ((limits[0]*limits[1])>=0)){
2103
    				System.out.println("parallelPositionMove(): ((limits[0]*limits[1])>=0), limits[0]="+limits[0]+", limits[1]="+limits[1]);
Andrey Filippov's avatar
Andrey Filippov committed
2104 2105 2106 2107 2108 2109 2110 2111 2112
    			}
    			return null;
    		}

    		if (
    				(targetMicrons<this.lensDistanceCalibration[0][indexMicrons]) ||
    				(targetMicrons>this.lensDistanceCalibration[last][indexMicrons])) {
				System.out.println("parallelPositionMove(): out of range, this.lensDistanceCalibration[0]["+indexMicrons+"]="+
						this.lensDistanceCalibration[0][indexMicrons]+
2113
						", this.lensDistanceCalibration["+last+"]["+indexMicrons+"]="+this.lensDistanceCalibration[last][indexMicrons]);
Andrey Filippov's avatar
Andrey Filippov committed
2114 2115 2116 2117 2118
    			return null ; // out of range
    		}
    		double x=
    			(targetMicrons - this.lensDistanceCalibration[0][indexMicrons])/
    			(this.lensDistanceCalibration[last][indexMicrons]-this.lensDistanceCalibration[0][indexMicrons]);
2119
    		if (debugLevel>1) System.out.println("parallelPositionMove(): x="+x);
Andrey Filippov's avatar
Andrey Filippov committed
2120 2121 2122
    		double newAveragePos=
    			this.lensDistanceCalibration[0][this.indexSteps]+
    			x*(this.lensDistanceCalibration[last][indexSteps]-this.lensDistanceCalibration[0][indexSteps]);
2123
    		if (debugLevel>1) System.out.println("parallelPositionMove(): newAveragePos="+newAveragePos);
Andrey Filippov's avatar
Andrey Filippov committed
2124 2125
			int [] position=getPosition(); // last position from the history ** here was the bug (made returning clone()
//			double [] resolutions=getCenterResolutions();
2126

Andrey Filippov's avatar
Andrey Filippov committed
2127 2128 2129
			double predictedAverage=getLensDistance(
					getCenterResolutions(), // {R-sharpness,G-sharpness,B-sharpness}
	    			false, //true, // return absolutely calibrated data
2130
	    			weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2131 2132 2133 2134
	    			weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
	    			debugLevel);
			if (Double.isNaN(predictedAverage)) return null;
			int move = (int) (newAveragePos-predictedAverage);
2135
    		if (debugLevel>1) System.out.println("parallelPositionMove(): calculated move of the center is "+move+" steps of 3 motors in parallel.");
Andrey Filippov's avatar
Andrey Filippov committed
2136 2137
			for (int i=0;i<3;i++) position[i]+=move;
    		return position; // "optimal" motor position
2138

Andrey Filippov's avatar
Andrey Filippov committed
2139
    	}
2140

Andrey Filippov's avatar
Andrey Filippov committed
2141
    	/**
2142
    	 * Calculates offsets (in microns) of each color best focus (if it is inside range)
Andrey Filippov's avatar
Andrey Filippov committed
2143 2144 2145 2146 2147 2148 2149 2150 2151 2152 2153 2154 2155 2156 2157 2158 2159 2160 2161 2162 2163 2164 2165 2166 2167 2168 2169 2170 2171 2172 2173 2174 2175 2176 2177 2178 2179 2180
    	 * @return
    	 */
    	public double [] colorFocusDistances(){
    		int [] bestIndices={0,0,0};
    		double [] bestResolution={0.0,0.0,0.0};
    		for (int i=0;i<this.lensDistanceCalibration.length;i++) {
    			double l2=    this.lensDistanceCalibration[i][this.indexY];
    			double rFrac2=this.lensDistanceCalibration[i][this.indexR];
    			double bFrac2=this.lensDistanceCalibration[i][this.indexB];
    			double [] resolutions={
    					Math.sqrt(rFrac2*l2),
    					Math.sqrt((1-rFrac2-bFrac2)*l2),
    					Math.sqrt(bFrac2*l2)};
    			for (int c=0;c< resolutions.length;c++) if (resolutions[c]>bestResolution[c]){
    				bestResolution[c]=resolutions[c];
    				bestIndices[c]=i;
    			}
    		}
    		double [] bestFocus ={Double.NaN,Double.NaN,Double.NaN}; // maximum outside of range
    		for (int c=0;c<bestResolution.length;c++) if ((bestIndices[c]>0) && (bestIndices[c]<(this.lensDistanceCalibration.length-1))) {
    			double [][] data= {
    					{bestIndices[c]-1,0.0},
    					{bestIndices[c],  0.0},
    					{bestIndices[c]+1,0.0}
    			};
    			for (int j=-1;j<2;j++){
        			double l2=    this.lensDistanceCalibration[bestIndices[c]+j][this.indexY];
        			double rFrac2=this.lensDistanceCalibration[bestIndices[c]+j][this.indexR];
        			double bFrac2=this.lensDistanceCalibration[bestIndices[c]+j][this.indexB];
        			double [] resolutions={
        					Math.sqrt(rFrac2*l2),
        					Math.sqrt((1-rFrac2-bFrac2)*l2),
        					Math.sqrt(bFrac2*l2)};
        			data[j+1][1]=resolutions[c];
    			}
    			double [] coeff=(new PolynomialApproximation()).polynomialApproximation1d(data, 2);
    			if ((coeff!=null) && (coeff[2]<0.0)) {
    				double x=(-coeff[1]/(2.0*coeff[2]))-bestIndices[c];
2181

Andrey Filippov's avatar
Andrey Filippov committed
2182 2183 2184 2185 2186 2187
    				bestFocus[c]= this.lensDistanceCalibration[bestIndices[c]][this.indexMicrons]+
    				x*(this.lensDistanceCalibration[bestIndices[c]+1][this.indexMicrons]-this.lensDistanceCalibration[bestIndices[c]][this.indexMicrons]);
    			}
    		}
    		return bestFocus;
    	}
2188

Andrey Filippov's avatar
Andrey Filippov committed
2189
     	/**
2190
    	 *
Andrey Filippov's avatar
Andrey Filippov committed
2191
    	 * @return array of lower and upper limits of the calibration range. None should be 0.0, that would meant that maximum
2192
    	 * was not covered by the scanned interval Returns null if there is no calibration
Andrey Filippov's avatar
Andrey Filippov committed
2193 2194 2195 2196 2197 2198 2199 2200 2201 2202
    	 */
    	public double [] absoluteCalibrationLimits(){
    		if (this.lensDistanceCalibration==null)  return null; // no calibration at all
    		if (Double.isNaN(this.lensDistanceCalibration[0][this.indexMicrons])) return null; // no absolute calibration
    		double [] limits={
    				this.lensDistanceCalibration[0][this.indexMicrons],
    				this.lensDistanceCalibration[this.lensDistanceCalibration.length-1][this.indexMicrons]
    		};
    		return limits;

2203 2204
    	}
/*
Andrey Filippov's avatar
Andrey Filippov committed
2205 2206 2207 2208 2209 2210 2211 2212 2213 2214 2215 2216 2217 2218 2219 2220 2221 2222 2223 2224 2225 2226 2227 2228 2229 2230 2231 2232 2233 2234 2235 2236 2237 2238 2239 2240 2241 2242 2243 2244 2245 2246 2247 2248 2249 2250 2251 2252 2253 2254 2255
    	public void calibrateFractions( // positive away from the lens!
    			int fromHistory, // first history data to consider
    			int toHistory,   // last history data to consider
    			int numPoints,
    			int polyDegree,  // polynomial degree to fit
//    			boolean sameTiltOnly,
    			boolean centerOnly, // only use center samples
    			double xMin,  // range of average motor coordinates to consider (should have the same tilt as last in history)
    			double xMax,
    			int debugLevel
    			){
    		if (fromHistory<0)fromHistory=0;          // already done this by the caller
    		if (toHistory>=size())toHistory=size()-1;

    		int indexR50=centerOnly?6:3; // 6 is average R50 in the center samples,3 - in all samples
    		int sameTiltSize=0;
//    		int [] origin=getPosition();
    		int [] origin=getPosition(toHistory);
    		boolean []sameTiltmask=new boolean[size()];
    		for (int i=0;i<sameTiltmask.length;i++) sameTiltmask[i]=false;
    		for (int i=fromHistory;i<=toHistory;i++){
    			sameTiltmask[i]=true;
    			sameTiltSize++;
    			int [] position=getPosition(i);
    			if (((position[0]-position[1])==(origin[0]-origin[1])) && ((position[2]-position[1])==(origin[2]-origin[1]))){
    				sameTiltmask[i]=true;
    				sameTiltSize++;
    			}
    		}

//    		System.out.println("solveSinglePoly(): sameTiltSize="+sameTiltSize+", size()="+size()+" debugLevel="+debugLevel);
    		if (debugLevel>1) System.out.println("solveSinglePoly(): sameTiltSize="+sameTiltSize+", size()="+size()+
    				" fromHistory="+fromHistory+" toHistory="+toHistory+" numPoints="+numPoints);
//    		double x0=(origin[0]+origin[1]+origin[2])/3.0;
			double x0=(xMax+xMin)/2;// center
			double scale=2.0/(xMax-xMin);
    		double [][] polyCoeff =new double[3][];
			PolynomialApproximation pa= new PolynomialApproximation(debugLevel);
    		for (int color=0;color<3;color++){
    			double [][] data=new double [size()][3]; //x,f,w
    			for (int i=0;i<size();i++) {
    				data[i][2]=0.0; // weight
    				if (sameTiltmask[i]){
    					double x=getPositionAverage(i);
    					if ((x<xMin) || (x>xMax)) continue; // out of range point
    					data[i][0]= scale*(x - x0); // interval center
    					double [][] metrics= this.history.get(i).getMetrics(0.0,0.0,0.0); // average is not used, any scales
    					double l2=
    						1.0/metrics[0][indexR50]/metrics[0][indexR50]+
    						1.0/metrics[1][indexR50]/metrics[1][indexR50]+
    						1.0/metrics[2][indexR50]/metrics[2][indexR50];
2256
    					if      (color==0) data[i][1]=1.0/metrics[0][indexR50]/metrics[0][indexR50]/l2; // red
Andrey Filippov's avatar
Andrey Filippov committed
2257 2258 2259 2260 2261 2262 2263 2264 2265 2266 2267 2268 2269 2270 2271 2272 2273 2274
    					else if (color==1) data[i][1]=1.0/metrics[2][indexR50]/metrics[2][indexR50]/l2; // blue
    					else               data[i][1]=l2; // Y-combined
 //   					data[i][1]=1.0/metrics[color][indexR50];// "sharpness"
    					data[i][2]=1.0;
    					if (debugLevel>1)	System.out.println("color=\t"+color+"\ti=\t"+i+"\tdata[i][0]=\t"+data[i][0]+"\tdata[i][1]=\t"+data[i][1]+"\tdata[i][2]=\t"+data[i][2]);
    				}
    			}
    			polyCoeff[color]=pa.polynomialApproximation1d(data, polyDegree);
    			if (debugLevel>1) {
    				String dbgStr="polyCoeff["+color+"]=";
    				for (int i=0;i<polyCoeff[color].length;i++) {
    					dbgStr+="\t"+polyCoeff[color][i];
    				}
    				System.out.println(dbgStr);
    			}
    		}
    		this.lensDistanceCalibration=new double [numPoints][8]; // fracR, fracG, fracB, weightR, weightG, weightB, centerSteps, microns (not assigned here)
			for (int n=0;n<numPoints;n++) {
2275
				this.lensDistanceCalibration[n][6]=xMin+n*(xMax-xMin)/(numPoints-1);
Andrey Filippov's avatar
Andrey Filippov committed
2276 2277 2278 2279 2280 2281 2282 2283 2284 2285 2286 2287 2288 2289 2290 2291 2292
				double x=(2.0*n)/(numPoints-1)-1.0;
				double [] xp=new double[polyDegree+1];
				xp[0]=1;
				for (int i=1;i<=polyDegree;i++) xp[i]=xp[i-1]*x;
//				double sumDerSq=0.0;
				for (int color=0;color<3;color++){
					this.lensDistanceCalibration[n][color]=0;
					this.lensDistanceCalibration[n][color+3]=0; // derivatives^2
					for (int i=0;i<=polyDegree;i++){
						this.lensDistanceCalibration[n][color]+=polyCoeff[color][i]*xp[i];
					}
					for (int i=1;i<=polyDegree;i++){
						this.lensDistanceCalibration[n][color+3]+=polyCoeff[color][i]*xp[i-1]*i; // derivative
					}
					this.lensDistanceCalibration[n][color+3]*=this.lensDistanceCalibration[n][color+3]; //squared
//					sumDerSq+=this.lensDistanceCalibration[n][color+3];
				}
2293
//				for (int color=0;color<3;color++)this.lensDistanceCalibration[n][color+3] /= sumDerSq;// should be non-zero, I hope
Andrey Filippov's avatar
Andrey Filippov committed
2294 2295 2296 2297 2298 2299 2300 2301 2302 2303
				this.lensDistanceCalibration[n][7]=Double.NaN; // not yet absolutely calibrated
				if (debugLevel>1) System.out.println(n+":\t"+
						IJ.d2s(this.lensDistanceCalibration[n][0],4)+"\t"+
						IJ.d2s(this.lensDistanceCalibration[n][1],4)+"\t"+
						IJ.d2s(this.lensDistanceCalibration[n][2],4)+"\t"+
						IJ.d2s(this.lensDistanceCalibration[n][3],4)+"\t"+
						IJ.d2s(this.lensDistanceCalibration[n][4],4)+"\t"+
						IJ.d2s(this.lensDistanceCalibration[n][5],4)+"\t"+
						IJ.d2s(this.lensDistanceCalibration[n][6],4));
			}
2304 2305


Andrey Filippov's avatar
Andrey Filippov committed
2306 2307 2308 2309
			//    	private double lensDistanceCalibrationYScale=0.0; // R-frac, B-frac, Y*scale

//    		return calib; // should be stored between calls
    	}
2310 2311 2312
 */


Andrey Filippov's avatar
Andrey Filippov committed
2313 2314 2315 2316 2317 2318 2319 2320 2321 2322 2323 2324 2325 2326 2327 2328 2329 2330 2331 2332 2333 2334 2335 2336 2337 2338 2339 2340 2341 2342 2343 2344 2345 2346 2347 2348 2349 2350 2351 2352 2353 2354
    	public double [][] calibrateFractions( // positive away from the lens!
    			int fromHistory, // first history data to consider
    			int toHistory,   // last history data to consider
//    			int numPoints,
    			int polyDegree,  // polynomial degree to fit
    			boolean centerOnly, // only use center samples
    			double xMin,  // range of average motor coordinates to consider (should have the same tilt as last in history)
    			double xMax,
    			int debugLevel
    			){
    		if (fromHistory<0)fromHistory=0;          // already done this by the caller
    		if (toHistory>=size())toHistory=size()-1;

    		int indexR50=centerOnly?6:3; // 6 is average R50 in the center samples,3 - in all samples
    		int sameTiltSize=0;
    		int [] origin=getPosition(toHistory);
    		boolean []sameTiltmask=new boolean[size()];
    		for (int i=0;i<sameTiltmask.length;i++) sameTiltmask[i]=false;
    		for (int i=fromHistory;i<=toHistory;i++){
    			sameTiltmask[i]=true;
    			sameTiltSize++;
    			int [] position=getPosition(i);
    			if (((position[0]-position[1])==(origin[0]-origin[1])) && ((position[2]-position[1])==(origin[2]-origin[1]))){
    				sameTiltmask[i]=true;
    				sameTiltSize++;
    			}
    		}
    		if (debugLevel>1) System.out.println("solveSinglePoly(): sameTiltSize="+sameTiltSize+", size()="+size()+
    				" fromHistory="+fromHistory+" toHistory="+toHistory);
			double x0=(xMax+xMin)/2;// center
			double scale=2.0/(xMax-xMin);
    		double [][] polyCoeff =new double[3][];
			PolynomialApproximation pa= new PolynomialApproximation(debugLevel);
    		for (int color=0;color<3;color++){
    			double [][] data=new double [size()][3]; //x,f,w
    			for (int i=0;i<size();i++) {
    				data[i][2]=0.0; // weight
    				if (sameTiltmask[i]){
    					double x=getPositionAverage(i);
    					if ((x<xMin) || (x>xMax)) continue; // out of range point
    					data[i][0]= scale*(x - x0); // interval center
    					double [][] metrics= this.history.get(i).getMetrics(0.0,0.0,0.0); // average is not used, any scales
2355
    					if (metrics==null) continue;
Andrey Filippov's avatar
Andrey Filippov committed
2356 2357 2358 2359
    					double l2=
    						1.0/metrics[0][indexR50]/metrics[0][indexR50]+
    						1.0/metrics[1][indexR50]/metrics[1][indexR50]+
    						1.0/metrics[2][indexR50]/metrics[2][indexR50];
2360
    					if      (color==0) data[i][1]=1.0/metrics[0][indexR50]/metrics[0][indexR50]/l2; // red
Andrey Filippov's avatar
Andrey Filippov committed
2361 2362 2363 2364 2365 2366 2367 2368 2369 2370 2371 2372 2373 2374 2375 2376 2377
    					else if (color==1) data[i][1]=1.0/metrics[2][indexR50]/metrics[2][indexR50]/l2; // blue
    					else               data[i][1]=l2; // Y-combined
 //   					data[i][1]=1.0/metrics[color][indexR50];// "sharpness"
    					data[i][2]=1.0;
    					if (debugLevel>1)	System.out.println("color=\t"+color+"\ti=\t"+i+"\tdata[i][0]=\t"+data[i][0]+"\tdata[i][1]=\t"+data[i][1]+"\tdata[i][2]=\t"+data[i][2]);
    				}
    			}
    			polyCoeff[color]=pa.polynomialApproximation1d(data, polyDegree);
    			if (debugLevel>1) {
    				String dbgStr="polyCoeff["+color+"]=";
    				for (int i=0;i<polyCoeff[color].length;i++) {
    					dbgStr+="\t"+polyCoeff[color][i];
    				}
    				System.out.println(dbgStr);
    			}
    		}
    		return polyCoeff;
2378 2379


Andrey Filippov's avatar
Andrey Filippov committed
2380 2381 2382 2383 2384 2385 2386 2387 2388 2389 2390 2391
			//    	private double lensDistanceCalibrationYScale=0.0; // R-frac, B-frac, Y*scale

//    		return calib; // should be stored between calls
    	}
    	public void makeCalibrationTable(
    			double [][] polyCoeff,
    			int numPoints,
//    			int polyDegree,  // polynomial degree to fit
    			double xMin,  // range of average motor coordinates to consider (should have the same tilt as last in history)
    			double xMax,
    			int debugLevel
    			){
2392
//    		debugLevel=2;
Andrey Filippov's avatar
Andrey Filippov committed
2393 2394 2395 2396 2397 2398 2399 2400 2401
   			int polyDegree=polyCoeff[0].length-1;  // polynomial degree
    		this.lensDistanceCalibration=new double [numPoints][8]; // fracR, fracG, fracB, weightR, weightG, weightB, centerSteps, microns (not assigned here)
    		if (debugLevel>1) {
    			for (int c=0;c<3;c++) {
    				System.out.println("makeCalibrationTable(), c="+c);
    				for (int n=0;n<polyDegree;n++) System.out.println(n+": "+polyCoeff[c][n]);
    			}
    		}
    		for (int n=0;n<numPoints;n++) {
2402
    			this.lensDistanceCalibration[n][6]=xMin+n*(xMax-xMin)/(numPoints-1);
Andrey Filippov's avatar
Andrey Filippov committed
2403 2404 2405 2406 2407 2408 2409 2410 2411 2412 2413 2414 2415 2416 2417 2418 2419
    			double x=(2.0*n)/(numPoints-1)-1.0;
    			double [] xp=new double[polyDegree+1];
    			xp[0]=1;
    			for (int i=1;i<=polyDegree;i++) xp[i]=xp[i-1]*x;
    			//		double sumDerSq=0.0;
    			for (int color=0;color<3;color++){
    				this.lensDistanceCalibration[n][color]=0;
    				this.lensDistanceCalibration[n][color+3]=0; // derivatives^2
    				for (int i=0;i<=polyDegree;i++){
    					this.lensDistanceCalibration[n][color]+=polyCoeff[color][i]*xp[i];
    				}
    				for (int i=1;i<=polyDegree;i++){
    					this.lensDistanceCalibration[n][color+3]+=polyCoeff[color][i]*xp[i-1]*i; // derivative
    				}
    				this.lensDistanceCalibration[n][color+3]*=this.lensDistanceCalibration[n][color+3]; //squared
    				//			sumDerSq+=this.lensDistanceCalibration[n][color+3];
    			}
2420
    			//		for (int color=0;color<3;color++)this.lensDistanceCalibration[n][color+3] /= sumDerSq;// should be non-zero, I hope
Andrey Filippov's avatar
Andrey Filippov committed
2421 2422 2423 2424 2425 2426 2427 2428 2429 2430 2431 2432 2433
    			this.lensDistanceCalibration[n][7]=Double.NaN; // not yet absolutely calibrated
    			if (debugLevel>1) System.out.println(n+":\t"+
    					IJ.d2s(this.lensDistanceCalibration[n][0],4)+"\t"+
    					IJ.d2s(this.lensDistanceCalibration[n][1],4)+"\t"+
    					IJ.d2s(this.lensDistanceCalibration[n][2],4)+"\t"+
    					IJ.d2s(this.lensDistanceCalibration[n][3],4)+"\t"+
    					IJ.d2s(this.lensDistanceCalibration[n][4],4)+"\t"+
    					IJ.d2s(this.lensDistanceCalibration[n][5],4)+"\t"+
    					IJ.d2s(this.lensDistanceCalibration[n][6],4));
    		}

    	}

2434

Andrey Filippov's avatar
Andrey Filippov committed
2435 2436 2437 2438 2439 2440 2441 2442 2443 2444 2445 2446 2447 2448 2449 2450 2451 2452 2453
    	/**
    	 *  Calculates worst ratio of the residual error to threshold. First parameters is focus or radial/tangential ratio, second and third - tilts
    	 *  Tilts should be 0, while the first is compared against "target"
    	 * @param weightRatioRedToGreen
    	 * @param weightRatioBlueToGreen
    	 * @param lensDistanceWeightK
    	 * @param lensDistanceWeightY
    	 * @param useRadialTangential
    	 * @param targetFarNear
    	 * @param targetMicrons
    	 * @param toleranceMicrons
    	 * @param toleranceTilt
    	 * @param toleranceThreshold
    	 * @param debugLevel
    	 * @return
    	 */
     	public double worstOverTolerance (
    			double weightRatioRedToGreen,
    			double weightRatioBlueToGreen,
2454
    			double lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2455 2456 2457 2458 2459 2460 2461 2462
    			double lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
				double targetMicrons,
				double toleranceMicrons,
				double toleranceTilt,
    			int debugLevel
    	){ return worstOverTolerance (
    			weightRatioRedToGreen,
    			weightRatioBlueToGreen,
2463
    			lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2464 2465 2466 2467 2468 2469 2470 2471 2472 2473 2474
    			lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
				false, //boolean useRadialTangential,
    			0.0,   //double targetFarNear,
				targetMicrons,
				toleranceMicrons,
				toleranceTilt,
    			debugLevel);
     	}
    	public double worstOverTolerance (
    			double weightRatioRedToGreen,
    			double weightRatioBlueToGreen,
2475
    			double lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2476 2477 2478 2479 2480 2481 2482 2483 2484 2485 2486 2487 2488 2489 2490 2491 2492 2493 2494 2495 2496
    			double lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
				boolean useRadialTangential,
    			double targetFarNear,
				double targetMicrons,
				double toleranceMicrons,
				double toleranceTilt,
    			int debugLevel
    	){
    		if (this.history.size()<1)  return Double.NaN;
    		if (!useRadialTangential){
    			// Make sure calibration is done and targetMicrons is within limits
    			double [] limits=absoluteCalibrationLimits();
    			if ((limits==null) || Double.isNaN(limits[0]) || Double.isNaN(limits[1])) return Double.NaN; // nonexistent/invalid calibration
    		}
    		double targetValue=useRadialTangential?targetFarNear:targetMicrons;
    		double [] lastAverageMetrics=this.history.get(this.history.size()-1).getMetrics(weightRatioRedToGreen,weightRatioBlueToGreen)[3];
    		double [] lastResolutions= this.history.get(this.history.size()-1).getCenterResolutions();
    		if (!useRadialTangential){
    			lastAverageMetrics[0]=getLensDistance(
    					lastResolutions, // {R-sharpness,G-sharpness,B-sharpness}
		    			true, // boolean absolute, // return absolutely calibrated data
2497
		    			lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2498 2499 2500 2501 2502 2503 2504 2505 2506 2507
		    			lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
		    			1 //debugLevel
		    			);
    		}
    		lastAverageMetrics[0]-=targetValue;
    		lastAverageMetrics[0]/=useRadialTangential?toleranceTilt:toleranceMicrons; //
    		lastAverageMetrics[1]/=toleranceTilt;
    		lastAverageMetrics[2]/=toleranceTilt;
    		return Math.max(Math.abs(lastAverageMetrics[0]),Math.max(Math.abs(lastAverageMetrics[1]),Math.abs(lastAverageMetrics[2])));
     	}
2508

Andrey Filippov's avatar
Andrey Filippov committed
2509 2510 2511
     	public int [] focusReadjustStep(
     			double targetMicrons, // target focal distance
     			double micronsFade, // reduce influence of far points
2512
    			double lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2513 2514 2515 2516 2517 2518 2519 2520 2521 2522 2523 2524 2525 2526 2527 2528 2529 2530 2531 2532 2533 2534 2535 2536 2537 2538 2539 2540
    			double lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			int debugLevel
     			){
			double [] limits=absoluteCalibrationLimits();
    		int last=this.history.size()-1;
			if (limits==null){
    			String msg="Failed: lens focus quality per colors vs. distance is not yet measured/calibrated";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
			}
			if (Double.isNaN(limits[0]) || Double.isNaN(limits[1])){
    			String msg="Failed: lens focus quality per colors vs. distance is measured, but absolute calibration was not performed.";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
			}
			if ((targetMicrons<limits[0]) || (targetMicrons>limits[1])){
    			String msg="Failed: requested distance ("+targetMicrons+"um) is not withing calibrated range "+limits[0]+".."+limits[1]+"um.";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
			}
			if (this.history.size()<1){
    			String msg="Failed: history is empty.";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
2541

Andrey Filippov's avatar
Andrey Filippov committed
2542 2543 2544 2545 2546 2547
			}
    		if (
    				(targetMicrons<this.lensDistanceCalibration[0][indexMicrons]) ||
    				(targetMicrons>this.lensDistanceCalibration[this.lensDistanceCalibration.length-1][indexMicrons])) {
				System.out.println("parallelPositionMove(): out of range, this.lensDistanceCalibration[0]["+indexMicrons+"]="+
						this.lensDistanceCalibration[0][indexMicrons]+
2548
						", this.lensDistanceCalibration["+(this.lensDistanceCalibration.length-1)+"]["+indexMicrons+"]="+this.lensDistanceCalibration[this.lensDistanceCalibration.length-1][indexMicrons]);
Andrey Filippov's avatar
Andrey Filippov committed
2549 2550 2551 2552 2553 2554 2555 2556 2557 2558 2559 2560 2561 2562 2563 2564 2565 2566 2567 2568 2569
    			return null ; // out of range
    		}
			if (debugLevel>1) {
				System.out.println("focusReadjustStep(), last="+last);
			}

    		int [] lastPosition=getPosition();
    		// use only same difference between motors
    		int [] diffs= {lastPosition[0]-lastPosition[2],lastPosition[1]-lastPosition[2]};
			double [][] data = new double[this.history.size()][3];
    		double k=(micronsFade>0.0)?0.5/micronsFade/micronsFade:1.0;
    		int numSameDiff=0;
    		double [] resolutions;
			for (int nSample=0;nSample<this.history.size();nSample++){
				int [] motors=       getPosition(nSample);
				if ((motors[0]-motors[2]==diffs[0]) && (motors[1]-motors[2]==diffs[1])){
					resolutions= this.history.get(nSample).getCenterResolutions();
					data [nSample][0]=(motors[0]+motors[1]+motors[2])/3.0;
					data [nSample][1]=getLensDistance(
							resolutions, // {R-sharpness,G-sharpness,B-sharpness}
							true, // boolean absolute, // return absolutely calibrated data
2570
							lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2571 2572 2573 2574 2575 2576 2577 2578 2579 2580 2581 2582
							lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
							debugLevel
					);
					if (micronsFade<=0) data [nSample][2]=1.0;
					else {
						double diffDistance=data [nSample][1]-targetMicrons;
						data [nSample][2]=Math.exp(-k*diffDistance*diffDistance);
					}
					if (nSample==(last)) data [nSample][2]=1.0; // last measurement full weight?
					numSameDiff++;
				} else {
					for (int j=0; j<data[nSample].length;j++) data [nSample][j]=0.0;
2583

Andrey Filippov's avatar
Andrey Filippov committed
2584 2585 2586 2587 2588 2589 2590
				}
			}
			double [] coeff=(new PolynomialApproximation( debugLevel)).polynomialApproximation1d(data, 1); // just linear approximation
			int move=0;
			if ((coeff.length==1)|| (coeff[1]==0.0)){ // single sample or several for the same motor position
				move=(int) Math.round((targetMicrons-data[last][1])/this.micronsPerStep);
			} else {
2591
				move=(int) Math.round((targetMicrons-coeff[0])/coeff[1] - data [last][0]);
Andrey Filippov's avatar
Andrey Filippov committed
2592 2593 2594 2595 2596 2597 2598 2599 2600 2601 2602 2603 2604 2605 2606 2607 2608 2609 2610 2611 2612 2613 2614 2615 2616 2617 2618 2619 2620 2621 2622 2623 2624 2625 2626 2627 2628 2629 2630 2631 2632 2633
				// linear approximate
			}
			if (debugLevel>1) {
				System.out.println("focusReadjustStep(), numSameDiff="+numSameDiff+" targetMicrons="+targetMicrons+" data [last][0]="+data [last][0]);
				for (int ii=0;ii<coeff.length;ii++){
					System.out.println("focusReadjustStep(), coeff["+ii+"]="+coeff[ii]);
				}
				System.out.println("focusReadjustStep(), move="+move);
			}
			int [] result={lastPosition[0]+move,lastPosition[1]+move,lastPosition[2]+move};
			return result;
     	}

     	public int [][] probeSequence(
     			int [] motors, // should not be null
     			boolean symmetrical,
     			double probe_M1M2M3, // sigma for decay for average of the 3 motors: (M1+M2+M3)/3
     			double probe_M3_M1M2, // sigma for decay for M3 opposite to M1 and M2: M3-(M1+M2)/2
     			double probe_M2_M1){ // sigma for decay for  M2 opposite to M1:    M2-M1
     		double [] orthoCenter=motorsToOrtho(motors);
    		double [][] seqShort={
    				{probe_M1M2M3,          0,                0        },
    				{     0,           probe_M3_M1M2,         0        },
    				{     0,                0,           probe_M2_M1   },
    				{-probe_M1M2M3/3, -probe_M3_M1M2/3, -probe_M2_M1/3 },
    				{     0,                0,                0}};
    		double [][] seqSymm={
    				{ probe_M1M2M3,      0,             0},
    				{-probe_M1M2M3,      0,             0},
    				{      0,       probe_M3_M1M2,      0},
    				{      0,      -probe_M3_M1M2,      0},
    				{      0,            0,        probe_M2_M1},
    				{      0,            0,       -probe_M2_M1},
    				{      0,            0,             0}};

     		double [][] seq=symmetrical?seqSymm:seqShort;
     		int [][] motorSeq=new int [seq.length][];
     		for (int i =0; i<motorSeq.length;i++) {
     			double [] ortho= {
     					orthoCenter[0]+seq[i][0],
     					orthoCenter[1]+seq[i][1],
     					orthoCenter[2]+seq[i][2]};
2634

Andrey Filippov's avatar
Andrey Filippov committed
2635 2636 2637 2638 2639
     			motorSeq[i]= motorsFromOrtho(ortho);
     		}
     		return motorSeq;
     	}
/**
2640
 * Approximate focal distance versus sensor temperature as a linera function
Andrey Filippov's avatar
Andrey Filippov committed
2641 2642 2643 2644 2645 2646
 * @param numSamples number of last samples from history to use
 * @param lensDistanceWeightK coefficients used for focal distance calculation
 * @param lensDistanceWeightY
 * @return array of 2 elements - {length at 0C, microns/degree}
 */
     	public double [] temperatureLinearApproximation(
2647
     			boolean useLMA,
Andrey Filippov's avatar
Andrey Filippov committed
2648
     			int numSamples,             // number of last samples from history to use, 0 - use all
2649
    			double lensDistanceWeightK, // used in distance calculation
Andrey Filippov's avatar
Andrey Filippov committed
2650 2651 2652 2653 2654 2655 2656 2657 2658 2659 2660
    			double lensDistanceWeightY // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
     			){

     		if (numSamples<=0) numSamples=this.history.size();
     		if (numSamples>this.history.size()) numSamples=this.history.size();
     		double [][] data =new double [numSamples][2]; // no weights
    		double [] resolutions;
    		int firstSample=this.history.size()-numSamples;

     		for (int nSample=0;nSample<numSamples;nSample++){
				data[nSample][0]=this.history.get(firstSample+nSample).getTemperature();
2661 2662 2663 2664 2665 2666 2667
				if (useLMA){
					data[nSample][1]=this.history.get(firstSample+nSample).getzTxTy()[0];
				} else {
					resolutions= this.history.get(firstSample+nSample).getCenterResolutions();
					data[nSample][1]=getLensDistance(
							resolutions, // {R-sharpness,G-sharpness,B-sharpness}
							true, // boolean absolute, // return absolutely calibrated data
2668
							lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
2669 2670 2671 2672
							lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
							1 //debugLevel
							);
				}
Andrey Filippov's avatar
Andrey Filippov committed
2673 2674 2675 2676 2677
     		}
			PolynomialApproximation pa= new PolynomialApproximation(debugLevel);
			double [] polyCoeff=pa.polynomialApproximation1d(data, 1); // just linear
			return polyCoeff;
     	}
2678 2679

     	public double [] temperatureLinearApproximation(
2680
     			double [][] ZTT,            // Z, tXx, tY, m1,m2,m3 found with LMA - may contain NaN
2681 2682
     			int numSamples             // number of last samples from history to use, 0 - use all
     			){
2683 2684 2685
//     		if (numSamples<=0)  this.history.size();
//     		if (numSamples>this.history.size()) numSamples=this.history.size();
//    		int firstSample=this.history.size()-numSamples;
2686

2687 2688 2689
     		if (numSamples<=0) numSamples=ZTT.length;
     		if (numSamples>ZTT.length) numSamples=ZTT.length;
    		int firstSample=ZTT.length-numSamples;
2690 2691
    		int numGoodSamples=0;
     		for (int nSample=0;nSample<numSamples;nSample++){
2692
     			if (ZTT[firstSample+nSample]!=null) numGoodSamples++;
2693 2694 2695
     		}
     		double [][] data =new double [numGoodSamples][2]; // no weights
     		int index=0;
2696 2697 2698 2699
     		for (int nSample=0;nSample<numSamples;nSample++) if (ZTT[firstSample+nSample]!=null) {
//				data[index][0]=this.history.get(firstSample+nSample).getTemperature();
				data[index][0]=ZTT[firstSample+nSample][3];
				data[index++][1]=ZTT[firstSample+nSample][0]; // Z
2700 2701 2702 2703 2704
     		}
			PolynomialApproximation pa= new PolynomialApproximation(debugLevel);
			double [] polyCoeff=pa.polynomialApproximation1d(data, 1); // just linear
			return polyCoeff;
     	}
2705 2706


Andrey Filippov's avatar
Andrey Filippov committed
2707
/**
2708
 *  will return same coordinates if tolerances are met
Andrey Filippov's avatar
Andrey Filippov committed
2709 2710 2711 2712 2713 2714 2715 2716 2717 2718 2719 2720 2721 2722 2723 2724 2725 2726 2727 2728 2729 2730 2731 2732 2733 2734
 * @param useTheBest
 * @param targetMicrons
 * @param micronsFade
 * @param tiltFade
 * @param sigma_M1M2M3
 * @param sigma_M3_M1M2
 * @param sigma_M2_M1
 * @param lensDistanceWeightK
 * @param lensDistanceWeightY
 * @param weightRatioRedToGreen
 * @param weightRatioBlueToGreen
 * @param toleranceMicrons
 * @param toleranceTilt
 * @param toleranceThreshold
 * @param maxStep
 * @param debugLevel
 * @return
 */
     	public int [] focusTiltStep(
     			boolean useTheBest, // start from the best sample, (false - from the last)
     			double targetMicrons, // target focal distance
     			double micronsFade, // reduce influence of far points
     			double tiltFade,    //goodTiltSigma
     			double sigma_M1M2M3, // sigma for decay for average of the 3 motors: (M1+M2+M3)/3
     			double sigma_M3_M1M2, // sigma for decay for M3 opposite to M1 and M2: M3-(M1+M2)/2
     			double sigma_M2_M1, // sigma for decay for  M2 opposite to M1:    M2-M1
2735
    			double lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2736 2737 2738
    			double lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			double weightRatioRedToGreen,  // used for tilt averaging
    			double weightRatioBlueToGreen, // used for tilt averaging
2739

Andrey Filippov's avatar
Andrey Filippov committed
2740 2741 2742 2743 2744 2745 2746 2747 2748 2749 2750 2751 2752 2753 2754 2755 2756 2757 2758 2759 2760 2761 2762 2763 2764 2765 2766 2767 2768 2769 2770 2771 2772
    			double toleranceMicrons,
    			double toleranceTilt,
    			double toleranceThreshold,
    			double maxStep,        // Maximal motor move for each step

    			int debugLevel
     			){
     		int debugThreshold=0;
			double [] limits=absoluteCalibrationLimits();
    		int last=this.history.size()-1;
			if (limits==null){
    			String msg="Failed: lens focus quality per colors vs. distance is not yet measured/calibrated";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
			}
			if (Double.isNaN(limits[0]) || Double.isNaN(limits[1])){
    			String msg="Failed: lens focus quality per colors vs. distance is measured, but absolute calibration was not performed.";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
			}
			if ((targetMicrons<limits[0]) || (targetMicrons>limits[1])){
    			String msg="Failed: requested distance ("+targetMicrons+"um) is not withing calibrated range "+limits[0]+".."+limits[1]+"um.";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
			}
			if (this.history.size()<1){
    			String msg="Failed: history is empty.";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
2773

Andrey Filippov's avatar
Andrey Filippov committed
2774 2775 2776 2777 2778 2779
			}
    		if (
    				(targetMicrons<this.lensDistanceCalibration[0][indexMicrons]) ||
    				(targetMicrons>this.lensDistanceCalibration[this.lensDistanceCalibration.length-1][indexMicrons])) {
				System.out.println("parallelPositionMove(): out of range, this.lensDistanceCalibration[0]["+indexMicrons+"]="+
						this.lensDistanceCalibration[0][indexMicrons]+
2780
						", this.lensDistanceCalibration["+(this.lensDistanceCalibration.length-1)+"]["+indexMicrons+"]="+this.lensDistanceCalibration[this.lensDistanceCalibration.length-1][indexMicrons]);
Andrey Filippov's avatar
Andrey Filippov committed
2781 2782 2783 2784 2785 2786 2787 2788 2789 2790
    			return null ; // out of range
    		}
			if (debugLevel>1) {
				System.out.println("focusTiltStep(), last="+last);
			}
// see if tolerances are met - return with unchanged position. If not met, but close - scale down step
    		int [] lastPosition=getPosition();
    		double  worstOver= worstOverTolerance (
        			weightRatioRedToGreen,
        			weightRatioBlueToGreen,
2791
        			lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2792 2793 2794 2795 2796 2797 2798 2799 2800 2801
        			lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    				targetMicrons,
    				toleranceMicrons,
    				toleranceTilt,
        			debugLevel);
	        	if (debugLevel>debugThreshold) {
	        		System.out.println("focusTiltStep() worstOver="+IJ.d2s(worstOver,5)+", m1="+lastPosition[0]+"m2="+lastPosition[1]+"m3="+lastPosition[2]);
	        	}

    		if (worstOver<1.0) return lastPosition; // no correction if worst is less than tolerance
2802 2803 2804


    		boolean scaleDownCorrection=(worstOver<toleranceThreshold);
Andrey Filippov's avatar
Andrey Filippov committed
2805 2806
    		double scaleDownValue=scaleDownCorrection?(worstOver/toleranceThreshold):1.0;

2807

Andrey Filippov's avatar
Andrey Filippov committed
2808 2809 2810 2811 2812 2813 2814 2815 2816 2817 2818 2819 2820 2821 2822 2823 2824 2825 2826
			double [][][] data = new double[this.history.size()][3][]; //{u,v,t},{f1(u,v,t), f2(u,v,t), f3(u,v,t)),{w}
    		double kF=(micronsFade>0.0)?0.5/micronsFade/micronsFade:1.0;
    		double kT=(tiltFade>0.0)?0.5/tiltFade/tiltFade:1.0;
    		double [] kP={
    				((sigma_M1M2M3>0.0)?(0.5/sigma_M1M2M3/sigma_M1M2M3):1.0),
    				((sigma_M3_M1M2>0.0)?(0.5/sigma_M3_M1M2/sigma_M3_M1M2):1.0),
    				((sigma_M2_M1>0.0)?(0.5/sigma_M2_M1/sigma_M2_M1):1.0)};
    		double [] resolutions;
    		double bestW=0.0;
    		int bestSample=0;
			for (int nSample=0;nSample<this.history.size();nSample++){
				int [] motors=       this.history.get(nSample).motorsPos;
				data[nSample][0]=motorsToOrtho(motors );
				data[nSample][1]=new double[3];
				data[nSample][2]=new double[1];
				resolutions= this.history.get(nSample).getCenterResolutions();
				data[nSample][1][0]=getLensDistance(
						resolutions, // {R-sharpness,G-sharpness,B-sharpness}
						true, // boolean absolute, // return absolutely calibrated data
2827
						lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2828 2829 2830 2831 2832 2833 2834 2835 2836 2837 2838 2839 2840 2841 2842 2843 2844 2845 2846 2847 2848 2849
						lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
						1 //debugLevel
				)-targetMicrons;
				double [] averageMetrics= this.history.get(nSample).getMetrics(weightRatioRedToGreen,weightRatioBlueToGreen)[3]; // 4-th line - weighted average for all colors (TODO: do averaging here!!)
				data [nSample][1][1]=averageMetrics[1];
				data [nSample][1][2]=averageMetrics[2];
				data [nSample][2][0]=1.0;
				if (micronsFade>0) data[nSample][2][0]*=Math.exp(-kF*data[nSample][1][0]*data[nSample][1][0]);
				if (tiltFade>0)    data[nSample][2][0]*=Math.exp(-kT*(data[nSample][1][1]*data[nSample][1][1]+data[nSample][1][2]*data[nSample][1][2]));
				if (data[nSample][2][0]>bestW){
					bestW=data[nSample][2][0];
					bestSample=nSample;
				}
				if (debugLevel>2){
					double [] oMotors=motorsToOrtho(motors );
					System.out.println ("#"+nSample+": ["+motors[0]+","+motors[1]+","+motors[2]+"]-> "+
							IJ.d2s(oMotors[0],1)+","+IJ.d2s(oMotors[1],1)+","+IJ.d2s(oMotors[2],1)); //+" == "+
//							IJ.d2s(data[nSample][0][0],1)+","+IJ.d2s(data[nSample][0][1],1)+","+IJ.d2s(data[nSample][0][2],1));
				}

			}
			if (!useTheBest) bestSample=last;
2850
    		int [] bestPosition=getPosition(bestSample); // either last or best position
Andrey Filippov's avatar
Andrey Filippov committed
2851

2852
			// Additionally decay by motor positions
Andrey Filippov's avatar
Andrey Filippov committed
2853 2854 2855 2856 2857 2858 2859 2860 2861 2862 2863 2864 2865
			double [] orthoBest=data[bestSample][0];
			for (int nSample=0;nSample<this.history.size();nSample++){
				double [] ortho=data[nSample][0].clone();
				double w=data[nSample][2][0]; // just for debugging
				for (int i=0;i<3;i++){
					ortho[i]-=orthoBest[i];
					if (kP[i]>0) data[nSample][2][0]*=Math.exp(-kP[i]*ortho[i]*ortho[i]);
				}
				if (debugLevel>2){
					System.out.println ("##"+nSample+": ["+IJ.d2s(data[nSample][0][0],1)+","+IJ.d2s(data[nSample][0][1],1)+","+IJ.d2s(data[nSample][0][2],1)+"] "+
							" fDistDiff="+IJ.d2s(data[nSample][1][0],3)+" tiltX="+IJ.d2s(data[nSample][1][1],3)+" tiltY="+IJ.d2s(data[nSample][1][2],3)+
							" weight="+data[nSample][2][0]+" ("+w+")");
				}
2866
			}
Andrey Filippov's avatar
Andrey Filippov committed
2867 2868
			double [] orthoSolution=(new PolynomialApproximation( debugLevel)).linear3d(data);
//			double [] orthoSolution=(new PolynomialApproximation(4)).linear3d(data);
2869
			int [] newPosition=motorsFromOrtho(orthoSolution);
Andrey Filippov's avatar
Andrey Filippov committed
2870 2871 2872 2873 2874 2875 2876 2877 2878 2879 2880 2881 2882 2883 2884 2885 2886 2887 2888 2889 2890 2891 2892 2893 2894 2895 2896 2897 2898 2899 2900 2901 2902 2903 2904 2905 2906 2907 2908 2909 2910 2911 2912 2913 2914 2915 2916
			if (debugLevel>1){
				System.out.println ("Ortho solution:   "+IJ.d2s(orthoSolution[0],1)+" : "+IJ.d2s(orthoSolution[1],1)+" : "+IJ.d2s(orthoSolution[2],1));
				System.out.println ("Motors solution:  "+newPosition[0]+" : "+newPosition[1]+" : "+newPosition[2]);
			}
			if (debugLevel>0){
				System.out.println ("Best known position:"+bestPosition[0]+" : "+bestPosition[1]+" : "+bestPosition[2]);
			}
			if (toleranceThreshold<=0) return newPosition; // just testing
// Scale and limit movement if needed.
			double [] moveFromBest={
					(newPosition[0]-bestPosition[0]),
					(newPosition[1]-bestPosition[1]),
					(newPosition[2]-bestPosition[2])};
    		if (scaleDownCorrection){
    			if (debugLevel>0) System.out.println("Scaling down correction as the residual error is under scaled tolerance, multiplying by "+IJ.d2s(scaleDownValue,3));
    			for (int i=0;i<moveFromBest.length;i++) moveFromBest[i]*=scaleDownValue;
    		}
    		double corrLength=Math.sqrt(moveFromBest[0]*moveFromBest[0]+moveFromBest[1]*moveFromBest[1]+moveFromBest[2]*moveFromBest[2]);
    		if (debugLevel>debugThreshold) System.out.println("corrLength="+corrLength+" maxStep="+maxStep);
    		if (corrLength>maxStep){
    			if (debugLevel>0) System.out.println("Scaling down correction to match maximal allowed movement ("+maxStep+" steps), multiplying by "+IJ.d2s(maxStep/corrLength,3));
    			for (int i=0;i<moveFromBest.length;i++) moveFromBest[i]*=maxStep/corrLength;
    		}

    		if (debugLevel>debugThreshold) System.out.println("motorCorrections="+IJ.d2s(moveFromBest[0],0)+", "+IJ.d2s(moveFromBest[1],0)+", "+IJ.d2s(moveFromBest[2],0));
    		int [] iNewPos=new int [moveFromBest.length];
    		for (int i=0;i<moveFromBest.length;i++)  iNewPos[i]=(int) Math.round(moveFromBest[i])+bestPosition[i];
    		if (debugLevel>debugThreshold) System.out.println("new position="+iNewPos[0]+", "+iNewPos[1]+", "+iNewPos[2]);
    		return iNewPos;
     	}

     	public double [] motorsToOrtho(int [] motors ){
     		double [] ortho={
     				((motors[0]+motors[1]+motors[2])/3.0),
     				(motors[2]-0.5*(motors[0]+motors[1])),
     				(motors[1]-motors[0])};
     		return ortho;

     	}

     	public int [] motorsFromOrtho(double [] ortho ){
     		int [] motors={
     				(int) Math.round(ortho[0] -ortho[1]/3.0 -ortho[2]/2.0),
     				(int) Math.round(ortho[0] -ortho[1]/3.0 +ortho[2]/2.0),
     				(int) Math.round(ortho[0] +2*ortho[1]/3.0)};
     		return motors;
     	}
2917 2918


Andrey Filippov's avatar
Andrey Filippov committed
2919 2920 2921
     	public int [] solveFocusing(
    			double weightRatioRedToGreen,
    			double weightRatioBlueToGreen,
2922
    			double lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2923 2924 2925 2926 2927 2928 2929 2930
    			double lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
				boolean useRadialTangential,
    			double targetFarNear,
				double targetMicrons,
				double toleranceMicrons,
				double toleranceTilt,
				double toleranceThreshold,
    			double motorsSigma,
2931
    			double believeLast,   // 0 - use 'honest' best fit, 1.0 - make each plane go through the last sample
Andrey Filippov's avatar
Andrey Filippov committed
2932 2933 2934 2935 2936 2937 2938 2939 2940 2941 2942 2943 2944
    			double maxStep,        // Maximal motor move for each step
    			double goodDistanceSigma, // to use only samples with small distance errors (<=0.0 - do not use)

    			int debugLevel
    	){
     		int debugThreshold=0;
    		if (this.history.size()<4) {
    			String msg="Failed: too few samples ("+this.history.size()+") to solve";
    			System.out.println(msg);
    			IJ.showMessage("Error",msg);
    			return null;
    		}
//    		double scaleDownValue=0.5; // make a parameter? multiply correction by this if residual error is smaller than scaled tolerance
2945

Andrey Filippov's avatar
Andrey Filippov committed
2946 2947 2948 2949 2950 2951 2952 2953 2954 2955 2956 2957 2958 2959 2960 2961 2962 2963 2964 2965 2966 2967 2968 2969 2970 2971 2972 2973 2974 2975 2976 2977
    		if (!useRadialTangential){
    			// Make sure calibration is done and targetMicrons is within limits
    			double [] limits=absoluteCalibrationLimits();
    			if (limits==null){
        			String msg="Failed: lens focus quality per colors vs. distance is not yet measured/calibrated";
        			System.out.println(msg);
        			IJ.showMessage("Error",msg);
        			return null;
    			}
    			if (Double.isNaN(limits[0]) || Double.isNaN(limits[1])){
        			String msg="Failed: lens focus quality per colors vs. distance is measured, but absolute calibration was not performed.";
        			System.out.println(msg);
        			IJ.showMessage("Error",msg);
        			return null;
    			}
    			if ((targetMicrons<limits[0]) || (targetMicrons>limits[1])){
        			String msg="Failed: requested distance ("+targetMicrons+"um) is not withing calibrated range "+limits[0]+".."+limits[1]+"um.";
        			System.out.println(msg);
        			IJ.showMessage("Error",msg);
        			return null;
    			}
    		}
    		double targetValue=useRadialTangential?targetFarNear:targetMicrons;

    		int [] lastPosition=this.history.get(this.history.size()-1).motorsPos;
    		double [] lastAverageMetrics=this.history.get(this.history.size()-1).getMetrics(weightRatioRedToGreen,weightRatioBlueToGreen)[3];
    		double [] lastResolutions= this.history.get(this.history.size()-1).getCenterResolutions();

    		if (!useRadialTangential){
    			lastAverageMetrics[0]=getLensDistance(
    					lastResolutions, // {R-sharpness,G-sharpness,B-sharpness}
		    			true, // boolean absolute, // return absolutely calibrated data
2978
		    			lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2979 2980 2981 2982 2983 2984 2985
		    			lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
		    			debugLevel
		    			);
    		}
    		double  worstOver= worstOverTolerance (
        			weightRatioRedToGreen,
        			weightRatioBlueToGreen,
2986
        			lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
2987 2988 2989 2990 2991 2992 2993 2994 2995 2996 2997 2998
        			lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    				useRadialTangential,
        			targetFarNear,
    				targetMicrons,
    				toleranceMicrons,
    				toleranceTilt,
        			debugLevel);
	        	if (debugLevel>debugThreshold) {
	        		System.out.println("solveFocusing() worstOver="+IJ.d2s(worstOver,5)+", m1="+lastPosition[0]+"m2="+lastPosition[1]+"m3="+lastPosition[2]);
	        	}

    		if (worstOver<1.0) return lastPosition; // no correction if worst is less than tolerance
2999
    		boolean scaleDownCorrection=(worstOver<toleranceThreshold);
Andrey Filippov's avatar
Andrey Filippov committed
3000 3001 3002
    		double scaleDownValue=scaleDownCorrection?1.0:(worstOver/toleranceThreshold);
    		/*
    		 * Approximating each of the 3 measured parameters (Far/near, tilt x and tilt y) with linear approximation in the vicinity of the last position
3003
    		 * For each parameter F(x,y,z)=A*x + B*y +C*z + D, using Gaussian weight function with sigma= motorsSigma
Andrey Filippov's avatar
Andrey Filippov committed
3004 3005 3006 3007 3008 3009 3010 3011 3012 3013 3014 3015 3016 3017 3018 3019 3020 3021 3022 3023
    		 */
    		double kExp=0.5/motorsSigma/motorsSigma;
    		double kGoodDistance=(goodDistanceSigma>0.0)?0.5/goodDistanceSigma/goodDistanceSigma:1.0;
    		double [][] aM3=new double [3][3]; // matrix with lines of A, B, C for each parameter, first line - far/nea, second - tiltX, 3-rd tiltY
    		double [][] aB3=new double [3][1];    // right side:  {targetFarNear-D(Far/Near), -D(tiltX), -D(tiltY)}
    		for (int parNum=0;parNum<aM3.length;parNum++){
    			double S0=0.0,SX=0.0,SY=0.0,SZ=0.0,SXY=0.0,SXZ=0.0,SYZ=0.0,SX2=0.0,SY2=0.0,SZ2=0.0,SF=0.0,SFX=0.0,SFY=0.0,SFZ=0.0;
    			for (int nSample=0;nSample<this.history.size();nSample++){
    				int [] motors=       this.history.get(nSample).motorsPos;
    				double [] averageMetrics= this.history.get(nSample).getMetrics(weightRatioRedToGreen,weightRatioBlueToGreen)[3]; // 4-th line - weighted average for all colors (TODO: do averaging here!!)
    	    		double [] resolutions= this.history.get(nSample).getCenterResolutions();
    				double [] xyz={
    						motors[0]-lastPosition[0],
    						motors[1]-lastPosition[1],
    						motors[2]-lastPosition[2]};
    				double w=Math.exp(-kExp*(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]));
    	    		if ((!useRadialTangential) && (parNum>0)){ // do not apply this to focal distance, only to tilts
    	    			averageMetrics[0]=getLensDistance(
    	    					resolutions, // {R-sharpness,G-sharpness,B-sharpness}
    			    			true, // boolean absolute, // return absolutely calibrated data
3024
    			    			lensDistanceWeightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
3025 3026 3027 3028 3029 3030 3031 3032
    			    			lensDistanceWeightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			    			debugLevel
    			    			);
    	    			if (goodDistanceSigma>0.0) {
    	    				double diffDistance=(averageMetrics[0]-targetValue);
    	    				w*=Math.exp(-kGoodDistance*diffDistance*diffDistance); // reduce weight of tilt measurements if di8stance is wrong
    	    			}
    	    		}
3033
    	    		//TODO replace with PolynomialApproximation class
Andrey Filippov's avatar
Andrey Filippov committed
3034 3035 3036 3037 3038 3039 3040 3041 3042 3043 3044 3045 3046 3047 3048 3049 3050 3051 3052 3053 3054 3055 3056 3057 3058 3059 3060 3061 3062 3063 3064 3065 3066 3067 3068 3069 3070 3071 3072 3073 3074 3075 3076 3077 3078 3079 3080 3081 3082 3083 3084 3085 3086 3087 3088 3089
    				S0+=w;
    				SX+=w*xyz[0];
    				SY+=w*xyz[1];
    				SZ+=w*xyz[2];
    				SXY+=w*xyz[0]*xyz[1];
    				SXZ+=w*xyz[0]*xyz[2];
    				SYZ+=w*xyz[1]*xyz[2];
    				SX2+=w*xyz[0]*xyz[0];
    				SY2+=w*xyz[1]*xyz[1];
    				SZ2+=w*xyz[2]*xyz[2];
    				SF+=w*averageMetrics[parNum];
    				SFX+=w*averageMetrics[parNum]*xyz[0];
    				SFY+=w*averageMetrics[parNum]*xyz[1];
    				SFZ+=w*averageMetrics[parNum]*xyz[2];
    			}
    			double [][] aM={
    					{SX2,SXY,SXZ,SX},
    					{SXY,SY2,SYZ,SY},
    					{SXZ,SYZ,SZ2,SZ},
    					{SX, SY, SZ, S0}};
    			double [][] aB={
    					{SFX},
    					{SFY},
    					{SFZ},
    					{SF}};
    			Matrix M=new Matrix(aM);
    			Matrix B=new Matrix(aB);
    			// Check for singular (near-singular) here
    			double [] abcd= M.solve(B).getColumnPackedCopy();
    			if (debugLevel>debugThreshold+1){
    				System.out.println(parNum+"M:");
    				M.print(10, 5);
    				System.out.println(parNum+"B:");
    				B.print(10, 5);
    				System.out.println(parNum+"A:");
    				M.solve(B).print(10, 5);
    			}
    			//believeLast
    			aM3[parNum][0]= abcd[0];
    			aM3[parNum][1]= abcd[1];
    			aM3[parNum][2]= abcd[2];
    			aB3[parNum][0]=believeLast*lastAverageMetrics[parNum] -abcd[3]*(1.0-believeLast);

    		}
    		aB3[0][0]+=targetValue;
    		Matrix M3=new Matrix(aM3);
    		Matrix B3=new Matrix(aB3);
    		double [] motorCorrections=M3.solve(B3).getColumnPackedCopy();
    		if (debugLevel>debugThreshold) {
    			System.out.println("M3:");
    			M3.print(10, 5);
    			System.out.println("B3:");
    			B3.print(10, 5);
    			System.out.println("A3:");
    			M3.solve(B3).print(10, 5);
    		}
3090

Andrey Filippov's avatar
Andrey Filippov committed
3091
    		//Limit maximal step (proportionally?)
3092
    		// Show correction here
Andrey Filippov's avatar
Andrey Filippov committed
3093 3094 3095 3096
    		if (scaleDownCorrection){
    			if (debugLevel>0) System.out.println("Scaling down correction as the residual error is under scaled tolerance");
    			for (int i=0;i<motorCorrections.length;i++) motorCorrections[i]*=scaleDownValue;
    		}
3097

Andrey Filippov's avatar
Andrey Filippov committed
3098 3099 3100 3101 3102 3103 3104 3105 3106 3107 3108 3109
    		double corrLength=Math.sqrt(motorCorrections[0]*motorCorrections[0]+motorCorrections[1]*motorCorrections[1]+motorCorrections[2]*motorCorrections[2]);
    		if (debugLevel>debugThreshold) System.out.println("corrLength="+corrLength+" maxStep="+maxStep);
    		if (corrLength>maxStep){
    			// Show that correction is reduced
    			for (int i=0;i<motorCorrections.length;i++) motorCorrections[i]*=maxStep/corrLength;
    		}
    		if (debugLevel>debugThreshold) System.out.println("motorCorrections="+IJ.d2s(motorCorrections[0],0)+", "+IJ.d2s(motorCorrections[1],0)+", "+IJ.d2s(motorCorrections[2],0));
    		int [] iNewPos=new int [motorCorrections.length];
    		for (int i=0;i<motorCorrections.length;i++)  iNewPos[i]=(int) Math.round(motorCorrections[i])+lastPosition[i];
    		if (debugLevel>debugThreshold) System.out.println("new position="+iNewPos[0]+", "+iNewPos[1]+", "+iNewPos[2]);
    		return iNewPos;
    	}
3110
    	// todo - add "probe around" - 6/3 points, +/- for each direction (fraction of sigma?)
Andrey Filippov's avatar
Andrey Filippov committed
3111 3112

    	public void list(
3113
    			boolean useLMA,
Andrey Filippov's avatar
Andrey Filippov committed
3114 3115 3116 3117 3118 3119 3120
    			String path,
    			String lensSerial, // if null - do not add average
        		String comment,
    			boolean showIndividualComponents,
    			boolean showIndividualSamples,
    			double weightRatioRedToGreen,
    			double weightRatioBlueToGreen,
3121
    			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
3122
    			double weightY){ // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
3123 3124 3125 3126
    		list (
    				useLMA,
    				path,
    				"", // serial; number
Andrey Filippov's avatar
Andrey Filippov committed
3127 3128 3129 3130 3131 3132
    				lensSerial, // if null - do not add average
    				comment,
    				showIndividualComponents,
    				showIndividualSamples,
    				weightRatioRedToGreen,
    				weightRatioBlueToGreen,
3133
    				weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
3134 3135 3136 3137 3138 3139 3140
    				weightY,
    				false,
    				Double.NaN,
    				Double.NaN,
    				Double.NaN,
    				Double.NaN);
    	}
3141
 /*
3142 3143 3144 3145 3146 3147 3148 3149 3150 3151
    	public void addCurrentHistory(
    			FocusingField focusingField)
    	{
    		for (int i=0;i<this.history.size();i++){
    			FocusingState focusingState=this.history.get(i);
    			focusingField.addSample(
    					focusingState.getTimestamp(),
    					focusingState.getTemperature(),
    					focusingState.motorsPos,
    					focusingState.getSamples());
3152
    		}
3153 3154

    	}
3155

3156
*/
3157 3158 3159 3160 3161 3162 3163 3164 3165
    	public void saveXML(
    			String path,
        		String serialNumber,
    			String lensSerial, // if null - do not add average
        		String comment,
    			double pX0,
    			double pY0,
    			double [][][] sampleCoord){ // x,y,r
    		String sep=" "; //",";
3166
    		ArrayList<String> nodeList=new ArrayList<String>();
3167 3168 3169 3170 3171 3172 3173 3174 3175 3176 3177 3178 3179 3180 3181 3182 3183 3184 3185 3186 3187 3188 3189 3190 3191 3192 3193 3194 3195 3196 3197
        	XMLConfiguration hConfig=new XMLConfiguration();
        	hConfig.setRootElementName("focusingHistory");
        	if (comment!=null)      hConfig.addProperty("comment",comment);
        	if (serialNumber!=null) hConfig.addProperty("serialNumber",serialNumber);
        	if (lensSerial!=null)   hConfig.addProperty("lensSerial",lensSerial);
        	 hConfig.addProperty("lens_center_x",pX0);
        	 hConfig.addProperty("lens_center_y",pY0);
        	 if ((sampleCoord!=null) && (sampleCoord.length>0) && (sampleCoord[0] != null) && (sampleCoord[0].length>0)){
        		 hConfig.addProperty("samples_x",sampleCoord[0].length);
        		 hConfig.addProperty("samples_y",sampleCoord.length);
        		 for (int i=0;i<sampleCoord.length;i++)
        			 for (int j=0;j<sampleCoord[i].length;j++){
//                		 double coord[] = {sampleCoord[i][j][0],sampleCoord[i][j][1]};
                		 hConfig.addProperty("sample_"+i+"_"+j,sampleCoord[i][j][0]+sep+sampleCoord[i][j][1]);
        			 }
        	 }
        	 hConfig.addProperty("measurements",this.history.size());
        	 for (int i=0;i<this.history.size();i++){
        		 FocusingState focusingState=this.history.get(i);
        		 String prefix="measurement_"+i+".";
        		 if (focusingState.getTimestamp()!=null)
        			 hConfig.addProperty(prefix+"timestamp",focusingState.getTimestamp());
        		 hConfig.addProperty(prefix+"temperature",focusingState.getTemperature());
        		 hConfig.addProperty(prefix+"motors",focusingState.motorsPos[0]+sep+focusingState.motorsPos[1]+sep+focusingState.motorsPos[2]); // array of 3
        		 double [][][][] samples=focusingState.getSamples();
        		 nodeList.clear();
        		 if ((samples!=null) && (samples.length>0) && (samples[0]!=null) && (samples[0].length>0)){
        			 for (int ii=0;ii<samples.length;ii++) for (int jj=0;jj<samples[ii].length;jj++){
        				 if ((samples[ii][jj]!=null) && (samples[ii][jj].length>0)) {
        					 String sdata=ii+sep+jj;
        					 for (int cc=0;cc<samples[ii][jj].length;cc++) {
3198
/*
3199 3200 3201 3202
        						 double rt=samples[ii][jj][cc][0]/Math.sqrt((1+samples[ii][jj][cc][1]*samples[ii][jj][cc][1])/2.0); // tangential;
        						 double rs=rt*samples[ii][jj][cc][1]; // saggital
        						 sdata += sep+rs+ // saggital
        								  sep+rt; // tangential
3203
*/
Andrey Filippov's avatar
Andrey Filippov committed
3204 3205 3206
        						 sdata += sep+samples[ii][jj][cc][0]+ // x2
        								 sep+samples[ii][jj][cc][1]+  // y2
        								 sep+samples[ii][jj][cc][2];  // xy
3207 3208 3209 3210 3211 3212 3213 3214 3215
//        						 double [] samples_st={
//        								 samples[ii][jj][cc][0], // saggital
//        								 samples[ii][jj][cc][0]/samples[ii][jj][cc][1] // tangential (was ratio)
//        						 };
        					 }
							 nodeList.add(sdata);
        				 }
        			 }
					 hConfig.addProperty(prefix+"sample",nodeList);
3216
        		 }
3217 3218 3219 3220 3221 3222 3223 3224 3225 3226 3227 3228 3229
        	 }
         	File file=new File (path);
         	BufferedWriter writer;
 			try {
 				writer = new BufferedWriter(new FileWriter(file));
 	        	hConfig.save(writer);
 			} catch (IOException e) {
 				// TODO Auto-generated catch block
 				e.printStackTrace();
 			} catch (ConfigurationException e) {
 				// TODO Auto-generated catch block
 				e.printStackTrace();
 			}
3230

3231
    	}
Andrey Filippov's avatar
Andrey Filippov committed
3232
     	public void list(
3233
        		boolean useLMA,
Andrey Filippov's avatar
Andrey Filippov committed
3234 3235 3236 3237 3238 3239 3240 3241
    			String path,
        		String serialNumber,
    			String lensSerial, // if null - do not add average
        		String comment,
    			boolean showIndividualComponents,
    			boolean showIndividualSamples,
    			double weightRatioRedToGreen,
    			double weightRatioBlueToGreen,
3242
    			double weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
3243 3244 3245 3246 3247 3248 3249 3250 3251 3252 3253 3254 3255 3256 3257 3258 3259 3260 3261 3262 3263 3264 3265 3266 3267 3268 3269 3270 3271 3272 3273
    			double weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
    			boolean justSummary,
    			double pX0,
    			double pY0,
    			double result_lastKT,
    			double result_allHistoryKT){
    		String shortHeader="timestamp\ttemperature\ttravel\tcenter\tMotor1\tMotor2\tMotor3\tDist\tFar/Near\tTilt X\tTilt Y\tR50%\tA50%\tB50%\tC-R50%\tResol\tResol-R\tResol-G\tResol-B\tC-Resol\tC-Resol-R\tC-Resol-G\tC-Resol-B";
    		String fullHeader= "timestamp\ttemperature\ttravel\tcenter\tdir1\tMotor1\tdir2\tMotor2\tdir3\tMotor3\tDist\tComponent\tFar/Near\tTilt X\tTilt Y\tR50%\tA50%\tB50%\tC-R50%\tResol\tC-Resol\tWeight(%)";
    		String header=showIndividualComponents?fullHeader:shortHeader;
    		if (justSummary){
    			header="Serial\tLens\t"+header+"\tCenter X\tCenter Y\tLast KT\tAll KT\tComment";
    		} else {
    			header="#\t"+header;
    		}
    		int sampleRows=0;
    		int sampleColumns=0;
    		if (showIndividualSamples && !justSummary) {
    			for (int i=0;i<this.history.size();i++){
    				FocusingState focusingState=this.history.get(i);
    				double [][][][] samples=focusingState.getSamples();
    				if (samples!=null){
    					sampleRows=samples.length;
    					sampleColumns=samples[0].length;
    					break;
    				}
    			}
    		}
    		if ((sampleRows==0) || (sampleColumns==0)) showIndividualSamples=false;
    		if (showIndividualSamples){
    			if (showIndividualComponents) { // additional 3 per-color rows
    				for (int i=0;i<sampleRows;i++) for (int j=0;j<sampleColumns;j++){
Andrey Filippov's avatar
Andrey Filippov committed
3274 3275
//    					header+="\tY"+i+"X"+j+"_R"+"\tY"+i+"X"+j+"_R/T";
    					header+="\tY"+i+"X"+j+"_X2"+"\tY"+i+"X"+j+"_Y2"+"\tY"+i+"X"+j+"_XY";
Andrey Filippov's avatar
Andrey Filippov committed
3276 3277 3278 3279
    				}

    			} else {
    				for (int i=0;i<sampleRows;i++) for (int j=0;j<sampleColumns;j++){ // single row fro all colors
Andrey Filippov's avatar
Andrey Filippov committed
3280 3281 3282 3283 3284 3285
//    					header+="\tY"+i+"X"+j+"_Red_R"+"\tY"+i+"X"+j+"_Red_R/T"+
//    					"\tY"+i+"X"+j+"_Green_R"+"\tY"+i+"X"+j+"_Green_R/T"+
//    					"\tY"+i+"X"+j+"_Blue_R"+"\tY"+i+"X"+j+"_Blue_R/T";
    					header+="\tY"+i+"X"+j+"_Red_X2"+"\tY"+i+"X"+j+"_Red_Y2"+"\tY"+i+"X"+j+"_Red_XY"+
    					"\tY"+i+"X"+j+"_Green_X2"+"\tY"+i+"X"+j+"_Green_Y2"+"\tY"+i+"X"+j+"_Green_XY"+
    					"\tY"+i+"X"+j+"_Blue_X2"+"\tY"+i+"X"+j+"_Blue_Y2"+"\tY"+i+"X"+j+"_Blue_XY";
Andrey Filippov's avatar
Andrey Filippov committed
3286 3287 3288 3289 3290 3291 3292 3293
    				}
    			}
    		}
    		StringBuffer sb = new StringBuffer();
    		int [] prevPos=null;
    		if (showIndividualComponents) {
    			for (int i=0;i<this.history.size();i++){
    				FocusingState focusingState=this.history.get(i);
3294
    				double [] zTxTy=useLMA?focusingState.getzTxTy():null;
Andrey Filippov's avatar
Andrey Filippov committed
3295 3296
    				double [][] metrics=focusingState.getMetrics(weightRatioRedToGreen,weightRatioBlueToGreen);
    				double [][] resolution=focusingState.getSharpness(weightRatioRedToGreen,weightRatioBlueToGreen);
3297
    				double dist= (zTxTy==null)?getLensDistance(
Andrey Filippov's avatar
Andrey Filippov committed
3298 3299
    						focusingState.getCenterResolutions(), // {R-sharpness,G-sharpness,B-sharpness}
    		    			true, //boolean absolute, // return absolutely calibrated data
3300
    		    			weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
Andrey Filippov's avatar
Andrey Filippov committed
3301
    		    			weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
3302 3303
    		    			1): //int debugLevel
    		    				zTxTy[0];
Andrey Filippov's avatar
Andrey Filippov committed
3304
    				double []   averageMetrics=metrics[3];
3305 3306 3307 3308 3309
    				if (zTxTy!=null){
        				averageMetrics=metrics[3].clone(); // to modify w/o changing original
    					averageMetrics[1]=zTxTy[1]; // tiltX
    					averageMetrics[2]=zTxTy[2]; // tiltY
    				}
Andrey Filippov's avatar
Andrey Filippov committed
3310 3311 3312 3313 3314
    				double []   averageResolution=resolution[3];
    				sb.append((i+1)+"\t");
    				String timestamp=focusingState.getTimestamp();
    				if (timestamp==null) timestamp="n/a";
    				sb.append(timestamp+"\t");
3315

Andrey Filippov's avatar
Andrey Filippov committed
3316
    				double sensorTemperature=focusingState.getTemperature();
3317
    				sb.append(IJ.d2s(sensorTemperature,2)+"\t");
Andrey Filippov's avatar
Andrey Filippov committed
3318 3319 3320 3321 3322 3323 3324 3325 3326 3327 3328 3329 3330 3331 3332 3333 3334 3335 3336 3337 3338 3339 3340 3341 3342 3343 3344 3345 3346 3347 3348 3349 3350 3351 3352 3353 3354 3355 3356 3357 3358 3359 3360
    				if (prevPos!=null) {
    					double [] xyz={focusingState.motorsPos[0]-prevPos[0],focusingState.motorsPos[1]-prevPos[1],focusingState.motorsPos[2]-prevPos[2]};
    					sb.append(""+IJ.d2s(Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]),0));
    				}
    				sb.append(  "\t"+IJ.d2s((focusingState.motorsPos[0]+focusingState.motorsPos[1]+focusingState.motorsPos[2])/3.0,1));
    				sb.append(  "\t"+(focusingState.motorsDir[0]?"(+)":"(-)")+"\t"+focusingState.motorsPos[0]);
    				sb.append(  "\t"+(focusingState.motorsDir[1]?"(+)":"(-)")+"\t"+focusingState.motorsPos[1]);
    				sb.append(  "\t"+(focusingState.motorsDir[2]?"(+)":"(-)")+"\t"+focusingState.motorsPos[2]);
    				if (Double.isNaN(dist)) sb.append("\t---"); else sb.append(  "\t"+IJ.d2s(dist,3));
    				sb.append("\t"+this.focusingState.getName(3));
    				sb.append(  "\t"+IJ.d2s(averageMetrics[0],3));
    				sb.append(  "\t"+IJ.d2s(averageMetrics[1],3));
    				sb.append(  "\t"+IJ.d2s(averageMetrics[2],3));
    				sb.append(  "\t"+IJ.d2s(averageMetrics[3],3));
    				sb.append(  "\t"+IJ.d2s(averageMetrics[4],3));
    				sb.append(  "\t"+IJ.d2s(averageMetrics[5],3));
    				sb.append(  "\t"+IJ.d2s(averageMetrics[6],3));
    				sb.append(  "\t"+IJ.d2s(averageResolution[0],3));
    				sb.append(  "\t"+IJ.d2s(averageResolution[1],3));
    				if (showIndividualSamples){
    					for (int ii=0;ii<sampleRows;ii++) for (int jj=0;jj<sampleColumns;jj++){
    						sb.append(  "\t-\t-");
    					}
    				}
					sb.append(  "\n");
    				for (int c=0;c<3;c++){
    					sb.append(  "\t\t\t\t\t\t\t\t\t\t\t"+focusingState.getName(c));
    					sb.append(  "\t"+IJ.d2s(metrics[c][0],3));
    					sb.append(  "\t"+IJ.d2s(metrics[c][1],3));
    					sb.append(  "\t"+IJ.d2s(metrics[c][2],3));
    					sb.append(  "\t"+IJ.d2s(metrics[c][3],3));
    					sb.append(  "\t"+IJ.d2s(metrics[c][4],3));
    					sb.append(  "\t"+IJ.d2s(metrics[c][5],3));
    					sb.append(  "\t"+IJ.d2s(metrics[c][6],3));
    					sb.append(  "\t"+IJ.d2s(resolution[c][0],3));
    					sb.append(  "\t"+IJ.d2s(resolution[c][1],3));
    					sb.append(  "\t"+IJ.d2s(100*metrics[c][7],1)+"%");
        				if (showIndividualSamples){
        					double [][][][] samples=focusingState.getSamples();
        					for (int ii=0;ii<sampleRows;ii++) for (int jj=0;jj<sampleColumns;jj++){
        						if ((samples[ii][jj]!=null) && (samples[ii][jj][c]!=null)) {
        							sb.append(  "\t"+IJ.d2s(samples[ii][jj][c][0],3));
        							sb.append(  "\t"+IJ.d2s(samples[ii][jj][c][1],3));
Andrey Filippov's avatar
Andrey Filippov committed
3361
        							sb.append(  "\t"+IJ.d2s(samples[ii][jj][c][2],3));
Andrey Filippov's avatar
Andrey Filippov committed
3362
        						} else {
Andrey Filippov's avatar
Andrey Filippov committed
3363
        							sb.append(  "\t---\t---\t---");
Andrey Filippov's avatar
Andrey Filippov committed
3364 3365 3366 3367 3368 3369 3370 3371 3372 3373 3374 3375 3376 3377 3378
        						}
        					}
        				}
    					sb.append(  "\n");
    				}
    				prevPos=focusingState.motorsPos;
    			}
    		} else { // calculate and add average values (not yet for individual above)
    			int numPars=23;
    			int numDistMeas=0;
    			double [] sums= new double [numPars];
    			for (int i=0;i<sums.length;i++) sums[i]=0.0;
    			for (int i=0;i<this.history.size();i++){
//    				int parIndex=0;
    				FocusingState focusingState=this.history.get(i);
3379
    				double [] zTxTy=useLMA?focusingState.getzTxTy():null;
Andrey Filippov's avatar
Andrey Filippov committed
3380
    				double [][] metrics=focusingState.getMetrics(weightRatioRedToGreen,weightRatioBlueToGreen);
Andrey Filippov's avatar
Andrey Filippov committed
3381 3382 3383 3384 3385 3386 3387
    				double [][] resolution={{Double.NaN,Double.NaN},{Double.NaN,Double.NaN},{Double.NaN,Double.NaN},{Double.NaN,Double.NaN}};
    				try {
    					resolution=focusingState.getSharpness(weightRatioRedToGreen,weightRatioBlueToGreen);
    				} catch (Exception e){
    					System.out.println("Failed to get resolution for history("+i+")");
    					continue; // skip the whole line
    				}
Andrey Filippov's avatar
Andrey Filippov committed
3388
    				double []   averageMetrics=metrics[3];
3389 3390 3391 3392 3393
    				if (zTxTy!=null){
        				averageMetrics=metrics[3].clone(); // to modify w/o changing original
    					averageMetrics[1]=zTxTy[1]; // tiltX
    					averageMetrics[2]=zTxTy[2]; // tiltY
    				}
Andrey Filippov's avatar
Andrey Filippov committed
3394 3395 3396 3397 3398 3399
    				double []   averageResolution=resolution[3];
    				if (!justSummary) sb.append((i+1)+"\t");
    				String timestamp=focusingState.getTimestamp();
    				if (timestamp==null) timestamp="n/a";
    				else sums[0]+=Double.parseDouble(timestamp); //0
    				if (!justSummary) sb.append(timestamp+"\t");
3400
    				double sensorTemperature=focusingState.getTemperature();
Andrey Filippov's avatar
Andrey Filippov committed
3401 3402 3403 3404 3405 3406 3407 3408 3409 3410 3411 3412 3413 3414 3415 3416 3417
    				if (!justSummary) sb.append(IJ.d2s(sensorTemperature,2)+"\t");
    				sums[1]+=sensorTemperature;
    				if (prevPos!=null) {
    					double [] xyz={focusingState.motorsPos[0]-prevPos[0],focusingState.motorsPos[1]-prevPos[1],focusingState.motorsPos[2]-prevPos[2]};
    					if (!justSummary) sb.append(""+IJ.d2s(Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]),0));
    					sums[2]+=Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]+xyz[2]*xyz[2]);//2
    				}
    				if (!justSummary) sb.append(  "\t"+IJ.d2s((focusingState.motorsPos[0]+focusingState.motorsPos[1]+focusingState.motorsPos[2])/3.0,1));
    				sums[3]+=(focusingState.motorsPos[0]+focusingState.motorsPos[1]+focusingState.motorsPos[2])/3.0; // 3
    				if (!justSummary) {
    					sb.append(  "\t"+focusingState.motorsPos[0]);
    					sb.append(  "\t"+focusingState.motorsPos[1]);
    					sb.append(  "\t"+focusingState.motorsPos[2]);
    				}
					sums[4]+=focusingState.motorsPos[0];
					sums[5]+=focusingState.motorsPos[1];
					sums[6]+=focusingState.motorsPos[2];
3418 3419 3420
					double dist= (zTxTy==null)?getLensDistance(
							focusingState.getCenterResolutions(), // {R-sharpness,G-sharpness,B-sharpness}
							true, //boolean absolute, // return absolutely calibrated data
3421
							weightK, // 0.0 - all 3 component errors are combined with the same weights. 1.0 - proportional to squared first derivatives
3422 3423 3424
							weightY, // R-frac, Y-frac have the same scale regardless of the sharpness, but not Y. This is to balance Y contribution
							1): //int debugLevel
								zTxTy[0];
Andrey Filippov's avatar
Andrey Filippov committed
3425 3426 3427 3428 3429 3430 3431 3432 3433 3434 3435 3436 3437 3438 3439 3440 3441 3442 3443 3444 3445 3446 3447 3448 3449 3450 3451 3452 3453 3454
    				if (Double.isNaN(dist)){
    					if (!justSummary) sb.append("\t---");
    				} else {
    					if (!justSummary) sb.append(  "\t"+IJ.d2s(dist,3));
    					numDistMeas++;
    					sums[7]+=dist; //may be NaN
    				}
    				if (!justSummary)  {
    					sb.append(  "\t"+IJ.d2s(averageMetrics[0],3));
    					sb.append(  "\t"+IJ.d2s(averageMetrics[1],3));
    					sb.append(  "\t"+IJ.d2s(averageMetrics[2],3));
    					sb.append(  "\t"+IJ.d2s(averageMetrics[3],3));
    					sb.append(  "\t"+IJ.d2s(averageMetrics[4],3));
    					sb.append(  "\t"+IJ.d2s(averageMetrics[5],3));
    					sb.append(  "\t"+IJ.d2s(averageMetrics[6],3));
    					sb.append(  "\t"+IJ.d2s(averageResolution[0],3));
    					sb.append(  "\t"+IJ.d2s(resolution[0][0],3));
    					sb.append(  "\t"+IJ.d2s(resolution[1][0],3));
    					sb.append(  "\t"+IJ.d2s(resolution[2][0],3));
    					sb.append(  "\t"+IJ.d2s(averageResolution[1],3));
    					sb.append(  "\t"+IJ.d2s(resolution[0][1],3));
    					sb.append(  "\t"+IJ.d2s(resolution[1][1],3));
    					sb.append(  "\t"+IJ.d2s(resolution[2][1],3));
        				if (showIndividualSamples){
        					double [][][][] samples=focusingState.getSamples();
        	    			for (int ii=0;ii<sampleRows;ii++) for (int jj=0;jj<sampleColumns;jj++){
        	    				for (int cc=0;cc<3;cc++) {
        	    					if ((samples[ii][jj]!=null) && (samples[ii][jj][cc]!=null)) {
        	        					sb.append(  "\t"+IJ.d2s(samples[ii][jj][cc][0],3));
        	        					sb.append(  "\t"+IJ.d2s(samples[ii][jj][cc][1],3));
Andrey Filippov's avatar
Andrey Filippov committed
3455
        	        					sb.append(  "\t"+IJ.d2s(samples[ii][jj][cc][2],3));
Andrey Filippov's avatar
Andrey Filippov committed
3456
        	    					} else {
Andrey Filippov's avatar
Andrey Filippov committed
3457
        	    						sb.append(  "\t---\t---\t---");
Andrey Filippov's avatar
Andrey Filippov committed
3458 3459 3460 3461 3462 3463 3464 3465 3466 3467 3468 3469 3470 3471 3472 3473 3474 3475 3476 3477 3478 3479 3480 3481 3482 3483 3484 3485 3486 3487 3488 3489 3490 3491 3492 3493 3494 3495 3496 3497 3498 3499 3500 3501 3502 3503 3504 3505 3506 3507 3508 3509
        	    					}
        	    				}
        	    			}
        				}
    					sb.append(  "\n");
    				}
					sums[ 8]+=averageMetrics[0];
					sums[ 9]+=averageMetrics[1];
					sums[10]+=averageMetrics[2];
					sums[11]+=averageMetrics[3];
					sums[12]+=averageMetrics[4];
					sums[13]+=averageMetrics[5];
					sums[14]+=averageMetrics[6];
					sums[15]+=averageResolution[0];
					sums[16]+=resolution[0][0];
					sums[17]+=resolution[1][0];
					sums[18]+=resolution[2][0];
					sums[19]+=averageResolution[1];
					sums[20]+=resolution[0][1];
					sums[21]+=resolution[1][1];
					sums[22]+=resolution[2][1];
    				prevPos=focusingState.motorsPos;
    			}
    			String averageData="";
    			// how averages
    			for (int i=0;i<sums.length;i++) sums[i]/=(i==7)?numDistMeas:this.history.size();
    			averageData+=String.format("LENS%S",lensSerial)+"\t";
    			averageData+=IJ.d2s(sums[ 0],3)+"\t"; //timestamp
    			averageData+=IJ.d2s(sums[ 1],2)+"\t"; // temperature
    			averageData+=IJ.d2s(sums[ 2],0)+"\t"; // travel
    			averageData+=IJ.d2s(sums[ 3],0)+"\t"; // center
    			averageData+=IJ.d2s(sums[ 4],0)+"\t"; // motorsPos[0];
    			averageData+=IJ.d2s(sums[ 5],0)+"\t"; // motorsPos[1];
    			averageData+=IJ.d2s(sums[ 6],0)+"\t"; // motorsPos[2];
    			averageData+=((numDistMeas>0)?IJ.d2s(sums[ 7],3):"---")+"\t"; // distance;
    			averageData+=IJ.d2s(sums[ 8],3)+"\t"; // averageMetrics[0];
    			averageData+=IJ.d2s(sums[ 9],3)+"\t"; // averageMetrics[1];
    			averageData+=IJ.d2s(sums[10],3)+"\t"; // averageMetrics[2];
    			averageData+=IJ.d2s(sums[11],3)+"\t"; // averageMetrics[3];
    			averageData+=IJ.d2s(sums[12],3)+"\t"; // averageMetrics[4];
    			averageData+=IJ.d2s(sums[13],3)+"\t"; // averageMetrics[5];
    			averageData+=IJ.d2s(sums[14],3)+"\t"; // averageMetrics[6];
    			averageData+=IJ.d2s(sums[15],3)+"\t"; // averageResolution[0];
    			averageData+=IJ.d2s(sums[16],3)+"\t"; // resolution[0][0];
    			averageData+=IJ.d2s(sums[17],3)+"\t"; // resolution[1][0];
    			averageData+=IJ.d2s(sums[18],3)+"\t"; // resolution[2][0];
    			averageData+=IJ.d2s(sums[19],3)+"\t"; // averageResolution[1];
    			averageData+=IJ.d2s(sums[20],3)+"\t"; // resolution[0][1];
    			averageData+=IJ.d2s(sums[21],3)+"\t"; // resolution[1][1];
    			averageData+=IJ.d2s(sums[22],3); // resolution[2][1];
    			if (justSummary){
    				averageData=serialNumber+"\t"+averageData;
3510
        			averageData+="\t"+ IJ.d2s(pX0,1); //
Andrey Filippov's avatar
Andrey Filippov committed
3511 3512
        			averageData+="\t"+ IJ.d2s(pY0,1); //
        			averageData+="\t"+ IJ.d2s(result_lastKT,3);
3513 3514
        			averageData+="\t"+ IJ.d2s(result_allHistoryKT,3);
        			averageData+="\t"+ comment; //
Andrey Filippov's avatar
Andrey Filippov committed
3515 3516 3517 3518 3519 3520 3521 3522 3523 3524 3525 3526 3527 3528 3529 3530 3531 3532 3533 3534 3535 3536 3537 3538 3539 3540 3541 3542 3543 3544 3545 3546
    			}
    			averageData+="\n";
    			// TODO: add data for individula samples?
    			sb.append(averageData);
    			if (path!=null){
    				System.out.println("Averaging measurements:");
    				System.out.println(header);
    				System.out.println(averageData);
    			}
    		}
    		if (path!=null) {
    			String footer= justSummary?"":("\nComment: "+comment+"\tdX=\t"+ IJ.d2s(pX0,1)+"\tdY=\t"+IJ.d2s(pY0,1));
    			CalibrationFileManagement.saveStringToFile (
						path,
						header+"\n"+sb.toString()+footer);
    		}else{
        		new TextWindow("History", header, sb.toString(), showIndividualComponents?1500:1700,500);
    		}
    	}

    	public void setLastProbed(){
    		if (this.history.size()==0){
    			String msg="List is empty";
    			IJ.showMessage("Error",msg);
    			throw new IllegalArgumentException (msg);
    		}
    		this.history.get(this.history.size()-1).isProbed=true;

    	}

    	public double distFromProbed(){
    		if (this.history.size()>0) return distFromProbed (null);
3547
    		else return Double.NaN;
Andrey Filippov's avatar
Andrey Filippov committed
3548 3549 3550
    	}
    	public double distFromProbed(int [] position){
    		double minDist=Double.NaN;
3551
    		if (this.history.size() <1) return Double.NaN;
Andrey Filippov's avatar
Andrey Filippov committed
3552
    		if (position==null) position= this.history.get(this.history.size()-1).motorsPos;
3553

Andrey Filippov's avatar
Andrey Filippov committed
3554 3555 3556 3557 3558 3559 3560 3561 3562 3563 3564 3565
    		if (this.history.size()>0){
    			for (int i=0;i<this.history.size();i++) if (this.history.get(i).isProbed) {
    				int [] probedPosition=this.history.get(i).motorsPos;
    				double dist=Math.sqrt(
    						(position[0]-probedPosition[0])*(position[0]-probedPosition[0])+
    						(position[1]-probedPosition[1])*(position[1]-probedPosition[1])+
    						(position[2]-probedPosition[2])*(position[2]-probedPosition[2]));
    				if (!(minDist<dist)){ //minDist may be NaN
    					minDist=dist;
    				}
    			}
    		}
3566
    		return minDist;
Andrey Filippov's avatar
Andrey Filippov committed
3567 3568
    	}

3569

Andrey Filippov's avatar
Andrey Filippov committed
3570 3571 3572
    	public void clear(){
    		this.history.clear();
    		this.currentSigma=Double.NaN;
3573
    		this.calibrationHistoryFromTo=null; // history indices, to-from (inclusive) used to calculate calibration
Andrey Filippov's avatar
Andrey Filippov committed
3574 3575 3576 3577 3578 3579 3580 3581 3582 3583 3584 3585 3586 3587 3588 3589 3590 3591 3592 3593 3594 3595 3596 3597
    		this.hysteresisHistoryFromTo=null; // history indices, to-from (inclusive) used to calculate hysteresis
    	}
    	public void add(
    			String sTimestamp,
    			double temperature,
    			int [] motors, // FocusingMotors.curpos
    			boolean [] lastDirs, // FocusingMotors.lastDirs
    			double[][] psfMetrics, // null OK
    			double [][][][] fullResults // null OK
    	){
    		this.history.add(new FocusingState(sTimestamp,temperature,motors,lastDirs,psfMetrics,fullResults));
    	}
    	/*
    	      	public void add(
    			String sTimestamp,
    			double temperature,
    			int [] motors, // FocusingMotors.curpos
    			boolean [] lastDirs, // FocusingMotors.lastDirs
    			double[][] psfMetrics // null OK
    	){
    		this.history.add(new FocusingState(sTimestamp,temperature,motors,lastDirs,psfMetrics));

    	}
        */
3598 3599


Andrey Filippov's avatar
Andrey Filippov committed
3600 3601 3602 3603 3604 3605 3606 3607 3608 3609 3610 3611 3612 3613 3614 3615 3616 3617 3618 3619 3620 3621 3622 3623 3624 3625 3626 3627 3628 3629 3630 3631 3632 3633 3634 3635 3636 3637 3638 3639 3640 3641 3642 3643 3644 3645 3646 3647
    	public void setLastMetrics(
    			double[][] psfMetrics
    	){
    		if (this.history.size()==0){
    			String msg="List is empty";
    			IJ.showMessage("Error",msg);
    			throw new IllegalArgumentException (msg);
    		}
    		this.history.get(this.history.size()-1).setMetrics(psfMetrics);
    	}


    	public class FocusingState{
    		public double temperature;
    		public int [] motorsPos=new int[3];
    		public boolean [] motorsDir=new boolean[3]; // true if was increasing
    		private double [][] psfMetricses=null;
    		boolean isProbed=false;
    		private int [] indices={1,5,2,6};
    		private String [] names={"red","green","blue","AVERAGE"};
    		public int indexFarNear=0;
    		public int indexTiltX=  1;
    		public int indexTiltY=  2;
    		public int indexR50All= 3;
    		public int indexA50All= 4;
    		public int indexB50All= 5;
    		public int indexR50Center= 6;
//    		public int indexA50Center= 7;
//    		public int indexB50Center= 8;
    		public String sTimestamp=null;
    		public double [][][][] fullResults=null;

    		public FocusingState(){}
    		public FocusingState(
    				String sTimestamp,
    				double temperature,
    				int [] motorsPos,
    				boolean [] motorsDir,
    				double [][] psfMetrics,
    				double [][][][] fullResults){
    			this.temperature=temperature;
    			this.motorsPos=motorsPos.clone();
    			this.motorsDir=motorsDir.clone();
    			this.fullResults=null;
    			if (fullResults!=null){ // deep clone
    				this.fullResults=new double [fullResults.length][][][];
    				for (int i=0;i<fullResults.length;i++){
    					if (fullResults[i]!=null){
3648
    						this.fullResults[i]=new double [fullResults[i].length][][];
Andrey Filippov's avatar
Andrey Filippov committed
3649 3650 3651 3652 3653 3654 3655 3656 3657 3658 3659 3660 3661
    						for (int j=0;j<fullResults[i].length;j++){
    	    					if (fullResults[i][j]!=null){
    	    						this.fullResults[i][j]=new double [fullResults[i][j].length][];
    	    						for (int c=0;c<fullResults[i][j].length;c++){
    	    							if (fullResults[i][j][c]!=null){
    	    								this.fullResults[i][j][c]=new double [fullResults[i][j][c].length]; // 2
    	    								for (int n=0;n<fullResults[i][j][c].length;n++){
    	    									this.fullResults[i][j][c][n]=fullResults[i][j][c][n];
    	    								}
    	    							} else this.fullResults[i][j][c]=null;
    	    						}
    	    					} else this.fullResults[i][j]=null;
    						}
3662
    					} else this.fullResults[i]=null;
Andrey Filippov's avatar
Andrey Filippov committed
3663 3664 3665 3666 3667 3668 3669
    				}
    			}

    			if (psfMetrics!=null) setMetrics(psfMetrics);
    			if (sTimestamp!=null) this.sTimestamp=sTimestamp;

    		}
3670
/*
Andrey Filippov's avatar
Andrey Filippov committed
3671 3672 3673 3674 3675 3676 3677 3678 3679 3680 3681 3682 3683 3684 3685 3686 3687 3688 3689 3690 3691 3692 3693 3694 3695 3696 3697 3698 3699 3700 3701 3702 3703 3704 3705 3706 3707 3708 3709 3710 3711 3712 3713 3714 3715 3716 3717 3718 3719 3720 3721 3722 3723 3724 3725 3726 3727 3728 3729 3730 3731 3732 3733 3734 3735 3736 3737 3738 3739 3740 3741 3742 3743 3744 3745 3746 3747 3748 3749 3750 3751 3752 3753 3754 3755 3756 3757 3758 3759 3760 3761 3762 3763 3764 3765 3766 3767 3768 3769 3770 3771 3772 3773 3774 3775 3776 3777 3778 3779 3780 3781 3782 3783 3784 3785 3786
    		public FocusingState(
    				String sTimestamp,
    				double temperature,
    				int [] motorsPos,
    				boolean [] motorsDir,
    				double [][] psfMetrics){
    			this.temperature=temperature;
    			this.motorsPos=motorsPos.clone();
    			this.motorsDir=motorsDir.clone();
    			this.fullResults=null;
    			if (psfMetrics!=null) setMetrics(psfMetrics);
    			if (sTimestamp!=null) this.sTimestamp=sTimestamp;

    		}
*/
    		public FocusingState(
    				String sTimestamp,
    				double temperature,
    				int [] motorsPos,
    				boolean [] motorsDir){
    			this.temperature=temperature;
    			this.motorsPos=motorsPos.clone();
    			this.motorsDir=motorsDir.clone();
    			if (sTimestamp!=null) this.sTimestamp=sTimestamp;
    		}

    		public FocusingState(
    				String sTimestamp,
    				double temperature,
    				int [] motorsPos,
    				FocusingState previousState){
    			this.temperature=temperature;
    			this.motorsPos=motorsPos.clone();
    			setDirsFromPrevious(previousState);
    			if (sTimestamp!=null) this.sTimestamp=sTimestamp;
    		}

    		public FocusingState(
    				String sTimestamp,
    				double temperature,
    				int [] motorsPos,
    				FocusingState previousState,
    				double [][] psfMetrics){
    			this.temperature=temperature;
    			this.motorsPos=motorsPos.clone();
    			setDirsFromPrevious(previousState);
    			if (psfMetrics!=null) setMetrics(psfMetrics);
    			if (sTimestamp!=null) this.sTimestamp=sTimestamp;
    		}

    		public FocusingState(
    				String sTimestamp,
    				double temperature,
    				int m1,
    				int m2,
    				int m3,
    				FocusingState previousState,
    				double [][] psfMetrics){
    			this.temperature=temperature;
    			int [] motorsPos={m1,m2,m3};
    			this.motorsPos=motorsPos.clone();
    			setDirsFromPrevious(previousState);
    			if (psfMetrics!=null) setMetrics(psfMetrics);
    			if (sTimestamp!=null) this.sTimestamp=sTimestamp;

    		}
    		public void setMetrics(double [][] psfMetrics){
    			if (psfMetrics!=null) {
    				this.psfMetricses=new double [psfMetrics.length][];
    				for (int i=0;i<psfMetrics.length;i++){
    					if (psfMetrics[i]!=null){
    						this.psfMetricses[i]=psfMetrics[i].clone();
    					} else this.psfMetricses[i]=null;
    				}
    			} else this.psfMetricses=null;

    		}
    		public void setDirsFromPrevious(
    				FocusingState previousState){
    			for (int i=0;i<motorsPos.length;i++){
    				if (previousState.motorsPos[i]==this.motorsPos[i]) this.motorsDir[i]=previousState.motorsDir[i];
    				else if (previousState.motorsPos[i]<this.motorsPos[i]) this.motorsDir[i]=true;
    				else this.motorsDir[i]=false;
    			}
    		}
    		public double getTemperature(){
    			return this.temperature;
    		}
    		public String getTimestamp(){
    			return this.sTimestamp;
    		}

    		public String getName(int i){
    			return this.names[i];
    		}

    		public int [] getMotors(){
    			return this.motorsPos.clone();
    		}
    		public double getMotorsAverage(){
    			return (this.motorsPos[0]+this.motorsPos[1]+this.motorsPos[2])/3.0;
    		}
    		public double [][] getSharpness(
    				double weightRatioRedToGreen,
    				double weightRatioBlueToGreen){
    			return getSharpness(
    					weightRatioRedToGreen,
    					1.0,
    					weightRatioBlueToGreen);
    		}
    		public double [][] getSharpness( // first index - color+average, second index 0 - all, 1 - center
    				double weightRed,
    				double weightGreen,
    				double weightBlue){
    			double [][] metrics=getMetrics(weightRed, weightGreen, weightBlue);
    			double [][]sharpness={
Andrey Filippov's avatar
Andrey Filippov committed
3787
    					{1.0/metrics[0][this.indexR50All],1.0/metrics[0][this.indexR50Center]}, //java.lang.NullPointerException
Andrey Filippov's avatar
Andrey Filippov committed
3788 3789 3790 3791 3792 3793 3794 3795 3796 3797 3798
    					{1.0/metrics[1][this.indexR50All],1.0/metrics[1][this.indexR50Center]},
    					{1.0/metrics[2][this.indexR50All],1.0/metrics[2][this.indexR50Center]},
    					{0.0,0.0}};
    			double [] weights={weightRed,weightGreen,weightBlue};
    			double w=0.0;
    			for (int i=0;i<weights.length;i++) w+=weights[i];
    			for (int i=0;i<weights.length;i++) weights[i]/=w;
    			for (int i=0;i<3;i++)for (int j=0;j<2;j++) sharpness[3][j]+=weights[i]*sharpness[i][j];
    			return sharpness;
    		}
    		public double [] getCenterResolutions(){
3799
    			try {
Andrey Filippov's avatar
Andrey Filippov committed
3800 3801 3802 3803 3804
    			double [] resolutions={
    					1.0/this.psfMetricses[this.indices[0]][this.indexR50Center],  // red, r50% center
    					1.0/this.psfMetricses[this.indices[1]][this.indexR50Center],  // green, r50% center
    					1.0/this.psfMetricses[this.indices[2]][this.indexR50Center]}; // blue, r50% center
    			return resolutions;
3805 3806 3807
    			} catch (Exception e){
    				return null;
    			}
Andrey Filippov's avatar
Andrey Filippov committed
3808
    		}
3809

Andrey Filippov's avatar
Andrey Filippov committed
3810 3811 3812 3813 3814 3815 3816 3817
    		public double [][] getMetrics(
    				double weightRatioRedToGreen,
    				double weightRatioBlueToGreen){
    			return getMetrics(
    					weightRatioRedToGreen,
    					1.0,
    					weightRatioBlueToGreen);
    		}
3818 3819 3820 3821 3822 3823 3824 3825 3826 3827 3828
    		// alternative mode (from aberration model) ising metrics[6][]
    		public double [] getzTxTy(){
    			if (this.psfMetricses==null){
    				return null;
    			}
    			if ((this.psfMetricses.length<7) || (this.psfMetricses[6]==null)) {
    				System.out.println("BUG? psfMetrics does not have line 6 with lens berration model z, tx, ty");
    				return null;
    			}
    			double [] zTxTy={this.psfMetricses[6][0],this.psfMetricses[6][1],this.psfMetricses[6][2]};
    			return zTxTy;
Andrey Filippov's avatar
Andrey Filippov committed
3829

3830
    		}
Andrey Filippov's avatar
Andrey Filippov committed
3831 3832 3833 3834 3835 3836 3837 3838 3839 3840
    		public double [][] getMetrics(
    				double weightRed,
    				double weightGreen,
    				double weightBlue){
    			boolean [] squared={false,false,false,false,true,true,false,false};
    			double [][] metrics=new double[4][];
    			double [] weights={weightRed,weightGreen,weightBlue};
    			double w=0.0;
    			for (int i=0;i<weights.length;i++) w+=weights[i];
    			for (int i=0;i<weights.length;i++) weights[i]/=w;
3841 3842 3843 3844 3845 3846
    			try {
    				metrics[3]=new double[this.psfMetricses[this.indices[0]].length];
    				for (int c=0;c<3;c++) metrics[c]= this.psfMetricses[this.indices[c]];
    			} catch (Exception e){
    				return null;
    			}
Andrey Filippov's avatar
Andrey Filippov committed
3847 3848 3849 3850 3851 3852 3853 3854 3855 3856 3857 3858 3859 3860 3861 3862 3863
    			for (int i=0; i<metrics[3].length;i++) {
    				metrics[3][i]=0.0;
    				if (squared[i]){
    					for (int c=0;c<3;c++) metrics[3][i]+=metrics[c][i]*metrics[c][i]*weights[c];
    					metrics[3][i]=Math.sqrt(metrics[3][i]);
    				} else for (int c=0;c<3;c++) metrics[3][i]+=metrics[c][i]*weights[c];
    			}
    			return metrics;
    		}
    		public double [][][][] getSamples(){
    			return this.fullResults;
    		}

    	}
    }


3864
    }