Implemented setInitialOrientationsCuas

src/main/java/com/elphel/imagej/cuas/Cuas.java

 package com.elphel.imagej.cuas;
 
+import java.awt.Rectangle;
+import java.io.File;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.concurrent.atomic.AtomicInteger;
+
+import com.elphel.imagej.cameras.CLTParameters;
+import com.elphel.imagej.common.PolynomialApproximation;
+import com.elphel.imagej.common.ShowDoubleFloatArrays;
+import com.elphel.imagej.gpu.GPUTileProcessor;
+import com.elphel.imagej.tileprocessor.ErsCorrection;
+import com.elphel.imagej.tileprocessor.ImageDtt;
+import com.elphel.imagej.tileprocessor.OpticalFlow;
+import com.elphel.imagej.tileprocessor.QuadCLT;
+import com.elphel.imagej.tileprocessor.TileProcessor;
+
 public class Cuas {
+    public static QuadCLT createCenterClt( // assuming cuas_rotation is true
+    		CLTParameters  clt_parameters,
+    		QuadCLT []     quadCLTs,
+    		int            ref_index, // where combo_dsi is
+    		int     []     range, // or null
+    		double  [][]   ref_combo_dsi, // DSI data for the reference scene (or null to read it from file)
+    		boolean        condition_dsi,
+    		int            sensor_mask, // -1 - all;
+    		int            debugLevel) {
+    	QuadCLT ref_clt = quadCLTs[ref_index];
+    	QuadCLT last_clt = quadCLTs[quadCLTs.length-1]; // save center with the latest timestamp
+    	if (range == null) {
+    		range = new int [] {0, quadCLTs.length-1};
+    	}
+    	if (ref_index <0) {
+    		ref_index = quadCLTs.length-1;
+    	}
+		double [][] center_ATR = CuasCenterLma.getCenterATR(
+				quadCLTs,    // QuadCLT [] quadCLTs,
+				ref_index,   //int        ref_index,
+				range,  //int     []     range,
+				debugLevel); // int debugLevel);
+		double [] cuas_atr = new double [] { center_ATR[0][0], center_ATR[0][1], center_ATR[0][2]};
+		if (ref_combo_dsi == null) {
+			ref_combo_dsi =quadCLTs[ref_index].restoreComboDSI(true); // also sets quadCLTs[ref_index].dsi and blue sky
+		}
+    	double [][] dls = { 
+    			ref_combo_dsi[OpticalFlow.COMBO_DSN_INDX_DISP],
+    			ref_combo_dsi[OpticalFlow.COMBO_DSN_INDX_LMA],
+    			ref_combo_dsi[OpticalFlow.COMBO_DSN_INDX_STRENGTH]
+    	};
+    	double [][] ds = new double [][] {dls[0],dls[2]}; // {disparity, strength}
+		if (condition_dsi) {
+			ds = OpticalFlow.conditionInitialDS( // only
+        			true,                // boolean        use_conf,       // use configuration parameters, false - use following  
+        			clt_parameters,      // CLTParameters  clt_parameters,
+        			dls,                 // double [][]    dls
+        			quadCLTs[ref_index], // QuadCLT        scene,
+        			debugLevel);			
+		}
+		// getRefPxPyD: pX, pY reference image X,Y correspomnding to the uniform grid of the center view
+        final double [][] ref_pXpYD = getRefPxPyD( // Use to interpolate disparity layers
+        		clt_parameters, // CLTParameters  clt_parameters,
+        		true,           // boolean        mode_cuas,            		
+        		null,           // Rectangle      fov_tiles,
+        		OpticalFlow.ZERO3,          // double []      stereo_xyz, // offset reference camera {x,y,z}
+        		cuas_atr,       // double []      stereo_atr_in, // offset reference orientation (cuas)
+        		ds[0],          // double []      ref_disparity,			
+        		ref_clt,        // QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+        		debugLevel);    // int            debugLevel)
+    	boolean save_weights = true; // always
+		float [][] center_clt_w = getTDComboSceneSequence(
+				clt_parameters,      // CLTParameters  clt_parameters,
+        		ref_pXpYD,           // double [][]    ref_pXpYD,    		
+				save_weights,        // boolean        save_weights, // output corresponding weights for each data
+				true,                // boolean        merge_all,
+				sensor_mask,         // int            sensor_mask,
+	    		null,                // Rectangle      fov_tiles,
+	    		OpticalFlow.ZERO3,               // double []      stereo_xyz, // offset reference camera {x,y,z}
+	    		cuas_atr,            // double []      stereo_atr_in, // offset reference orientation (cuas)
+	    		null, // ds[0],      // double []      ref_disparity,	// may be null if ref_pXpYD != null	
+	    		quadCLTs,            // QuadCLT []     quadCLTs,
+	    		ref_clt,             // QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+	    		debugLevel);         // int            debugLevel)
+		// remove weights
+		float [][] center_clt = save_weights? new float [center_clt_w.length/2][]: center_clt_w;
+		int [][] center_clt_num = save_weights? new int [center_clt_w.length/2][]: null;
+		if (save_weights) { // remove second half (weights as integer numbers)
+			for (int i = 0; i< center_clt.length; i++) {
+				center_clt[i] = center_clt_w[i];
+				center_clt_num[i] = new int [center_clt_w[i + center_clt.length].length];
+				for (int j = 0; j < center_clt_num[i].length; j++) {
+					center_clt_num[i][j] = (int) Math.round(center_clt_w[i + center_clt.length][j]);
+				}
+			}
+		}
+		int [] wh_c = ref_clt.getWHC(false);
+		//final int tile_size_td = 4 * GPUTileProcessor.DTT_SIZE * GPUTileProcessor.DTT_SIZE;
+		final int width_clt = wh_c[0] * 2; // to make image dimensions similar
+		
+		String center_name =  QuadCLT.getCenterDirName(last_clt.getImageName()); // make name from last timestamp, not reference
+		String ref_dir_path = ref_clt.getX3dDirectory(center_name);
+		File cdir = new File(ref_dir_path);
+		QuadCLT center_CLT = new QuadCLT(ref_clt, center_name); //null
+		cdir.mkdirs();
+		center_CLT.setImagePath(cdir.getPath());
+		center_CLT.setCenterClt(
+				center_clt[0], // float [] clt,
+				center_clt_num[0], // int []   clt_num,
+				width_clt); // int      clt_width);
+		TileProcessor tp =         quadCLTs[quadCLTs.length - 1].getTileProcessor();
+		final int transform_size = tp.getTileSize();
+		final int tilesX =         tp.getTilesX();
+		final int tilesY =         tp.getTilesY();
+		final double [][] center_combo_dsi = interpolateDSI(
+				ref_pXpYD,      // final double [][] ref_pXpYD,
+				ref_combo_dsi, // combo_dsi_cli,  // final double [][] dsi_in,
+				transform_size, // final int         transform_size,
+				tilesX,         // final int         tilesX,
+				tilesY);        //final int         tilesY);
+		center_CLT.setDSIFromCombo(center_combo_dsi); // reformat
+//		center_CLT.setDSI(center_combo_dsi); // WRONG!
+		String rslt_suffix = "-INTER-INTRA";
+		rslt_suffix += (clt_parameters.correlate_lma?"-LMA":"-NOLMA");
+		/*
+		// fixing NaN in strengths. It is uses to return RMS in Not needed - NaN was from Arrays.fill(combo_dsn_final[i], Double.NaN);
+		// OpticalFlow:10348
+		for (int slice: OpticalFlow.COMBO_DSN_NONNAN) { // new int[] {COMBO_DSN_INDX_STRENGTH,COMBO_DSN_INDX_STRENGTH_BG}) {
+			if (center_combo_dsi[slice] != null) {
+				for (int i = 0; i <center_combo_dsi[slice].length; i++) {
+					if (Double.isNaN(center_combo_dsi[slice][i])) {
+						center_combo_dsi[slice][i] = 0.0;
+					}
+				}
+			}
+		}
+		*/
+
+		center_CLT.saveDoubleArrayInModelDirectory( // error
+				rslt_suffix,                  // String      suffix,
+				OpticalFlow.COMBO_DSN_TITLES, // combo_dsn_titles_full, // null,          // String []   labels, // or null
+				center_combo_dsi,             // dbg_data,         // double [][] data,
+				tilesX,                       // int         width,
+				tilesY);                      // int         height)
+//		center_CLT.saveDSI();
+		center_CLT.saveCenterClt();
+
+		
+		/*
+		boolean save_clt = true;
+		boolean save_dsi = true;
+		boolean save_in_ref = true;
+		boolean save_in_last = true;
+		if (save_clt) {
+			final int height_clt = combo_seq_clt[0].length/width_clt;
+			String [] clt_titles = new String [combo_seq_clt_w.length];
+			for (int i = 0; i < combo_seq_clt.length; i++) {
+				clt_titles[i] = "chn-"+i;
+				if (save_weights) {
+					clt_titles[i+combo_seq_clt.length] = "weight-"+i;
+				}
+			}
+			String suffix_clt = "-clt_combo";
+			if (save_weights) {
+				suffix_clt+="_weights";
+			}
+			if (save_in_ref) {
+			ref_clt.saveFloatArrayInModelDirectory( // error
+					suffix_clt,     // String      suffix,
+					clt_titles,       // combo_dsn_titles_full, // null,          // String []   labels, // or null
+					combo_seq_clt_w,    // dbg_data,         // float [][] data,
+					width_clt,        // int         width,
+					height_clt);      // int         height)
+			}
+			if (save_in_last) {
+				quadCLTs[quadCLTs.length - 1].saveFloatArrayInModelDirectory( // error
+						suffix_clt,     // String      suffix,
+						clt_titles,       // combo_dsn_titles_full, // null,          // String []   labels, // or null
+						combo_seq_clt_w,    // dbg_data,         // float [][] data,
+						width_clt,        // int         width,
+						height_clt);      // int         height)
+			}
+		}
+		
+		if (save_dsi) {
+//			final double [][] combo_dsi_cli = ref_clt.restoreComboDSI (true);
+			TileProcessor tp =         quadCLTs[quadCLTs.length - 1].getTileProcessor();
+			final int transform_size = tp.getTileSize();
+			final int tilesX =         tp.getTilesX();
+			final int tilesY =         tp.getTilesY();
+			final double [][] center_dsi = interpolateDSI(
+					ref_pXpYD,      // final double [][] ref_pXpYD,
+					ref_dsi, // combo_dsi_cli,  // final double [][] dsi_in,
+					transform_size, // final int         transform_size,
+					tilesX,         // final int         tilesX,
+					tilesY);        //final int         tilesY);
+			
+			String rslt_suffix = "-CENTER-INTER-INTRA";
+			// use quadCLTs[quadCLTs.length - 1].restoreComboDSI ("-CENTER",true); to restore
+			
+			rslt_suffix += (clt_parameters.correlate_lma?"-LMA":"-NOLMA");
+	    	if (save_in_ref) {
+	    		ref_clt.saveDoubleArrayInModelDirectory( // error
+	    				rslt_suffix,           // String      suffix,
+	    				OpticalFlow.COMBO_DSN_TITLES, // combo_dsn_titles_full, // null,          // String []   labels, // or null
+	    				center_dsi,       // dbg_data,         // double [][] data,
+	    				tilesX,                // int         width,
+	    				tilesY);               // int         height)
+	    	}
+	    	if (save_in_last) {
+	    		quadCLTs[quadCLTs.length - 1].saveDoubleArrayInModelDirectory( // error
+	    				rslt_suffix,           // String      suffix,
+	    				OpticalFlow.COMBO_DSN_TITLES, // combo_dsn_titles_full, // null,          // String []   labels, // or null
+	    				center_dsi,       // dbg_data,         // double [][] data,
+	    				tilesX,                // int         width,
+	    				tilesY);               // int         height)
+	    	}
+		}
+		*/
+		return center_CLT;
+    }
+	public static double [][] interpolateDSI( // uses combo_dsi, not this.dsi
+			final double [][] ref_pXpYD,
+			final double [][] combo_dsi_in,
+			final int         transform_size,
+			final int         tilesX,
+			final int         tilesY){
+		final double [][] combo_dsi = new double [combo_dsi_in.length][];
+		final int tiles = tilesY * tilesX;
+		//COMBO_DSN_NONNAN
+		boolean [] non_nan = new boolean[combo_dsi_in.length];
+		for (int indx:OpticalFlow.COMBO_DSN_NONNAN) {
+			if (indx < non_nan.length) {
+				non_nan[indx] = true;
+			}
+		}
+		for (int slice = 0; slice < combo_dsi_in.length; slice++) if (combo_dsi_in[slice] != null) {
+			combo_dsi[slice] = new double [tiles];
+			if (!non_nan[slice]) {
+				Arrays.fill(combo_dsi[slice], Double.NaN);// improve
+			}
+		}
+		final Thread[] threads = ImageDtt.newThreadArray(ImageDtt.THREADS_MAX);
+		final AtomicInteger ai = new AtomicInteger(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+//					PolynomialApproximation pa = new PolynomialApproximation();		
+					for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) if (ref_pXpYD[nTile] != null){
+						double [] pXpYD = ref_pXpYD[nTile];
+						double px = pXpYD[0];
+						double py = pXpYD[1];
+						double tx = px/transform_size-0.5;
+						double ty = py/transform_size-0.5;
+						int tx0 = (int) Math.floor(tx);
+						int ty0 = (int) Math.floor(ty);
+						int tx1 = tx0+1;
+						int ty1 = ty0+1;
+						double kx = tx-tx0;
+						double ky = ty-ty0;
+						if (tx0 < 0) {
+							tx0 = 0;
+							tx1 = tx0+1;
+							kx = 0;
+						} else if (tx1 >= tilesX) {
+							tx1 = tilesX - 1;
+							tx0 = tx1-1;
+							kx = 1.0;
+						}
+						if (ty0 < 0) {
+							ty0 = 0;
+							ty1 = ty0+1;
+							ky = 0;
+						} else if (ty1 >= tilesY) {
+							ty1 = tilesY - 1;
+							ty0 = ty1-1;
+							ky = 1.0;
+						}
+						for (int slice = 0; slice < combo_dsi_in.length; slice++) 	if (combo_dsi_in[slice] != null) {
+							double [][] d4 = {
+									{combo_dsi_in[slice][tilesX*ty0 + tx0],combo_dsi_in[slice][tilesX*ty0 + tx1]},
+									{combo_dsi_in[slice][tilesX*ty1 + tx0],combo_dsi_in[slice][tilesX*ty1 + tx1]},
+							};
+							int dn = 0;
+							for (int i = 0; i < 2; i++) { 
+								for (int j = 0; j < 2; j++) {
+									if (Double.isNaN(d4[i][j])) {
+										dn |= 1 << (2*i + j);
+									}
+								}
+							}
+							if (dn == 15) continue; // all NaNs
+							if (dn != 0) { // fill NaNs
+								for (int i = 0; i < 2; i++) { 
+									for (int j = 0; j < 2; j++) {
+										if ((dn & (1 << (2*i + j))) != 0) { // is NaN
+											int hnb =  (1 << (2*i + 1 -j)); // horizontal neighbor bit
+											int vnb =  (1 << (2 - 2*i +j)); // horizontal neighbor
+											if ((dn & hnb) == 0) { // non-NaN
+												if ((dn & vnb) == 0) { // non-NaN
+													// use average of 2 neighbors, ignore diagonal
+													d4[i][j] = 0.5*(d4[i][1-j]+d4[1-i][j]);
+												} else { // if ((dn & vnb) == 0) {
+													d4[i][j] = d4[i][1-j]; // use horizontal neighbor
+												}
+											} else { // if ((dn & hnb) == 0)
+												if ((dn & vnb) == 0) {  // non-NaN
+													d4[i][j] = d4[1-i][j]; // use vertical neighbor
+												} else { // if ((dn & vnb) == 0) {
+													int dnb =  (1 << (2 - 2*i + 1 - j)); // diagonal neighbor
+													d4[i][j] = d4[1-i][1-j]; // copy diagonal
+												}
+											}
+										}
+									}
+								}
+							}
+							combo_dsi[slice][nTile] =
+									(1-kx)*(1-ky)*d4[0][0] +
+									(  kx)*(1-ky)*d4[0][1] +
+									(1-kx)*(  ky)*d4[1][0] +
+									(  kx)*(  ky)*d4[1][1];
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+		return combo_dsi;
+	}
+	
+    public static float [][] getTDComboSceneSequence( // never used
+    		CLTParameters  clt_parameters,
+    		final boolean  save_weights, // output corresponding weights for each data
+    		final boolean  merge_all,
+    		final int      sensor_mask,
+    		Rectangle      fov_tiles,
+    		double []      stereo_xyz, // offset reference camera {x,y,z}
+    		double []      stereo_atr, // offset reference orientation (cuas)
+    		double []      ref_disparity,			
+    		QuadCLT []     quadCLTs,
+    		QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+//    		int            ref_index,
+    		int            debugLevel) {
+//    	double [] stereo_atr =        (stereo_atr_in != null)? stereo_atr_in: ZERO3; // maybe later play with rotated camera
+    	boolean mode_cuas =   (stereo_atr!=null) && ((stereo_atr[0] != 0) || (stereo_atr[1] != 0) || (stereo_atr[2] != 0));
+        double [][] ref_pXpYD = getRefPxPyD(
+        		clt_parameters, // CLTParameters  clt_parameters,
+        		mode_cuas,      // boolean        mode_cuas,            		
+//        		save_weights,   // final boolean  save_weights, // output corresponding weights for each data
+        		fov_tiles,      // Rectangle      fov_tiles,
+        		stereo_xyz,     // double []      stereo_xyz, // offset reference camera {x,y,z}
+        		stereo_atr,     // double []      stereo_atr_in, // offset reference orientation (cuas)
+        		ref_disparity,  // double []      ref_disparity,			
+        		refCLT,         // QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+        		debugLevel);    // int            debugLevel)
+        
+        return getTDComboSceneSequence(
+        		clt_parameters, // CLTParameters  clt_parameters,
+        		ref_pXpYD,      // double [][]    ref_pXpYD,    		
+        		save_weights,   // final boolean  save_weights, // output corresponding weights for each data
+        		merge_all,      // final boolean  merge_all,
+        		sensor_mask,    // final int      sensor_mask,
+        		fov_tiles,      // Rectangle      fov_tiles,
+        		stereo_xyz,     // double []      stereo_xyz, // offset reference camera {x,y,z}
+        		stereo_atr,     // double []      stereo_atr_in, // offset reference orientation (cuas)
+        		ref_disparity,  // double []      ref_disparity,			
+        		quadCLTs,       // QuadCLT []     quadCLTs,
+        		refCLT,         // QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+        		debugLevel);    // int            debugLevel);
+    }
+ 	
+    public static float [][] getTDComboSceneSequence(
+    		CLTParameters  clt_parameters,
+    		double [][]    ref_pXpYD,    		
+    		final boolean  save_weights, // output corresponding weights for each data
+    		final boolean  merge_all,
+    		final int      sensor_mask,
+    		Rectangle      fov_tiles,
+    		double []      stereo_xyz, // offset reference camera {x,y,z}
+    		double []      stereo_atr_in, // offset reference orientation (cuas)
+    		double []      ref_disparity,			
+    		QuadCLT []     quadCLTs,
+    		QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+    		int            debugLevel) {
+        int dbg_scene = -95;
+    	double [] stereo_atr = (stereo_atr_in != null)? stereo_atr_in: OpticalFlow.ZERO3; // maybe later play with rotated camera
+    	boolean mode_cuas =    (stereo_atr[0] != 0) || (stereo_atr[1] != 0) || (stereo_atr[2] != 0);
+        double [][] ref_pXpYD_or_null = mode_cuas ? ref_pXpYD : null; // debugging cuas mode keeping old
+		boolean mb_en =               clt_parameters.imp.mb_en && (fov_tiles==null);
+		double  mb_tau =              clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
+		double  mb_max_gain =         clt_parameters.imp.mb_max_gain; // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
+		ErsCorrection ers_reference = refCLT.getErsCorrection();
+        //  quadCLTs[0].getCltLength(boolean use_ref)
+        int sc0 = -1;
+		for (int nscene =  0; nscene < quadCLTs.length ; nscene++) if (quadCLTs[nscene] != null){
+			sc0 = nscene;
+			break;
+		}
+		
+		final int num_slices = merge_all? 1 : quadCLTs[sc0].getNumSensors();
+        final float [][] sumFclt = new float [num_slices * (save_weights? 2:1)][quadCLTs[sc0].getCltSize(false)];
+        final int [][] numAcc = new int [num_slices][sumFclt[0].length];
+        // next two to improve multithreading performance
+        final int tile_size_td = 4 * GPUTileProcessor.DTT_SIZE * GPUTileProcessor.DTT_SIZE;
+        final int tiles_td = sumFclt[0].length/tile_size_td;
+        final int tiles_td_all = tiles_td * num_slices; // usually sumFclt.length==1
+        
+		final Thread[] threads = ImageDtt.newThreadArray(ImageDtt.THREADS_MAX);
+		final AtomicInteger ai = new AtomicInteger(0);
+		for (int nscene =  0; nscene < quadCLTs.length ; nscene++) if (quadCLTs[nscene] != null){
+			if (nscene== dbg_scene) {
+				System.out.println("renderSceneSequence(): nscene = "+nscene);
+			}
+			String ts = quadCLTs[nscene].getImageName();
+			double []   scene_xyz = OpticalFlow.ZERO3;
+			double []   scene_atr = OpticalFlow.ZERO3;
+//			if (nscene != ref_index) { // Check even for raw, so video frames will match in all modes 
+			if (quadCLTs[nscene] != refCLT) { // Check even for raw, so video frames will match in all modes 
+				scene_xyz = ers_reference.getSceneXYZ(ts);
+				scene_atr = ers_reference.getSceneATR(ts);
+				if ((scene_atr==null) || (scene_xyz == null)) {
+					continue;
+				}
+				double []   scene_ers_xyz_dt = ers_reference.getSceneErsXYZ_dt(ts);
+				double []   scene_ers_atr_dt = ers_reference.getSceneErsATR_dt(ts);
+				quadCLTs[nscene].getErsCorrection().setErsDt(
+						scene_ers_xyz_dt, // double []    ers_xyz_dt,
+						scene_ers_atr_dt); // double []    ers_atr_dt)(ers_scene_original_xyz_dt);
+			}
+			if (!mode_cuas && (stereo_xyz != null)) { // offset all, including reference scene - now always, it is never null
+				double [][] combo_xyzatr = ErsCorrection.combineXYZATR(
+						stereo_xyz,  // double [] reference_xyz,
+						stereo_atr,  // double [] reference_atr, 
+						scene_xyz,   // double [] scene_xyz,
+						scene_atr);  // double [] scene_atr) 
+				scene_xyz = combo_xyzatr[0];
+				scene_atr = combo_xyzatr[1];
+			}
+			int sm = merge_all? -1: sensor_mask;
+			float [][] fclt = null;
+			double [][] dxyzatr_dt = null;
+			// should get velocities from HashMap at reference scene from timestamp , not re-calculate.
+			if (mb_en) {
+//				dxyzatr_dt = getVelocities(
+//						quadCLTs, // QuadCLT []     quadCLTs,
+//						nscene);  // int            nscene)
+				dxyzatr_dt = new double[][] { // for all, including ref
+					quadCLTs[nscene].getErsCorrection().getErsXYZ_dt(),
+					quadCLTs[nscene].getErsCorrection().getErsATR_dt()};				
+			}
+			
+			
+			if (mb_en && (dxyzatr_dt != null)) {
+				double [][] motion_blur = OpticalFlow.getMotionBlur(
+						refCLT, // quadCLTs[ref_index],   // QuadCLT        ref_scene,
+						quadCLTs[nscene],      // QuadCLT        scene,         // can be the same as ref_scene
+						ref_pXpYD,             // double [][]    ref_pXpYD,     // here it is scene, not reference!
+						scene_xyz,             // double []      camera_xyz,
+						scene_atr,             // double []      camera_atr,
+						dxyzatr_dt[0],         // double []      camera_xyz_dt,
+						dxyzatr_dt[1],         // double []      camera_atr_dt,
+						0,                     // int            shrink_gaps,  // will gaps, but not more that grow by this
+						debugLevel);           // int            debug_level)
+				fclt = QuadCLT.getTDCombo(
+						sm,                  // final int         sensor_mask,
+						merge_all,           // final boolean     merge_channels,
+						null,                // final Rectangle   full_woi_in,      // show larger than sensor WOI (or null)
+						clt_parameters,      // CLTParameters     clt_parameters,
+						ref_disparity,       // double []         disparity_ref,
+						ref_pXpYD_or_null,   // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
+						// motion blur compensation 
+						mb_tau,              // double            mb_tau,      // 0.008; // time constant, sec
+						mb_max_gain,         // double            mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
+						motion_blur,         // double [][]       mb_vectors,  //
+						scene_xyz,           // final double []   scene_xyz, // camera center in world coordinates
+						scene_atr,           // final double []   scene_atr, // camera orientation relative to world frame
+						quadCLTs[nscene],    // final QuadCLT     scene,
+						refCLT, // quadCLTs[ref_index], // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
+						true,                // final boolean     show_nan,
+						QuadCLT.THREADS_MAX, // int               threadsMax,
+						debugLevel);         // final int         debugLevel) 
+
+			} else {
+				fclt = QuadCLT.getTDCombo(
+						sm,                  // final int         sensor_mask,
+						merge_all,           // final boolean     merge_channels,
+						null,                // final Rectangle   full_woi_in,      // show larger than sensor WOI (or null)
+						clt_parameters,      // CLTParameters     clt_parameters,
+						ref_disparity,       // double []         disparity_ref,
+						ref_pXpYD_or_null,   // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
+						// motion blur compensation 
+						0,                   // double            mb_tau,      // 0.008; // time constant, sec
+						0,                   // double            mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
+						null,                // double [][]       mb_vectors,  //
+						scene_xyz,           // final double []   scene_xyz, // camera center in world coordinates
+						scene_atr,           // final double []   scene_atr, // camera orientation relative to world frame
+						quadCLTs[nscene],    // final QuadCLT     scene,
+						refCLT, // quadCLTs[ref_index], // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
+						true,                // final boolean     show_nan,
+						QuadCLT.THREADS_MAX, // int               threadsMax,
+						debugLevel);         // final int         debugLevel) 
+			}
+			final float [][] ffclt = fclt;
+			ai.set(0);
+			for (int ithread = 0; ithread < threads.length; ithread++) {
+				threads[ithread] = new Thread() {
+					public void run() {
+						for (int nTileAll = ai.getAndIncrement(); nTileAll < tiles_td_all; nTileAll = ai.getAndIncrement()) {
+							int ntile = nTileAll % tiles_td;
+							int nsens = nTileAll / tiles_td;
+							int indx0 = ntile * tile_size_td;
+							int indx1 = indx0 + tile_size_td;
+							for (int indx = indx0; indx < indx1; indx++) {
+								float d =ffclt[nsens][indx];
+								if (!Float.isNaN(d)){
+									sumFclt[nsens][indx] += d;
+									numAcc[nsens][indx]++ ;
+								}
+							}
+						}
+					}
+				};
+			}		      
+			ImageDtt.startAndJoin(threads);
+		}
+		// calculate averages
+		ai.set(0);
+		for (int ithread = 0; ithread < threads.length; ithread++) {
+			threads[ithread] = new Thread() {
+				public void run() {
+					for (int nTileAll = ai.getAndIncrement(); nTileAll < tiles_td_all; nTileAll = ai.getAndIncrement()) {
+						int ntile = nTileAll % tiles_td;
+						int nsens = nTileAll / tiles_td;
+						int indx0 = ntile * tile_size_td;
+						int indx1 = indx0 + tile_size_td;
+						for (int indx = indx0; indx < indx1; indx++) {
+							if (numAcc[nsens][indx] > 0) {
+								sumFclt[nsens][indx]/=numAcc[nsens][indx];
+							} else {
+								sumFclt[nsens][indx] = Float.NaN;
+							}
+							if (save_weights) {
+								sumFclt[nsens + num_slices][indx] = numAcc[nsens][indx];
+							}
+						}
+					}
+				}
+			};
+		}		      
+		ImageDtt.startAndJoin(threads);
+    	return sumFclt;
+    }
+    
+    public static double [][] getRefPxPyD(
+    		CLTParameters  clt_parameters,
+    		boolean        mode_cuas,
+    		Rectangle      fov_tiles,
+    		double []      stereo_xyz, // offset reference camera {x,y,z}
+    		double []      stereo_atr_in, // offset reference orientation (cuas)
+    		double []      ref_disparity,			
+    		QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+    		int            debugLevel) {
+    	double [] stereo_atr =        (stereo_atr_in != null)? stereo_atr_in: OpticalFlow.ZERO3; // maybe later play with rotated camera
+        double [][] ref_pXpYD;
+        if (mode_cuas) { //  && (dbg_scene > 0)) {
+			int         around =              2;
+			double      around_sigma =        4.0;
+			int         num_virtual_refines = 2;
+			String      debugSuffix=         null; // "virtual";
+			ref_pXpYD= transformFromVirtual(
+					ref_disparity, // double []         disparity_ref,
+					stereo_xyz, // final double []   scene_xyz,         // camera center in world (reference) coordinates
+					stereo_atr, // final double []   scene_atr,         // camera orientation relative to world (reference) frame
+					refCLT, // quadCLTs[ref_index], // final QuadCLT     reference_QuadClt,
+					around, // final int         around,            // 2 search around for interpolation
+					around_sigma, // final double      sigma,
+					num_virtual_refines, // final int         num_refines,
+					debugSuffix); // final String      debugSuffix)
+			refCLT.getErsCorrection().setupERS();
+			System.out.println("getRefPxPyD(): Calculated virtual_PxPyD");
+        } else {
+        	ref_pXpYD = OpticalFlow.transformToScenePxPyD( // now should work with offset ref_scene
+        			fov_tiles,            // final Rectangle [] extra_woi,    // show larger than sensor WOI (or null)
+        			ref_disparity,        // final double []   disparity_ref, // invalid tiles - NaN in disparity
+        			OpticalFlow.ZERO3, // stereo_xyz, // ZERO3,                // final double []   scene_xyz, // camera center in world coordinates
+        			OpticalFlow.ZERO3, // stereo_atr, // ZERO3,                // final double []   scene_atr, // camera orientation relative to world frame
+        			refCLT, // quadCLTs[ref_index],  // final QuadCLT     scene_QuadClt,
+        			refCLT, // quadCLTs[ref_index],  // final QuadCLT     reference_QuadClt, // now - may be null - for testing if scene is rotated ref
+        			QuadCLT.THREADS_MAX);          // int               threadsMax)
+        }
+        return ref_pXpYD;
+    }
+
+    
+	public static double [][] transformFromVirtual(
+			double []         disparity_ref,
+			final double []   scene_xyz,         // camera center in world (reference) coordinates
+			final double []   scene_atr,         // camera orientation relative to world (reference) frame
+			final QuadCLT     reference_QuadClt,
+			final int         around,            // 2 search around for interpolation
+			final double      sigma,
+			final int         num_refines,
+			final String      debugSuffix){
+		final boolean debug = debugSuffix != null;
+		final double sigma2 = 2*sigma*sigma;
+    	final double normal_damping = 0.001; // pull to horizontal if not enough data 
+    	final double [] damping = new double [] {normal_damping, normal_damping};
+		TileProcessor tp =         reference_QuadClt.getTileProcessor();
+		final int tilesX =         tp.getTilesX();
+		final int tilesY =         tp.getTilesY();
+		final int tiles =          tilesY*tilesX;
+		final int tilesX_around =  tilesX + 2 * around; 
+		final int tilesY_around =  tilesY + 2 * around; 
+		final int tiles_around =   tilesY_around * tilesX_around;
+		final int transform_size = tp.getTileSize();
+		final Thread[] threads = ImageDtt.newThreadArray(ImageDtt.THREADS_MAX);
+		final AtomicInteger ai = new AtomicInteger(0);
+		int initial_capacity = 4;
+		final int min_samples = 1;
+		
+		final String debugTitle=debug?reference_QuadClt.getImageName()+"-"+debugSuffix : null;
+		String [] debug_frame_titles = {"X", "Y", "D"};
+		String [] debug_titles = new String[num_refines+1];
+		for (int i = 0; i <= num_refines; i++) {
+			debug_titles[i]=i+"";
+		}
+		final double [][][] debug_data = debug ? new double [3][num_refines+1][tiles]:null;
+		if (debug) {
+			for (int f = 0; f < debug_data.length; f++) {
+				for (int i = 0; i < debug_data[f].length; i++) {
+					Arrays.fill(debug_data[f][i], Double.NaN);
+				}
+			}
+		}
+
+		final ArrayList<ArrayList<Integer>> interp_list = new ArrayList<ArrayList<Integer>>(tiles_around);
+		for (int nTile = 0; nTile < tiles_around; nTile++) {
+			interp_list.add(new ArrayList<Integer>(initial_capacity));
+		}
+		// create uniform grid for initial interpolations
+		final double [][] reference_pXpYD = new double [tiles][];
+		for (int nrefine = 0; nrefine <= num_refines; nrefine++) {
+			if (nrefine == 0) {
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) if (!Double.isNaN(disparity_ref[nTile])){
+								double disparity = disparity_ref[nTile];
+								int tileY = nTile / tilesX;  
+								int tileX = nTile % tilesX;
+								double centerX = (tileX + 0.5) * transform_size; //  - shiftX;
+								double centerY = (tileY + 0.5) * transform_size; //  - shiftY;
+								if (disparity < 0) {
+									disparity = 1.0* disparity; // 0.0;
+								}
+								reference_pXpYD[nTile] = new double[] {centerX, centerY, disparity};
+							}
+						}
+					};
+				}		      
+				ImageDtt.startAndJoin(threads);
+				
+			} else {
+				// get scene pXpYD corresponding to reference_pXpYD
+				final double [][] scene_pXpYD = OpticalFlow.transformToScenePxPyD(
+						reference_pXpYD,    // final double [][] reference_pXpYD,// invalid tiles - NaN in disparity. Should be no nulls, no NaN disparity
+						scene_xyz,          // final double []   scene_xyz,         // camera center in world (reference) coordinates
+						scene_atr,          // final double []   scene_atr,         // camera orientation relative to world (reference) frame
+						reference_QuadClt, // final QuadCLT     reference_QuadClt) //
+						null); // final QuadCLT     scene_QuadClt) //
+
+				// deep clone reference_pXpYD
+				final double [][] reference_pXpYD_next = new double [tiles][];
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							for (int nTile = ai.getAndIncrement(); nTile < tiles_around; nTile = ai.getAndIncrement()) {
+								interp_list.get(nTile).clear(); // clear each list
+							}
+						}
+					};
+				}		      
+				ImageDtt.startAndJoin(threads);
+				// not multithreaded 
+				for (int nTile = 0; nTile < tiles; nTile++) if (scene_pXpYD[nTile] != null) {
+					int tileX = (int) Math.floor (scene_pXpYD[nTile][0]/transform_size);
+					int tileY = (int) Math.floor (scene_pXpYD[nTile][1]/transform_size);
+					int tileX_around = tileX + around; 
+					int tileY_around = tileY + around;
+					if ((tileX_around >= 0) && (tileY_around >= 0) && (tileX_around < tilesX_around) && (tileY_around < tilesY_around)) {
+						int tile_around = tileX_around + tilesX_around * tileY_around;
+						interp_list.get(tile_around).add(nTile);
+					}
+				}
+				// interpolate reference_pXpYD_next 
+				ai.set(0);
+				for (int ithread = 0; ithread < threads.length; ithread++) {
+					threads[ithread] = new Thread() {
+						public void run() {
+							PolynomialApproximation pa = new PolynomialApproximation();		
+							for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) {
+								int tileY = nTile / tilesX;  
+								int tileX = nTile % tilesX;
+								double centerX = (tileX + 0.5) * transform_size; //  - shiftX;
+								double centerY = (tileY + 0.5) * transform_size; //  - shiftY;
+								int tileX_around = tileX + around; 
+								int tileY_around = tileY + around;
+								int num_samples = 0;
+								for (int dy = -around; dy <= around; dy++) {
+									for (int dx = -around; dx <= around; dx++) {
+										int tile_around = (tileX_around + dx) + tilesX_around * (tileY_around + dy);
+										num_samples+=interp_list.get(tile_around).size();
+									}
+								}
+								int mindx = 0;
+								if (num_samples >= min_samples) {
+									// inter/extrapolate with regularization (for few samples) and weights
+									double [][][] mdata = new double [num_samples][3][];
+									for (int dy = -around; dy <= around; dy++) {
+										for (int dx = -around; dx <= around; dx++) {
+											int tile_around = (tileX_around + dx) + tilesX_around * (tileY_around + dy);
+											for (int nst:interp_list.get(tile_around)) {
+												double [] ref_xyd = reference_pXpYD[nst];
+												double [] scene_xyd = scene_pXpYD[nst];
+												double x = scene_xyd[0]-centerX;
+												double y = scene_xyd[1]-centerY;
+												double w = Math.exp (-(x*x + y*y)/sigma2);
+												mdata[mindx][0] = new double [2];
+												mdata[mindx][0][0] =  x;
+												mdata[mindx][0][1] =  y;
+												mdata[mindx][1] = new double [2]; // [3];
+												mdata[mindx][1][0] =  ref_xyd[0]; 
+												mdata[mindx][1][1] =  ref_xyd[1]; 
+//												mdata[mindx][1][2] =  ref_xyd[2];
+												mdata[mindx][2] = new double [1];
+												mdata[mindx][2][0] =  w;
+												mindx++;
+											}
+										}
+									}
+									double[][] approx2d = pa.quadraticApproximation(
+											mdata,
+											true,       // boolean forceLinear,  // use linear approximation
+											damping,    // double [] damping, null OK
+											-1);        // debug level
+									double px = approx2d[0][2];
+									double py = approx2d[1][2];
+									// interpolate disparity from disparity_ref
+									double tx = px/transform_size-0.5;
+									double ty = py/transform_size-0.5;
+									int tx0 = (int) Math.floor(tx);
+									int ty0 = (int) Math.floor(ty);
+									int tx1 = tx0+1;
+									int ty1 = ty0+1;
+									double kx = tx-tx0;
+									double ky = ty-ty0;
+									if (tx0 < 0) {
+										tx0 = 0;
+										tx1 = tx0+1;
+										kx = 0;
+									} else if (tx1 >= tilesX) {
+										tx1 = tilesX - 1;
+										tx0 = tx1-1;
+										kx = 1.0;
+									}
+									
+									if (ty0 < 0) {
+										ty0 = 0;
+										ty1 = ty0+1;
+										ky = 0;
+									} else if (ty1 >= tilesY) {
+										ty1 = tilesY - 1;
+										ty0 = ty1-1;
+										ky = 1.0;
+									}
+									double d = 
+											(1-kx)*(1-ky)*disparity_ref[tx0 + tilesX*ty0] +
+											(  kx)*(1-ky)*disparity_ref[tx1 + tilesX*ty0] +
+											(1-kx)*(  ky)*disparity_ref[tx0 + tilesX*ty1] +
+											(  kx)*(  ky)*disparity_ref[tx1 + tilesX*ty1];
+									reference_pXpYD_next[nTile] = new double[] {px, py, d};
+								}
+							}
+						}
+					};
+				}		      
+				ImageDtt.startAndJoin(threads);
+				System.arraycopy(
+						reference_pXpYD_next,
+						0,
+						reference_pXpYD,
+						0,
+						tiles);
+				// TODO: refine disparity from disparity_ref interpolation
+			} // if (nrefine == 0) else
+			if (debug) {
+				for (int nTile = 0; nTile < tiles; nTile++) if (reference_pXpYD[nTile] != null) {
+					debug_data[0][nrefine][nTile] = reference_pXpYD[nTile][0];
+					debug_data[1][nrefine][nTile] = reference_pXpYD[nTile][1];
+					debug_data[2][nrefine][nTile] = reference_pXpYD[nTile][2];
+				}
+			}
+		}
+		if (debug) {
+			ShowDoubleFloatArrays.showArraysHyperstack(
+					debug_data,         // double[][][] pixels, 
+					tilesX,             // int          width, 
+					debugTitle,         // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
+					debug_titles,       // String []    titles, // all slices*frames titles or just slice titles or null
+					debug_frame_titles, // String []    frame_titles, // frame titles or null
+					true);              // boolean      show)
+			
+		}		
+		
+		return reference_pXpYD; // ref_pXpYD_0; // 
+	}
 
 }

src/main/java/com/elphel/imagej/cuas/CuasCenterLma.java

@@ -84,6 +84,44 @@ public class CuasCenterLma {
 		return rslt;
 	}
 	
+	public static double [][] getCenterATR(
+			double [][] scenes_atr,
+			int []     range,
+			int        debugLevel) {
+		boolean [] param_select = new boolean[PARAMETER_NAMES.length];
+		Arrays.fill(param_select, true);
+		CuasCenterLma cuasCenterLma = new CuasCenterLma(
+				param_select, // boolean [] param_select,
+				scenes_atr,     // double [][] scenes_atr,
+				range,        // int []     range,
+				debugLevel);  // int debugLevel)
+		double lambda =             0.1;
+		double lambda_scale_good =  0.5;
+		double lambda_scale_bad =   8.0;
+		double lambda_max =       100;
+		double rms_diff =           0.001;
+		int    num_iter =          20;
+		int lmaResult = cuasCenterLma.runLma(
+				lambda,           // double lambda,           // 0.1
+				lambda_scale_good,// double lambda_scale_good,// 0.5
+				lambda_scale_bad, // double lambda_scale_bad, // 8.0
+				lambda_max,       // double lambda_max,       // 100
+				rms_diff,         // double rms_diff,         // 0.001
+				num_iter,         // int    num_iter,         // 20
+				debugLevel);      // int    debug_level)
+		double [][] rslt = {cuasCenterLma.getCenter(),cuasCenterLma.getRadius()};
+		if (debugLevel > -3) {
+			System.out.println("lmaResult =" +lmaResult+" iterations, RMSE ="+
+					cuasCenterLma.getRMS()+" ("+cuasCenterLma.getInitialRMS()+")");
+			System.out.println("azimuth_center = "+ rslt[0][0]);
+			System.out.println("   tilt_center = "+ rslt[0][1]);
+			System.out.println("  average roll = "+ rslt[0][2]);
+			System.out.println("azimuth_radius = "+ rslt[1][0]);
+			System.out.println("   tilt_radius = "+ rslt[1][1]);
+		}
+		return rslt;
+	}
+	
 	public CuasCenterLma(
 			boolean [] param_select,
 			QuadCLT [] quadCLTs,
@@ -98,6 +136,21 @@ public class CuasCenterLma {
 				range,        // int []     range,
 				debugLevel);  //int debugLevel)
 	}
+
+	public CuasCenterLma(
+			boolean [] param_select,
+			double [][] scenes_atr,
+			int []     range,
+			int        debugLevel) {
+		prepareLMA(
+				param_select, // boolean [] param_select,
+				scenes_atr,     // double [][] scenes_atr,
+				range,        // int []     range,
+				debugLevel);  //int debugLevel)
+	}
+	
+	
+	
 	public int prepareLMA(
 			boolean [] param_select,
 			QuadCLT [] quadCLTs,
@@ -121,6 +174,29 @@ public class CuasCenterLma {
         	String ts = quadCLTs[nscene].getImageName();
         	scenes_atr[nscene] = ers_reference.getSceneATR(ts);
         }
+        return prepareLMA(
+        		param_select, // boolean [] param_select,
+        		scenes_atr, // double [][] scenes_atr,
+    			new int [] {earliest_scene,last_scene}, // int []     range,
+    			debugLevel); // int        debugLevel);
+	}
+	
+	public int prepareLMA(
+			boolean [] param_select,
+			double [][] scenes_atr,
+			int []     range,
+			int        debugLevel) {
+		earliest_scene = range[0];
+		last_scene =     range[1];
+        int num_scenes = last_scene- earliest_scene + 1;
+        for (int nscene = last_scene; nscene >= earliest_scene; nscene--) {
+			// just checking it is not isolated
+        	if (scenes_atr[nscene] == null) {
+        		earliest_scene = nscene + 1;
+        		break;
+        	}
+        }
+        
 		full_parameters_vector = new double[PARAM_PHASE + num_scenes];
 		sel_par_rindex = new int[full_parameters_vector.length];
 		Arrays.fill(full_parameters_vector, Double.NaN);

src/main/java/com/elphel/imagej/gpu/GpuQuad.java

@@ -1891,14 +1891,24 @@ public class GpuQuad{ // quad camera description
 			return;
 		}
 		
-		
-		if (wh == null) {
-			wh = new int[] {img_width, img_height};
-		}
 		if (!rectilinear) {
 			setConvolutionKernels(false); // set kernels if they are not set already
 			setBayerImages(false); // set Bayer images if this.quadCLT instance has new ones
 		}
+
+		int [] wh1 = handleWH(
+				wh,         // int []  wh_in,
+				ref_scene); // boolean ref_scene);
+		if (wh == null) { // so wh will not change if was not null
+			wh = wh1;
+		}
+		int tilesX =  wh[0] / GPUTileProcessor.DTT_SIZE;
+		int tilesY =  wh[1] / GPUTileProcessor.DTT_SIZE;
+
+		/*
+		if (wh == null) {
+			wh = new int[] {img_width, img_height};
+		}
 		// kernel parameters: pointer to pointers
 		int tilesX =  wh[0] / GPUTileProcessor.DTT_SIZE;
 		int tilesY =  wh[1] / GPUTileProcessor.DTT_SIZE;
@@ -1954,6 +1964,8 @@ public class GpuQuad{ // quad camera description
 				gpu_clt_wh = wh.clone();
 			}
 		}
+		*/
+		
 		CUdeviceptr gpu_clt_selected = ref_scene ? gpu_clt_ref : gpu_clt;
 		int [] GridFullWarps =    {1, 1, 1};
 		int [] ThreadsFullWarps = {1, 1, 1};
@@ -2036,7 +2048,80 @@ public class GpuQuad{ // quad camera description
 			System.out.println("======execConvertDirect("+ref_scene+", "+erase_clt+")");
 		}
 	}
-
+	
+	/**
+	 * Handle image width and height, allocate/re-allocate GPU memory if size changed
+	 * @param wh     input {width, height}, use sensor dimensions if null. 
+	 * @param ref_scene use reference buffer, if false - main scene buffer
+	 * @return updated width, height pair, non-null
+	 */
+	public int [] handleWH(
+			int []  wh,
+			boolean ref_scene) {
+		
+//		int [] wh = (wh_in != null) ? wh_in :  (new int[] {img_width, img_height});
+		if (wh == null) {
+			wh = new int[] {img_width, img_height};
+		}
+		// kernel parameters: pointer to pointers
+		int tilesX =  wh[0] / GPUTileProcessor.DTT_SIZE;
+		int tilesY =  wh[1] / GPUTileProcessor.DTT_SIZE;
+		// De-allocate if size mismatch, allocate if needed. Now it is the only place where clt is allocated
+		if (ref_scene) {
+			if ((gpu_clt_ref_wh != null) && ((gpu_clt_ref_wh[0] != wh[0]) || (gpu_clt_ref_wh[1] != wh[1]))) {
+				for (int ncam = 0; ncam < num_cams; ncam++) {
+					cuMemFree (gpu_clt_ref_h[ncam]);
+				}
+				cuMemFree (gpu_clt_ref);
+				gpu_clt_ref = null;
+				gpu_clt_ref_wh = null;
+			}
+			if (gpu_clt_ref == null) { // Allocate memory, create pointers for reference scene TD representation
+				long [] gpu_clt_ref_l = new long [num_cams];
+				gpu_clt_ref_h =    new CUdeviceptr[num_cams];
+				for (int ncam = 0; ncam < num_cams; ncam++) {
+					gpu_clt_ref_h[ncam] = new CUdeviceptr();
+					cuMemAlloc(gpu_clt_ref_h[ncam],
+							   tilesY * tilesX * num_colors * 4 * GPUTileProcessor.DTT_SIZE * GPUTileProcessor.DTT_SIZE * Sizeof.FLOAT );
+				}
+				gpu_clt_ref =  new CUdeviceptr();
+				cuMemAlloc(gpu_clt_ref, num_cams * Sizeof.POINTER);
+				for (int ncam = 0; ncam < num_cams; ncam++) {
+					gpu_clt_ref_l[ncam] = GPUTileProcessor.getPointerAddress(gpu_clt_ref_h[ncam]);
+				}
+				cuMemcpyHtoD(gpu_clt_ref, Pointer.to(gpu_clt_ref_l), num_cams * Sizeof.POINTER);
+				gpu_clt_ref_wh = wh.clone();
+			}
+		} else { // same for main (not ref) memory
+			if ((gpu_clt_wh != null) && ((gpu_clt_wh[0] != wh[0]) || (gpu_clt_wh[1] != wh[1]))) {
+				for (int ncam = 0; ncam < num_cams; ncam++) {
+					cuMemFree (gpu_clt_h[ncam]);
+				}
+				cuMemFree (gpu_clt);
+				gpu_clt = null;
+				gpu_clt_wh = null;
+			}
+			if (gpu_clt == null) { // Allocate memory, create pointers for reference scene TD representation
+				long [] gpu_clt_l = new long [num_cams];
+				gpu_clt_h =    new CUdeviceptr[num_cams];
+				for (int ncam = 0; ncam < num_cams; ncam++) {
+					gpu_clt_h[ncam] = new CUdeviceptr();
+					cuMemAlloc(gpu_clt_h[ncam],
+							   tilesY * tilesX * num_colors * 4 * GPUTileProcessor.DTT_SIZE * GPUTileProcessor.DTT_SIZE * Sizeof.FLOAT );
+				}
+				gpu_clt =  new CUdeviceptr();
+				cuMemAlloc(gpu_clt, num_cams * Sizeof.POINTER);
+				for (int ncam = 0; ncam < num_cams; ncam++) {
+					gpu_clt_l[ncam] = GPUTileProcessor.getPointerAddress(gpu_clt_h[ncam]);
+				}
+				cuMemcpyHtoD(gpu_clt, Pointer.to(gpu_clt_l), num_cams * Sizeof.POINTER);
+				gpu_clt_wh = wh.clone();
+			}
+		}
+		return wh;
+	}
+	
+	
 	public boolean reAllocateClt(
 			int []  wh,
 			boolean ref_scene) {

src/main/java/com/elphel/imagej/tileprocessor/ImageDtt.java

@@ -1339,6 +1339,7 @@ public class ImageDtt extends ImageDttCPU {
 	/**
 	 * Convert reference scene to FD and save result in extra GPU array for the future interscene correlation
 	 * Geometry correction and images will come from gpuQuad instance - 
+	 * @param fclt load TD data and bypass execConvertDirect. If fclt.length==1 load same data to all sensor channels
 	 * @param erase_clt erase CLT (<0 - do not erase, 0 - erase to 0.0, >0 - erase to NaN). Needed only for later IMCLT
 	 *                  end rendering images. NaN produces sharp, distinct borders; 0f - blended
 	 * @param wh if null, will uses sensor dimensions. Otherwise {width, height} in pixels
@@ -1353,6 +1354,7 @@ public class ImageDtt extends ImageDttCPU {
 	 * @param globalDebugLevel
 	 */
 	public void setReferenceTD(
+			final float [][]          fclt, 
 			final int                 erase_clt,
 			final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 			final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
@@ -1384,7 +1386,26 @@ public class ImageDtt extends ImageDttCPU {
 		gpuQuad.updateTasks(
 				tp_tasks,
 				false); // boolean use_aux    // while is it in class member? - just to be able to free
-		gpuQuad.execConvertDirect(use_reference_buffer, wh, erase_clt); // put results into a "reference" buffer
+		if (fclt != null) {
+			gpuQuad.handleWH( // allocate/reallocate GPU memory, that was normally done by gpuQuad.execConvertDirect()
+					  wh, // int []  wh,
+					  use_reference_buffer); // boolean ref_scene)
+			int sensor_mask = -1;
+			int num_cams = getNumSensors();
+			boolean merge_channels = (fclt.length == 1);
+			for (int ncam = 0; ncam < num_cams; ncam++) if (((1 << ncam) & sensor_mask) != 0){
+				int src_cam = merge_channels? 0: ncam;
+				if (src_cam < fclt.length) {
+					gpuQuad.setCltData( // for testing only
+							ncam, // int      ncam,
+							fclt[src_cam], // float [] fclt, //
+							true); // boolean  use_ref);
+				}
+			}
+		} else {
+			gpuQuad.execConvertDirect(use_reference_buffer, wh, erase_clt); // put results into a "reference" buffer
+		}
+		
 	}
 
 	public void setRectilinearReferenceTD(

src/main/java/com/elphel/imagej/tileprocessor/Interscene.java

@@ -38,6 +38,7 @@ import com.elphel.imagej.cameras.ColorProcParameters;
 import com.elphel.imagej.cameras.EyesisCorrectionParameters;
 import com.elphel.imagej.common.DoubleGaussianBlur;
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
+import com.elphel.imagej.cuas.CuasCenterLma;
 import com.elphel.imagej.gpu.GpuQuad;
 import com.elphel.imagej.gpu.TpTask;
 import com.elphel.imagej.ims.Did_ins_1;
@@ -182,12 +183,15 @@ public class Interscene {
 			int adjusted_scene_index =  ref_index - 1;
 			for (; adjusted_scene_index >= 0 ; adjusted_scene_index--) {
 				// to include ref scene photometric calibration
-				quadCLTs[adjusted_scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT(
-						set_channels[adjusted_scene_index].set_name,
-						clt_parameters,
-						colorProcParameters, //
-						threadsMax,
-						debugLevel-2);
+				// added 05.21.2024: skip if was spawn already
+				if (quadCLTs[adjusted_scene_index] == null) {
+					quadCLTs[adjusted_scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT(
+							set_channels[adjusted_scene_index].set_name,
+							clt_parameters,
+							colorProcParameters, //
+							threadsMax,
+							debugLevel-2);
+				}
 				String ts = quadCLTs[adjusted_scene_index].getImageName();
 				if ((ts == null) || (ers_reference.getSceneXYZATR(ts) == null)){
 					break;
@@ -905,12 +909,15 @@ public class Interscene {
 		int []    fail_reason = new int[1];   // null or int[1]: 0 - OK, 2 - LMA, 3 - min, 4 - max
 		for (int scene_index =  ref_index - 1; scene_index >= 0 ; scene_index--) {
 			// to include ref scene photometric calibration
-			quadCLTs[scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT(
-					set_channels[scene_index].set_name,
-					clt_parameters,
-					colorProcParameters, //
-					threadsMax,
-					debugLevel-2);
+			// added 05.21.2024: skip if was spawn already
+			if (quadCLTs[scene_index] == null) { 
+				quadCLTs[scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT(
+						set_channels[scene_index].set_name,
+						clt_parameters,
+						colorProcParameters, //
+						threadsMax,
+						debugLevel-2);
+			}
 		} // split cycles to remove output clutter
 		
 		
@@ -972,7 +979,7 @@ public class Interscene {
 		if (debugLevel > -3) {
 			System.out.println("setInitialOrientationsIms(): compensate_ims_rotation="+compensate_ims_rotation);
 			System.out.println("setInitialOrientationsIms(): will "+((quat_corr==null)? "NOT ":"")+" correct orientation offset");
-			QuadCLTCPU.showQuatCorr(quat_corr, enu_corr_metric);
+			QuadCLTCPU.showQuatCorr(quat_corr, enu_corr_metric); // enu_corr_metric - linear correction
 		}
 		boolean test_quat_corr = debugLevel > 1000;
 		if (test_quat_corr) {
@@ -983,21 +990,6 @@ public class Interscene {
 					debugLevel);    // final int           debugLevel);
 		}
 		double [][][] ims_xyzatr;
-		
-		// TODO: Remove when done!
-//		if ((nscene == 57) || (nscene == 58)) {
-			
-//		}
-		/*
-		double [][][] ims_xyzatr_PIMU= QuadCLT.integratePIMU(
-				clt_parameters, // final CLTParameters clt_parameters,
-				quadCLTs,       // final QuadCLT[]     quadCLTs,
-				ref_index,      // final int           ref_index,
-				null,           // double [][][]       dxyzatr,
-	    		0, // final int           early_index,
-	    		(quadCLTs.length -1) // int           last_index,
-	    		);
-       */
 		if (inertial_only) {
 			ims_xyzatr = QuadCLT.integratePIMU(
 						clt_parameters, // final CLTParameters clt_parameters,
@@ -1005,9 +997,10 @@ public class Interscene {
 						ref_index,      // final int           ref_index,
 						null,           // double [][][]       dxyzatr,
 			    		0, // final int           early_index,
-			    		(quadCLTs.length -1) // int           last_index,
+			    		(quadCLTs.length -1), // int           last_index,
 //			    		(new double[2][3]) //final double [][] pimu_offsets // null - will not subtract from velocities
-			    		);
+						debugLevel      // final int           debugLevel
+					);
 		} else {
 			ims_xyzatr = quadCLTs[ref_index].getXyzatrIms(
 					clt_parameters, // CLTParameters clt_parameters,
@@ -1431,6 +1424,438 @@ public class Interscene {
 		return earliest_scene;
 	}
 
+	public static int setInitialOrientationsCuas(
+			final CLTParameters          clt_parameters,
+			final double [][]            center_xyzatr, // center_xyzatr[0]=ZERO3 - for this sequence
+			final QuadCLT                center_CLT,    // contains center CLT and DSI
+			final boolean                compensate_ims_rotation,
+			final boolean                inertial_only, 			
+			int                          min_num_scenes,
+			final ColorProcParameters    colorProcParameters,
+			final QuadCLT[]              quadCLTs, //
+			final QuadCLT.SetChannels [] set_channels,
+			final boolean                batch_mode,
+			int []                       start_ref_pointers, // [0] - earliest valid scene, [1] ref_index. Set if non-null
+			final boolean                updateStatus,
+			final int                    debugLevel) {
+		final int last_index =      start_ref_pointers[1];
+		final int earliest_scene = start_ref_pointers[0];
+		RMSEStats rmse_stats = new RMSEStats();
+		RMSEStats rmse_stats_metric = clt_parameters.imp.eig_use ? (new RMSEStats()):null;
+		double  min_ref_str =         clt_parameters.imp.min_ref_str;
+		boolean ref_need_lma =        clt_parameters.imp.ref_need_lma;
+		double  min_ref_frac=         clt_parameters.imp.min_ref_frac;
+
+		boolean ref_smooth =          clt_parameters.imp.ref_smooth;  //true;   // use SfM filtering if available 
+		boolean lock_position =       true; // clt_parameters.imp.lock_position;
+    	boolean ims_use =             clt_parameters.imp.ims_use;
+    	
+    	if (!ims_use) {
+    		System.out.println("setInitialOrientationsCuas(): IMS use is disabled");
+    		return -1;
+    	}
+		if (clt_parameters.ilp.ilma_3d) {
+			System.out.println("*** setInitialOrientationsCuas(): clt_parameters.ilp.ilma_3d is TRUE. You may want to disable it when using IMS");
+		}
+
+		double [] lma_rms = new double[clt_parameters.imp.eig_use? 6 : 4]; // [2];
+
+		int tilesX =      center_CLT.getTileProcessor().getTilesX();
+//		int tilesY =      center_CLT.getTileProcessor().getTilesY();
+        int tile_size =   center_CLT.getTileProcessor().getTileSize();
+		double min_offset =         clt_parameters.imp.min_offset; // 0.0; // 
+		double max_offset =         clt_parameters.imp.max_rel_offset * tilesX * tile_size;
+		double max_roll =           clt_parameters.imp.max_roll_deg*Math.PI/180.0;
+		double max_zoom_diff =      clt_parameters.imp.max_zoom_diff;
+		boolean    lma_use_Z =      clt_parameters.imp.lma_use_Z;          // true;       // lpf x and y, re-adjust X,Y,Z,A,T,R with pull for X,Y. Disables  
+		boolean    lma_use_R =      clt_parameters.imp.lma_use_R;          // true;       // lpf x and y, re-adjust X,Y,Z,A,T,R with pull for X,Y. Disables  
+		if (debugLevel > -3) {
+			System.out.println("setInitialOrientationsCuas(): lock_position="+lock_position+", lma_use_Z="+lma_use_Z+", lma_use_R="+lma_use_R);
+	 	}
+		boolean[]    param_select = 
+				ErsCorrection.getParamSelect( // ZR - always
+						lma_use_Z && !lock_position, // boolean use_Z,
+						lma_use_R, // boolean use_R,
+						!lock_position, // true,      // boolean use_XY, May change with oorient preference
+						lock_position, // false,     // boolean use_AT, May change with oorient preference
+						false,     // boolean use_ERS);//clt_parameters.ilp.ilma_lma_select;
+						false,     // boolean use_ERS_tilt);
+						false);    // boolean use_ERS_roll);
+		
+		int []    fail_reason = new int[1];   // null or int[1]: 0 - OK, 2 - LMA, 3 - min, 4 - max
+		for (int scene_index =  last_index; scene_index >= 0 ; scene_index--) {
+			// to include ref scene photometric calibration
+			// added 05.21.2024: skip if was spawn already
+			if (quadCLTs[scene_index] == null) { 
+//				quadCLTs[scene_index] = quadCLTs[last_index].spawnNoModelQuadCLT(
+				quadCLTs[scene_index] = center_CLT.spawnNoModelQuadCLT(
+						set_channels[scene_index].set_name,
+						clt_parameters,
+						colorProcParameters, //
+						ImageDtt.THREADS_MAX, // threadsMax,
+						debugLevel-2);
+			}
+		}
+//		ErsCorrection ers_reference = quadCLTs[ref_index].getErsCorrection();		
+		ErsCorrection ers_center =    center_CLT.getErsCorrection();		
+		int debug_scene =            -57; // -15; last before
+		boolean debug2 =             !batch_mode; // false; // true;
+		boolean [] reliable_ref =     null;
+		boolean use_lma_dsi =         clt_parameters.imp.use_lma_dsi;
+		double [][] center_dsi =      center_CLT.dsi;
+        double [] reduced_strength =  new double[1];
+        // use combo if second pass?
+		boolean is_aux = center_CLT.isAux();
+		boolean [] blue_sky =     QuadCLT.getBooleanBlueSky(center_dsi[is_aux?TwoQuadCLT.DSI_BLUE_SKY_AUX:TwoQuadCLT.DSI_BLUE_SKY_MAIN]);
+    	double [] disparity =     center_dsi[is_aux?TwoQuadCLT.DSI_DISPARITY_AUX:TwoQuadCLT.DSI_DISPARITY_MAIN];
+    	double [] disparity_lma = center_dsi[is_aux?TwoQuadCLT.DSI_DISPARITY_AUX_LMA:TwoQuadCLT.DSI_DISPARITY_MAIN_LMA];
+    	double [] strength =      center_dsi[is_aux?TwoQuadCLT.DSI_STRENGTH_AUX:TwoQuadCLT.DSI_STRENGTH_MAIN];
+    	double [] center_disparity = use_lma_dsi? disparity_lma:disparity; //CHANGE?
+		if (min_ref_str > 0.0) {
+			reliable_ref = QuadCLT.getReliableTiles(
+					min_ref_str,     // double     min_strength,
+	    			min_ref_frac,     // double     min_ref_frac,
+	    			ref_need_lma,        // boolean    needs_lma,
+	    			false,       // boolean    sfm_filter, 
+					0, // sfm_minmax,  // double  sfm_minmax, // minimal value of the SfM gain maximum to consider available
+					0, // sfm_fracmax, // double  sfm_fracmax,// minimal fraction of the SfM maximal gain
+					0, // sfm_fracall, // double  sfm_fracall,// minimal relative area of the SfM-enabled tiles (do not apply filter if less)						
+	    			reduced_strength, // double []  reduced_strength, // if not null will return >0 if had to reduce strength (no change if did not reduce)
+	    			disparity_lma,    // double []  disparity_lma,
+	    			strength,         // double []  strength,
+	    			null, // sfm_gain,         // double []  sfm_gain,
+	    			blue_sky,         // boolean [] blue_sky,
+	    			debugLevel);      // int        debugLevel)
+			if (reduced_strength[0] > 0) {
+				use_lma_dsi = false; // too few points
+			}
+			if (debug2) {
+				double [] dbg_img = new double [reliable_ref.length]; // null
+				for (int i = 0; i < dbg_img.length; i++) {
+					dbg_img[i] = reliable_ref[i]?1:0;
+				}
+				ShowDoubleFloatArrays.showArrays(
+						dbg_img,
+						center_CLT.getTileProcessor().getTilesX(),
+						center_CLT.getTileProcessor().getTilesY(),
+						"reliable_center_ref");
+			}
+		}
+		
+		double max_z_change = Double.NaN; // only applicable for drone images
+		double [][][] scenes_dxyzatr = QuadCLTCPU.getDxyzatrPIMU( // does not depend on correction
+				clt_parameters, // CLTParameters clt_parameters,
+				quadCLTs); //,       // QuadCLT[]     quadCLTs,
+		//*********** Get IMS correction from reference scene *******************
+		double [] quat_corr = compensate_ims_rotation? quadCLTs[last_index].getQuatCorr() : null; //
+///		double [] quat_corr = compensate_ims_rotation? center_CLT.getQuatCorr() : null; //
+		double [] enu_corr_metric = quadCLTs[last_index].getENUCorrMetric();
+///		double [] enu_corr_metric = center_CLT.getENUCorrMetric();
+		// quat_corr may still be null if does not exist
+		
+		if (debugLevel > -3) {
+			System.out.println("setInitialOrientationsCuas(): compensate_ims_rotation="+compensate_ims_rotation);
+			System.out.println("setInitialOrientationsCuas(): will "+((quat_corr==null)? "NOT ":"")+" correct orientation offset");
+			QuadCLTCPU.showQuatCorr(quat_corr, enu_corr_metric); // enu_corr_metric - linear correction
+		}
+
+		double [][][] ims_xyzatr;
+		if (inertial_only) {
+			ims_xyzatr = QuadCLT.integratePIMU(
+						clt_parameters, // final CLTParameters clt_parameters,
+						quadCLTs,       // final QuadCLT[]     quadCLTs,
+//						ref_index,      // final int           ref_index,
+						last_index,      // final int           ref_index,
+						null,           // double [][][]       dxyzatr,
+			    		0, // final int           early_index,
+			    		(quadCLTs.length -1), // int           last_index,
+						debugLevel      // final int           debugLevel
+					);
+		} else {
+			ims_xyzatr = quadCLTs[last_index].getXyzatrIms(
+///				ims_xyzatr = center_CLT.getXyzatrIms( // use last index
+					clt_parameters, // CLTParameters clt_parameters,
+					quadCLTs, // QuadCLT[]  quadCLTs,
+					quat_corr, // double []  quat_corr, // only applies to rotations - verify!
+					debugLevel) ; // int        debugLevel)
+		}
+		if (true) { //  lock_position) { // zero linear offsets
+			for (int i = 0; i < quadCLTs.length; i++) for (int j = 0; j < 3; j++){
+				ims_xyzatr[i][0][j]=0.0;
+			}
+		}
+		// re-calculate center for this sequence
+		double [][] ims_atr = new double [ims_xyzatr.length][];
+		for (int i = 0; i < quadCLTs.length; i++) if (ims_xyzatr[i] != null) {
+			ims_atr[i] = ims_xyzatr[i][1];
+		}
+		int [] range = {0, last_index};
+		double [] center_atr =  CuasCenterLma.getCenterATR(
+				ims_atr,         // double [][] ims_atr,
+			     range,          // int []     range,
+			     debugLevel)[0]; // int        debugLevel)
+		
+		double [][][] scenes_xyzatr =      new double [quadCLTs.length][][]; // previous scene relative to the next one
+////		scenes_xyzatr[last_index] =         new double[2][3]; // all zeros
+		double [][] last_corr_xyzatr = {ZERO3,{-center_atr[0],-center_atr[1],-center_atr[2]}};
+
+		int margin =          clt_parameters.imp.margin;
+		boolean mb_en =       clt_parameters.imp.mb_en;
+		double  mb_tau =      clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
+		double  mb_max_gain = clt_parameters.imp.mb_max_gain_inter; // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
+		double [] reg_weights = clt_parameters.ilp.ilma_regularization_weights;
+		final double max_rms =     clt_parameters.imp.eig_use? clt_parameters.imp.eig_max_rms: clt_parameters.imp.max_rms;
+		float [][] fclt = new float [][] {center_CLT.getCenterClt()};
+		// set center tasks and CLT
+		double [][]  pXpYD_center = OpticalFlow.transformToScenePxPyD( // full size - [tilesX*tilesY], some nulls
+				null,               // final Rectangle [] extra_woi,    // show larger than sensor WOI (or null)
+				center_disparity,   // dls[0],  // final double []   disparity_ref, // invalid tiles - NaN in disparity (maybe it should not be masked by margins?)
+				ZERO3,              // final double []   scene_xyz, // camera center in world coordinates
+				ZERO3,              // final double []   scene_atr, // camera orientation relative to world frame
+				center_CLT,         // final QuadCLT     scene_QuadClt,
+				center_CLT);        // final QuadCLT     reference_QuadClt)
+		TpTask [][] tp_tasks_center = setReferenceGPU (
+				clt_parameters,     // CLTParameters      clt_parameters,			
+				center_CLT,         // QuadCLT            ref_scene,
+				center_disparity,   // null      // double []          ref_disparity, // null or alternative reference disparity
+				pXpYD_center,       // double [][]        ref_pXpYD,
+				fclt, 			    // final float [][]          fclt, 
+				reliable_ref,       // final boolean []   selection, // may be null, if not null do not  process unselected tiles
+				margin,             // final int          margin,
+				// motion blur compensation 
+				mb_tau,             // double             mb_tau,      // 0.008; // time constant, sec
+				mb_max_gain,        // double             mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
+				null,               // double [][]        mb_vectors,  // now [2][ntiles]; // center does not move
+				debugLevel);        // int                debug_level)
+		
+		for (int scene_index =  last_index ; scene_index >= earliest_scene ; scene_index--) {
+			if (scene_index == debug_scene) {
+				System.out.println("scene_index = "+scene_index);
+				System.out.println("scene_index = "+scene_index);
+			}
+			QuadCLT scene_QuadClt = quadCLTs[scene_index];
+			double [][] pose_ims = ims_xyzatr[scene_index];
+			
+			// Apply correction from the earlier processed scenes (closer to reference)
+			double [][] initial_pose = ErsCorrection.combineXYZATR(
+					last_corr_xyzatr, // current correction
+					pose_ims);        // raw IMS prediction
+			scenes_xyzatr[scene_index] = new double [][]{
+					initial_pose[0].clone(),
+					initial_pose[1].clone()
+			};
+			// Refine with LMA
+			//boost_zoom_short
+			double max_z_change_scaled = max_z_change;
+			double [] min_max = {min_offset, max_offset, 0.0} ; // {min, max, actual rms)
+			// Will be used in prepareLMA()
+			quadCLTs[scene_index].getErsCorrection().setErsDt( // set for ref also (should be set before non-ref!)
+					QuadCLTCPU.scaleDtToErs(
+							clt_parameters,
+							scenes_dxyzatr[scene_index])); // cam_dxyzatr));		
+// Trying new version with motion blur and single-setting of the reference frame
+			boolean adjust_OK = false;
+			boolean disable_ers = true;
+			double [][] mb_vectors_scene1 = null;
+			if (mb_en) {
+				double [][] scene1_xyzatr_dt = QuadCLTCPU.scaleDtFromErs(
+						clt_parameters,
+						quadCLTs[scene_index].getErsCorrection().getErsXYZATR_dt());
+				
+				mb_vectors_scene1 = OpticalFlow.getMotionBlur(
+						center_CLT,            // QuadCLT        ref_scene,
+						quadCLTs[scene_index], // QuadCLT        scene,         // can be the same as ref_scene
+						pXpYD_center,          // double [][]    ref_pXpYD,     // here it is scene, not reference!
+						initial_pose[0],       // scene1_xyzatr[0],     // double []      camera_xyz,
+						initial_pose[1],       // double []      camera_atr,
+						scene1_xyzatr_dt[0],   // double []      camera_xyz_dt,
+						scene1_xyzatr_dt[1],   // double []      camera_atr_dt,
+						0,                     // int            shrink_gaps,  // will gaps, but not more that grow by this
+						debugLevel);           // int            debug_level)
+			}
+			
+			scenes_xyzatr[scene_index] = adjustPairsLMAInterscene( // assumes reference GPU data is already set - once for multiple scenes
+					clt_parameters,        // CLTParameters  clt_parameters,
+					true,                  // boolean        initial_adjust,
+					false,                 // boolean        fpn_disable,   // disable fpn filter if images are known to be too close
+					disable_ers,           // boolean        disable_ers,
+					min_max,               // double []      min_max,       // null or pair of minimal and maximal offsets
+					fail_reason,           // int []         fail_reason,   // null or int[1]: 0 - OK, 1 - LMA, 2 - min, 3 - max
+					center_CLT,            // QuadCLT        ref_QuadClt,   // reference scene for ref_disparity, zero xyz and atr
+					center_disparity,      // double []      ref_disparity, // null or alternative reference disparity
+					center_CLT,            // QuadCLT        first_QuadClt, // First in comparison pair
+					pXpYD_center,          // double [][]    pXpYD_ref,     // pXpYD for the reference scene
+					reliable_ref,          // boolean []     reliable_ref, // null or bitmask of reliable reference tiles
+					tp_tasks_center[0],    // TpTask[]       tp_tasks_ref,  // only (main if MB correction) tasks for FPN correction
+					quadCLTs[scene_index], // QuadCLT        scene_QuadClt,
+					initial_pose,          // double [][]    camera_xyzatr,
+					initial_pose,          // double [][]    scene_xyzatr_pull,  // if both are not null, specify target values to pull to
+					new double[2][3],      // double [][]    ref_xyzatr, 
+					param_select,          // boolean[]      param_select,
+					reg_weights,           // double []      param_regweights,
+					lma_rms,               // double []      rms_out, // null or double [2]
+					max_rms,               // double         max_rms,
+					// motion blur compensation 
+					mb_tau,                // double         mb_tau,      // 0.008; // time constant, sec
+					mb_max_gain,           // double         mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
+					mb_vectors_scene1,     // double [][]    mb_vectors,  // now [2][ntiles];
+					debugLevel);           // int            debug_level)
+			
+			adjust_OK = scenes_xyzatr[scene_index] != null;
+    		if (adjust_OK) { // check only for initial orientation, do not check on readjustments
+    			/// No roll or zoom failure should be possible in cuas mode, keep just as assert
+				if (Math.abs(scenes_xyzatr[scene_index][1][2]) > max_roll) {
+					fail_reason[0] = FAIL_REASON_ROLL;
+					adjust_OK = false;
+				}
+				if (max_zoom_diff > 0) { // ignore if set to 
+					if (Math.abs(scenes_xyzatr[scene_index][0][2]) > max_z_change) {
+						if (Math.abs(scenes_xyzatr[scene_index][0][2]) > max_z_change_scaled) {
+							fail_reason[0] = FAIL_REASON_ZOOM;
+							adjust_OK = false;
+						} else {
+							System.out.println("Z-change "+Math.abs(scenes_xyzatr[scene_index][0][2])+
+							"m exceeds limit of "+max_z_change+"m, but as it is beginning/end of the series, limit is raised to "+
+							max_z_change_scaled+"m");
+						}
+					}
+				}
+			}
+			// update	last_corr_xyzatr, // current correction
+			double [][] corr_diff = ErsCorrection.combineXYZATR(
+					scenes_xyzatr[scene_index],                // new
+					ErsCorrection.invertXYZATR(initial_pose)); //old
+			last_corr_xyzatr = ErsCorrection.combineXYZATR (
+					last_corr_xyzatr,
+					corr_diff);
+			// not used here with IMS: after_spiral == false
+			
+			double [][] adjusted_xyzatr_dt = QuadCLTCPU.scaleDtFromErs(
+					clt_parameters,
+					quadCLTs[scene_index].getErsCorrection().getErsXYZATR_dt(
+							));
+			ers_center.addScene(scene_QuadClt.getImageName(),
+					scenes_xyzatr[scene_index][0],
+					scenes_xyzatr[scene_index][1],
+					adjusted_xyzatr_dt[0], // ZERO3, // ers_scene.getErsXYZ_dt(),		
+					adjusted_xyzatr_dt[1] // ZERO3 // ers_scene.getErsATR_dt()		
+					);
+			rmse_stats.add(lma_rms[0]);
+			if ((rmse_stats_metric != null) && (lma_rms.length >=4)) {
+				rmse_stats_metric.add(lma_rms[4]);
+			}
+			if (debugLevel > -3) {
+				System.out.println("Pass multi scene "+scene_index+" (of "+ quadCLTs.length+") "+
+						center_CLT.getImageName() + "/" + scene_QuadClt.getImageName()+
+		                " Done. Weight = "+lma_rms[2]+", number="+lma_rms[3]);
+		        System.out.print("RMS="+lma_rms[0]+
+		                ", maximal so far was "+rmse_stats.getMax()+", average was "+rmse_stats.getAverage());
+		        if ((rmse_stats_metric != null) && (lma_rms.length >=6)) {
+		            System.out.print(". Pixel RMS="+lma_rms[4]+
+		                    "pix, maximal so far was "+rmse_stats_metric.getMax()+
+		                    "pix, average was "+rmse_stats_metric.getAverage()+"pix.");
+		        }
+		        System.out.println();
+			}
+		} // for (int scene_index =  ref_index - 1; scene_index >= 0 ; scene_index--)
+		
+		if ((last_index - earliest_scene + 1) < min_num_scenes) {
+			System.out.println("Total number of useful scenes = "+(last_index - earliest_scene + 1)+
+					" < "+min_num_scenes+". Scrapping this series.");
+			if (start_ref_pointers != null) {
+				start_ref_pointers[0] = earliest_scene;
+			}
+			return -1;
+		}
+		
+		if (earliest_scene > 0) {
+			System.out.print("setInitialOrientationsCuas(): Not all scenes matched, earliest useful scene = "+earliest_scene+
+					" (total number of scenes = "+(last_index - earliest_scene + 1)+
+					"). Maximal RMSE was "+rmse_stats.getMax()+", average was "+rmse_stats.getAverage());
+	        if ((rmse_stats_metric != null) && (lma_rms.length >=4)) {
+	            System.out.print(". Maximal metric so far was "+rmse_stats_metric.getMax()+
+	                    "pix, average was "+rmse_stats_metric.getAverage()+"pix.");
+	        }
+	        System.out.println();
+			
+			
+		} else if (debugLevel > -4) {
+			System.out.print("All multi scene passes are Done. Maximal RMSE was "+rmse_stats.getMax()+
+					", average was "+rmse_stats.getAverage());
+	        if ((rmse_stats_metric != null) && (lma_rms.length >=4)) {
+	            System.out.print(". Maximal metric so far was "+rmse_stats_metric.getMax()+
+	                    "pix, average was "+rmse_stats_metric.getAverage()+"pix.");
+	        }
+	        System.out.println();
+		}
+		center_CLT.set_orient(1); // first orientation
+		center_CLT.set_accum(0);  // reset accumulations ("build_interscene") number
+		center_CLT.saveInterProperties( // save properties for interscene processing (extrinsics, ers, ...)  // null pointer
+				null, // String path,             // full name with extension or w/o path to use x3d directory
+				debugLevel+1);
+		// Add to log
+    	StringBuffer sb = new StringBuffer();
+    	sb.append(new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime())+"\n");
+		sb.append("Finished setInitialOrientationsCuas():\n");
+    	sb.append("getNumOrient()= "+center_CLT.getNumOrient()+" of "+clt_parameters.imp.min_num_orient+"\n");
+    	sb.append("getNumAccum()=  "+center_CLT.getNumAccum()+ " of "+clt_parameters.imp.min_num_interscene+"\n");
+    	sb.append("earliest_scene= "+earliest_scene+"\n");
+    	sb.append("Maximal RMSE=   "+rmse_stats.getMax()+"\n");
+    	sb.append("Average RMSE=   "+rmse_stats.getAverage()+"\n");
+    	if (rmse_stats_metric != null) {
+	    	sb.append("Maximal RMSE=   "+rmse_stats_metric.getMax()+"pix\n");
+	    	sb.append("Average RMSE=   "+rmse_stats_metric.getAverage()+"pix\n");
+    	}
+    	sb.append("------------------------\n\n");
+    	center_CLT.saveStringInModelDirectory(sb.toString(), QuadCLT.ORIENTATION_LOGS_SUFFIX); // String  suffix)
+		return earliest_scene;
+	}
+	
+	
+	
+	public static double [][] getXyzatrImsCenter(
+			final CLTParameters          clt_parameters,
+			final boolean                compensate_ims_rotation,
+			final boolean                inertial_only, 			
+			final QuadCLT                refClt,
+			final QuadCLT[]              quadCLTs, //
+			final int                    debugLevel) {
+//		boolean lock_position = true;
+		double [] quat_corr = (compensate_ims_rotation && (refClt != null))? refClt.getQuatCorr() : null; //
+//		double [] enu_corr_metric = (refClt != null) ? refClt.getENUCorrMetric() : null;
+		double [][][] ims_xyzatr;
+		if (inertial_only) {
+			ims_xyzatr = QuadCLT.integratePIMU(
+						clt_parameters, // final CLTParameters clt_parameters,
+						quadCLTs,       // final QuadCLT[]     quadCLTs,
+						quadCLTs.length - 1, // ref_index,      // final int           ref_index,
+						null,           // double [][][]       dxyzatr,
+			    		0, // final int           early_index,
+			    		(quadCLTs.length -1), // int           last_index,
+//			    		(new double[2][3]) //final double [][] pimu_offsets // null - will not subtract from velocities
+						debugLevel      // final int           debugLevel
+					);
+		} else {
+			ims_xyzatr = refClt.getXyzatrIms(
+					clt_parameters, // CLTParameters clt_parameters,
+					quadCLTs,       // QuadCLT[]  quadCLTs,
+					quat_corr,      //  here! quat_corr, // double []  quat_corr, // only applies to rotations - verify!
+					debugLevel) ;   // int        debugLevel)
+		}
+		double [][] ims_atr = new double [ims_xyzatr.length][];
+		for (int i = 0; i < quadCLTs.length; i++) if (ims_xyzatr[i] != null) {
+			ims_atr[i] = ims_xyzatr[i][1];
+		}
+		int [] range = {0, quadCLTs.length-1};
+		double [][] center_atr =  CuasCenterLma.getCenterATR(
+				ims_atr,      // double [][] ims_atr,
+			     range,       // int []     range,
+			     debugLevel); // int        debugLevel)
+		
+		return center_atr;
+	}
+	
+
 	public static void manualIMSCorrection (
 			double [][][] xyz_atr,
 			int           corr_indx, // 57 correct scenes LE this
@@ -1470,9 +1895,10 @@ public class Interscene {
 				ref_index,      // final int           ref_index,
 				null,           // double [][][]       dxyzatr,
 	    		0, // final int           early_index,
-	    		(quadCLTs.length -1) // int           last_index,
+	    		(quadCLTs.length -1), // int           last_index,
 //	    		(new double[2][3]) //final double [][]   offsets // null - will not subtract from velocities
-	    		);
+				debugLevel      // final int           debugLevel
+				);
     	double [][][] dxyzatr =     QuadCLT.getDxyzatrPIMU(
         		clt_parameters,
         		quadCLTs); //,
@@ -1621,12 +2047,15 @@ public class Interscene {
 		int []    fail_reason = new int[1];   // null or int[1]: 0 - OK, 2 - LMA, 3 - min, 4 - max
 		for (int scene_index =  ref_index - 1; scene_index >= 0 ; scene_index--) {
 			// to include ref scene photometric calibration
-			quadCLTs[scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT(
-					set_channels[scene_index].set_name,
-					clt_parameters,
-					colorProcParameters, //
-					threadsMax,
-					debugLevel-2);
+			// added 05.21.2024: skip if was spawn already
+			if (quadCLTs[scene_index] == null) {
+				quadCLTs[scene_index] = quadCLTs[ref_index].spawnNoModelQuadCLT(
+						set_channels[scene_index].set_name,
+						clt_parameters,
+						colorProcParameters, //
+						threadsMax,
+						debugLevel-2);
+			}
 		} // split cycles to remove output clutter
 		ErsCorrection ers_reference = quadCLTs[ref_index].getErsCorrection();		
 		int debug_scene = -15;
@@ -2149,6 +2578,8 @@ public class Interscene {
 		boolean copy_pull_current = clt_parameters.imp.copy_pull_current; //false; // true;
 		boolean restore_imu =       clt_parameters.imp.restore_imu;       //false; // restore imu omega-az and omega-tl, freeze ERS, adjust X,Y,Z,Az,Tl,Rl
 		boolean lock_position =      clt_parameters.imp.lock_position;
+		double  ref_smooth_diff =    clt_parameters.imp.ref_smooth_diff; // 0.75;  // minimal fraction of the SfM maximal gain
+//		boolean ref_smooth =         clt_parameters.imp.ref_smooth;  //true;   // use SfM filtering if available 
 		boolean[] param_select =  configured_lma? clt_parameters.ilp.ilma_lma_select :
 			ErsCorrection.getParamSelect(
 					lma_xyzatr || use_Z,   // boolean use_Z,
@@ -2793,7 +3224,7 @@ public class Interscene {
 						fpn_pairs[ipair][0],                  // int            nscene1,       // compares to nscene0
 						interscene_ref_disparity,// double []      ref_disparity, // null or alternative reference disparity
 						reliable_ref,            // boolean []     reliable_ref, // null or bitmask of reliable reference tiles
-    					false,               // boolean        smooth_disparity, // smooth disparity (according to clt_parameters)+update reliable_ref if true
+						false,               // boolean        smooth_disparity, // smooth disparity (according to clt_parameters)+update reliable_ref if true
 						scenes_xyzatr[fpn_pairs[ipair][1]],// double [][]    scene0_xyzatr,
 						scenes_xyzatr[fpn_pairs[ipair][0]],            // double [][]    scene1_xyzatr,
 						avg_z,                   //  double         average_z,     
@@ -2870,7 +3301,7 @@ public class Interscene {
 						test_pair[0],                  // int            nscene1,       // compares to nscene0
 						interscene_ref_disparity,// double []      ref_disparity, // null or alternative reference disparity
 						reliable_ref,            // boolean []     reliable_ref, // null or bitmask of reliable reference tiles
-    					false,               // boolean        smooth_disparity, // smooth disparity (according to clt_parameters)+update reliable_ref if true
+						false,               // boolean        smooth_disparity, // smooth disparity (according to clt_parameters)+update reliable_ref if true
 						scenes_xyzatr[test_pair[1]],// double [][]    scene0_xyzatr,
 						scenes_xyzatr[test_pair[0]],            // double [][]    scene1_xyzatr,
 						avg_z,                   //  double         average_z,     
@@ -3016,11 +3447,14 @@ public class Interscene {
 		double  ref_sigma =          clt_parameters.imp.ref_sigma;  //10.0;   // minimal value of the SfM gain maximum to consider available
 		double  ref_smooth_diff =    clt_parameters.imp.ref_smooth_diff; // 0.75;  // minimal fraction of the SfM maximal gain
 		boolean apply_nan = true;
+		if (clt_parameters.imp.ref_smooth_always) {
+			smooth_disparity=true; // **********************************************
+		}
+		
 //	    boolean use_lma_dsi =          clt_parameters.imp.use_lma_dsi;
 		if (ref_disparity == null) { 
 			ref_disparity = quadCLTs[ref_index].getDLS()[use_lma_dsi?1:0];
 		}
-		smooth_disparity=true;
 		if (smooth_disparity) {
 			double [] dbg_disparity =  (debugLevel > 2) ? ref_disparity.clone() : null;
 			boolean [] dbg_reliable = ((dbg_disparity != null) && (reliable_ref != null)) ?  reliable_ref.clone(): null;
@@ -3108,6 +3542,7 @@ public class Interscene {
 					quadCLTs[nscene0],  // QuadCLT            ref_scene,
 					ref_disparity,// null      // double []          ref_disparity, // null or alternative reference disparity
 					pXpYD_scene0,       // double [][]        ref_pXpYD,
+					null, 			            // final float [][]          fclt, 
 					reliable_ref,       // final boolean []   selection, // may be null, if not null do not  process unselected tiles
 					margin,             // final int          margin,
 					// motion blur compensation 
@@ -3134,7 +3569,7 @@ public class Interscene {
 					0,                     // int            shrink_gaps,  // will gaps, but not more that grow by this
 					debugLevel);           // int            debug_level)
 		}
-// avg_z		
+// avg_z	
 		double pull_offset = 0;
 		if (!Double.isNaN(average_z) && (scene1_xyzatr_pull != null) && (scene1_xyzatr_pull[0] != null) && (scene1_xyzatr_pull[1] != null)){
 			pull_offset = quadCLTs[ref_index].estimateAverageShift(
@@ -3321,6 +3756,7 @@ public class Interscene {
 				reference_QuadClt,  // QuadCLT            ref_scene,
 				ref_disparity,      // double []          ref_disparity, // null or alternative reference disparity
 				pXpYD_ref,          // double [][]        ref_pXpYD,
+				null, 			    // final float [][]          fclt, 
 				reliable_ref,       // final boolean []   selection, // may be null, if not null do not  process unselected tiles
 				margin,             // final int          margin,
 				// motion blur compensation 
@@ -3780,6 +4216,7 @@ public class Interscene {
     			reference_QuadClt,  // QuadCLT            ref_scene,
     			ref_disparity,      // double []          ref_disparity, // null or alternative reference disparity
     			pXpYD_ref,          // double [][]        ref_pXpYD,
+				null, 			    // final float [][]          fclt, 
     			reliable_ref,       // final boolean []   selection, // may be null, if not null do not  process unselected tiles
     			margin,             // final int          margin,
     			// motion blur compensation 
@@ -4288,7 +4725,7 @@ public class Interscene {
 					debug_level);               // final int                 globalDebugLevel);
 		}
 		if (show_tptask_scene) {
-			for (int i = 0; i < tp_tasks.length; i++) if (tp_tasks_ref[i] != null){
+			for (int i = 0; i < tp_tasks.length; i++) if (tp_tasks[i] != null){
 				String dbg_title = scene.getImageName()+"-scn-"+i;
 				TpTask.showTpTask(
 						tp_tasks[i], // TpTask[] tp_tasks,
@@ -5606,6 +6043,7 @@ public class Interscene {
 	 * @param ref_scene     reference scene QuadCLT instance
 	 * @param ref_disparity either alternative disparity array or null to use it from the reference scene itself
 	 * @param ref_pXpYD - precalculated or null (will be calculated)
+	 * @param fclt load TD data and bypass execConvertDirect. If fclt.length==1 load same data to all sensor channels
 	 * @param selection optional selection to ignore unselected tiles)
 	 * @param margin do not use tiles with centers closer than this to the edges. Measured in pixels.
 	 * @param mb_tau Sensor time constant in seconds (only needed if mb_vectors != null)
@@ -5619,6 +6057,7 @@ public class Interscene {
 			QuadCLT            ref_scene,
 			double []          ref_disparity, // null or alternative reference disparity
 			double [][]        ref_pXpYD,
+			final float [][]   fclt, 
 			final boolean []   selection, // may be null, if not null do not  process unselected tiles
 			final int          margin,
 			// motion blur compensation 
@@ -5707,6 +6146,7 @@ public class Interscene {
 					QuadCLT.THREADS_MAX);                  // final int                 threadsMax)  // maximal number of threads to launch
 			ref_scene.saveQuadClt(); // to re-load new set of Bayer images to the GPU (do nothing for CPU) and Geometry
 			image_dtt.setReferenceTD( // tp_tasks_ref will be updated
+					fclt, 			            // final float [][]          fclt, 
 					erase_clt,
 					null,                       // final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 					clt_parameters.img_dtt,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
@@ -5845,8 +6285,9 @@ public class Interscene {
 				ref_index,      // final int           ref_index,
 				null,           // double [][][]       dxyzatr,
 				earliest_scene, // final int           early_index,
-	    		(quadCLTs.length -1) // int           last_index,
-	    		);
+	    		(quadCLTs.length -1), // int           last_index,
+				debugLevel      // final int           debugLevel
+				);
 		int cent_index =        earliest_scene + (quadCLTs.length - earliest_scene) / 2;
 		double [][] cent_xyzatr = pimu_xyzatr[cent_index];
 //		int            ref_index,

src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java

@@ -253,7 +253,8 @@ public class IntersceneMatchParameters {
 
 	
 	public  boolean center_reference =   false;
-	public  boolean lock_position =      false; // camera position is not changing, only oriantation (cuas mode)
+	public  boolean lock_position =      false; // camera position is not changing, only orientation (cuas mode)
+	public  boolean cuas_rotation =      false; // assumes lock_position, camera is rotating around the center
 	public  boolean manual_correction =  false; // once used for foliage to merge two sequences
 	public  boolean overlap_sequences =  false; // overlap sequences: scan down from the previous center  
 	public  boolean reset_photometric =  true;  // reset photometric calibration - once for each new series
@@ -564,7 +565,10 @@ min_str_neib_fpn	0.35
 	public  double  min_ref_str_lma =               0.8;   // 0.22;  // For orientations: use only tiles of the reference scene DSI_MAIN is stronger  
 	public  double  min_ref_frac =                  0.2;   // 0.22;  if fraction number of reliable tiles is less than this, use best possible
 	
-	public  boolean ref_smooth =                    false; // smooth reference disparity for initial matching
+	// was overwritten in code to always be true ****
+	public  boolean ref_smooth =                    true; // false; // smooth reference disparity for initial matching
+	public  boolean ref_smooth_always =             true; // in older code it was hard-wired true inside adjustDiffPairsLMAInterscene
+//	ref_smooth_always, when true mimics old behavior
 	public  double  ref_sigma =                     2.0;   // Gaussian sigma to smooth reference disparity for initial matching
 	public  double  ref_smooth_diff =               0.15;  // discard smooth disparity if it differs by more from measured
 	
@@ -1247,7 +1251,9 @@ min_str_neib_fpn	0.35
 				"True for mapping with post-processing. If false, use the last scene as a reference (driving).");
 		gd.addCheckbox ("Lock camera position",                      this.lock_position,
 				"Only camera orientation changes, position is locked (cUAS mode).");
-		gd.addCheckbox ("Manual correction",                      this.manual_correction,
+		gd.addCheckbox ("CUAS mode",                                 this.cuas_rotation,
+				"Assumes lock_position, camera is rotating around the center.");
+		gd.addCheckbox ("Manual correction",                         this.manual_correction,
 				"Once used for foliage to merge two sequences.");
 		gd.addCheckbox ("Overlap sequences",                         this.overlap_sequences,
 				"Overlap sequences by half - start (down) from the previous center.");
@@ -1786,6 +1792,8 @@ min_str_neib_fpn	0.35
 		gd.addMessage  ("Smooth reference disparity for initial matching");
 		gd.addCheckbox ("Smooth reference disparity",                this.ref_smooth,
 				"Smooth reference disparity for initial matching.");
+		gd.addCheckbox ("Smooth always",                             this.ref_smooth_always,
+				"In older code it was hard-wired true inside adjustDiffPairsLMAInterscene. Checked mimics that legacy behavior.");
 		gd.addNumericField("Reference disparity sigma",              this.ref_sigma, 5,7,"pix",
 				"Gaussian sigma to smooth reference disparity for initial matching.");
 		gd.addNumericField("Max smooth offset",                      this.ref_smooth_diff, 5,7,"pix",
@@ -2501,6 +2509,8 @@ min_str_neib_fpn	0.35
         
 		this.center_reference =               gd.getNextBoolean();
 		this.lock_position =                  gd.getNextBoolean();
+		this.cuas_rotation =                  gd.getNextBoolean();
+		this.lock_position |= this.cuas_rotation; 
 		this.manual_correction =              gd.getNextBoolean();
 		this.overlap_sequences =              gd.getNextBoolean();
 		this.reset_photometric =              gd.getNextBoolean();
@@ -2758,6 +2768,7 @@ min_str_neib_fpn	0.35
 		this.min_ref_str_lma =          gd.getNextNumber();
 		this.min_ref_frac =             gd.getNextNumber();
 		this.ref_smooth =               gd.getNextBoolean();
+		ref_smooth_always =             gd.getNextBoolean();		
 		this.ref_sigma =                gd.getNextNumber();
 		this.ref_smooth_diff =          gd.getNextNumber();
 		this.sfm_filter =               gd.getNextBoolean();
@@ -3301,6 +3312,7 @@ min_str_neib_fpn	0.35
 
 		properties.setProperty(prefix+"center_reference",              this.center_reference + "");            // boolean
 		properties.setProperty(prefix+"lock_position",                 this.lock_position + "");               // boolean
+		properties.setProperty(prefix+"cuas_rotation",                 this.cuas_rotation + "");               // boolean
 		properties.setProperty(prefix+"manual_correction",             this.manual_correction + "");           // boolean
 		properties.setProperty(prefix+"overlap_sequences",             this.overlap_sequences + "");           // boolean
 		properties.setProperty(prefix+"reset_photometric",             this.reset_photometric + "");           // boolean
@@ -3569,6 +3581,7 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"min_ref_frac",         this.min_ref_frac+"");        // double
 		
 		properties.setProperty(prefix+"ref_smooth",           this.ref_smooth+"");          // boolean
+		properties.setProperty(prefix+"ref_smooth_always",    this.ref_smooth_always+"");   // boolean
 		properties.setProperty(prefix+"ref_sigma",            this.ref_sigma+"");           // double
 		properties.setProperty(prefix+"ref_smooth_diff",      this.ref_smooth_diff+"");     // double
 		
@@ -4083,6 +4096,8 @@ min_str_neib_fpn	0.35
 		
 		if (properties.getProperty(prefix+"center_reference")!=null)              this.center_reference=Boolean.parseBoolean(properties.getProperty(prefix+"center_reference"));		
 		if (properties.getProperty(prefix+"lock_position")!=null)                 this.lock_position=Boolean.parseBoolean(properties.getProperty(prefix+"lock_position"));		
+		if (properties.getProperty(prefix+"cuas_rotation")!=null)                 this.cuas_rotation=Boolean.parseBoolean(properties.getProperty(prefix+"cuas_rotation"));
+		this.lock_position |= this.cuas_rotation; 
 		if (properties.getProperty(prefix+"manual_correction")!=null)             this.manual_correction=Boolean.parseBoolean(properties.getProperty(prefix+"manual_correction"));		
 		if (properties.getProperty(prefix+"overlap_sequences")!=null)             this.overlap_sequences=Boolean.parseBoolean(properties.getProperty(prefix+"overlap_sequences"));		
 		if (properties.getProperty(prefix+"reset_photometric")!=null)             this.reset_photometric=Boolean.parseBoolean(properties.getProperty(prefix+"reset_photometric"));		
@@ -4356,6 +4371,7 @@ min_str_neib_fpn	0.35
 		if (properties.getProperty(prefix+"min_ref_frac")!=null)         this.min_ref_frac=Double.parseDouble(properties.getProperty(prefix+"min_ref_frac"));
 
 		if (properties.getProperty(prefix+"ref_smooth")!=null)           this.ref_smooth=Boolean.parseBoolean(properties.getProperty(prefix+"ref_smooth"));		
+		if (properties.getProperty(prefix+"ref_smooth_always")!=null)    this.ref_smooth_always=Boolean.parseBoolean(properties.getProperty(prefix+"ref_smooth_always"));		
 		if (properties.getProperty(prefix+"ref_sigma")!=null)            this.ref_sigma=Double.parseDouble(properties.getProperty(prefix+"ref_sigma"));
 		if (properties.getProperty(prefix+"ref_smooth_diff")!=null)      this.ref_smooth_diff=Double.parseDouble(properties.getProperty(prefix+"ref_smooth_diff"));
 
@@ -4882,7 +4898,8 @@ min_str_neib_fpn	0.35
         imp.patt_save_subdir              = this.patt_save_subdir;
         
 		imp.center_reference              = this.center_reference;
-		imp.lock_position                 = this.lock_position;
+		imp.lock_position                 = this.lock_position || this.cuas_rotation;
+		imp.cuas_rotation                 = this.cuas_rotation;
 		imp.manual_correction             = this.manual_correction;
 		imp.overlap_sequences             = this.overlap_sequences;
 		imp.reset_photometric             = this.reset_photometric;
@@ -5146,6 +5163,7 @@ min_str_neib_fpn	0.35
 		imp.min_ref_frac          = this.min_ref_frac;
 		
 		imp.ref_smooth =            this.ref_smooth;
+		imp.ref_smooth_always =     this.ref_smooth_always;
 		imp.ref_sigma =             this.ref_sigma;
 		imp.ref_smooth_diff =       this.ref_smooth_diff;
 		

src/main/java/com/elphel/imagej/tileprocessor/OpticalFlow.java

@@ -54,6 +54,7 @@ import com.elphel.imagej.common.DoubleGaussianBlur;
 import com.elphel.imagej.common.PolynomialApproximation;
 import com.elphel.imagej.common.ShowDoubleFloatArrays;
 import com.elphel.imagej.correction.CorrectionColorProc;
+import com.elphel.imagej.cuas.Cuas;
 import com.elphel.imagej.cuas.CuasCenterLma;
 import com.elphel.imagej.gpu.GPUTileProcessor;
 import com.elphel.imagej.gpu.GpuQuad;
@@ -97,6 +98,7 @@ public class OpticalFlow {
 	public static int COMBO_DSN_INDX_GROUND =     13; // Smooth ground surface w/o vegetation
 	public static int COMBO_DSN_INDX_TERRAIN =    14; // Smooth ground surface w/o vegetation
 	
+	public static int [] COMBO_DSN_NONNAN = {COMBO_DSN_INDX_STRENGTH, COMBO_DSN_INDX_STRENGTH_BG,COMBO_DSN_INDX_SFM_GAIN};
 // move to Interscene class?	
 	// interscene adjustments failure reasons.
 	public static final int FAIL_REASON_LMA =            1; // LMA failed
@@ -122,7 +124,6 @@ public class OpticalFlow {
 			return "unknown scene-matching failure reason="+fr;
 		}
 	}
-	
 	public static double [] ZERO3 = {0.0,0.0,0.0};
 	public static double  LINE_ERR = 0.1;
 	public static int     THREADS_MAX =          100;  // maximal number of threads to launch
@@ -3025,232 +3026,7 @@ public class OpticalFlow {
 		return pXpYD;
 	}
 	
-	
-	
-	public static double [][] transformFromVirtual(
-			double []         disparity_ref,
-			final double []   scene_xyz,         // camera center in world (reference) coordinates
-			final double []   scene_atr,         // camera orientation relative to world (reference) frame
-			final QuadCLT     reference_QuadClt,
-			final int         around,            // 2 search around for interpolation
-			final double      sigma,
-			final int         num_refines,
-			final String      debugSuffix){
-		final boolean debug = debugSuffix != null;
-		final double sigma2 = 2*sigma*sigma;
-    	final double normal_damping = 0.001; // pull to horizontal if not enough data 
-    	final double [] damping = new double [] {normal_damping, normal_damping};
-		TileProcessor tp =         reference_QuadClt.getTileProcessor();
-		final int tilesX =         tp.getTilesX();
-		final int tilesY =         tp.getTilesY();
-		final int tiles =          tilesY*tilesX;
-		final int tilesX_around =  tilesX + 2 * around; 
-		final int tilesY_around =  tilesY + 2 * around; 
-		final int tiles_around =   tilesY_around * tilesX_around;
-		final int transform_size = tp.getTileSize();
-		final Thread[] threads = ImageDtt.newThreadArray(ImageDtt.THREADS_MAX);
-		final AtomicInteger ai = new AtomicInteger(0);
-		int initial_capacity = 4;
-		final int min_samples = 1;
-		
-		final String debugTitle=debug?reference_QuadClt.getImageName()+"-"+debugSuffix : null;
-		String [] debug_frame_titles = {"X", "Y", "D"};
-		String [] debug_titles = new String[num_refines+1];
-		for (int i = 0; i <= num_refines; i++) {
-			debug_titles[i]=i+"";
-		}
-		final double [][][] debug_data = debug ? new double [3][num_refines+1][tiles]:null;
-		if (debug) {
-			for (int f = 0; f < debug_data.length; f++) {
-				for (int i = 0; i < debug_data[f].length; i++) {
-					Arrays.fill(debug_data[f][i], Double.NaN);
-				}
-			}
-		}
-
-		final ArrayList<ArrayList<Integer>> interp_list = new ArrayList<ArrayList<Integer>>(tiles_around);
-		for (int nTile = 0; nTile < tiles_around; nTile++) {
-			interp_list.add(new ArrayList<Integer>(initial_capacity));
-		}
-		// create uniform grid for initial interpolations
-		final double [][] reference_pXpYD = new double [tiles][];
-		for (int nrefine = 0; nrefine <= num_refines; nrefine++) {
-			if (nrefine == 0) {
-				ai.set(0);
-				for (int ithread = 0; ithread < threads.length; ithread++) {
-					threads[ithread] = new Thread() {
-						public void run() {
-							for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) if (!Double.isNaN(disparity_ref[nTile])){
-								double disparity = disparity_ref[nTile];
-								int tileY = nTile / tilesX;  
-								int tileX = nTile % tilesX;
-								double centerX = (tileX + 0.5) * transform_size; //  - shiftX;
-								double centerY = (tileY + 0.5) * transform_size; //  - shiftY;
-								if (disparity < 0) {
-									disparity = 1.0* disparity; // 0.0;
-								}
-								reference_pXpYD[nTile] = new double[] {centerX, centerY, disparity};
-							}
-						}
-					};
-				}		      
-				ImageDtt.startAndJoin(threads);
-				
-			} else {
-				// get scene pXpYD corresponding to reference_pXpYD
-				final double [][] scene_pXpYD = transformToScenePxPyD(
-						reference_pXpYD,    // final double [][] reference_pXpYD,// invalid tiles - NaN in disparity. Should be no nulls, no NaN disparity
-						scene_xyz,          // final double []   scene_xyz,         // camera center in world (reference) coordinates
-						scene_atr,          // final double []   scene_atr,         // camera orientation relative to world (reference) frame
-						reference_QuadClt, // final QuadCLT     reference_QuadClt) //
-						null); // final QuadCLT     scene_QuadClt) //
-
-				// deep clone reference_pXpYD
-				final double [][] reference_pXpYD_next = new double [tiles][];
-				ai.set(0);
-				for (int ithread = 0; ithread < threads.length; ithread++) {
-					threads[ithread] = new Thread() {
-						public void run() {
-							for (int nTile = ai.getAndIncrement(); nTile < tiles_around; nTile = ai.getAndIncrement()) {
-								interp_list.get(nTile).clear(); // clear each list
-							}
-						}
-					};
-				}		      
-				ImageDtt.startAndJoin(threads);
-				// not multithreaded 
-				for (int nTile = 0; nTile < tiles; nTile++) if (scene_pXpYD[nTile] != null) {
-					int tileX = (int) Math.floor (scene_pXpYD[nTile][0]/transform_size);
-					int tileY = (int) Math.floor (scene_pXpYD[nTile][1]/transform_size);
-					int tileX_around = tileX + around; 
-					int tileY_around = tileY + around;
-					if ((tileX_around >= 0) && (tileY_around >= 0) && (tileX_around < tilesX_around) && (tileY_around < tilesY_around)) {
-						int tile_around = tileX_around + tilesX_around * tileY_around;
-						interp_list.get(tile_around).add(nTile);
-					}
-				}
-				// interpolate reference_pXpYD_next 
-				ai.set(0);
-				for (int ithread = 0; ithread < threads.length; ithread++) {
-					threads[ithread] = new Thread() {
-						public void run() {
-							PolynomialApproximation pa = new PolynomialApproximation();		
-							for (int nTile = ai.getAndIncrement(); nTile < tiles; nTile = ai.getAndIncrement()) {
-								int tileY = nTile / tilesX;  
-								int tileX = nTile % tilesX;
-								double centerX = (tileX + 0.5) * transform_size; //  - shiftX;
-								double centerY = (tileY + 0.5) * transform_size; //  - shiftY;
-								int tileX_around = tileX + around; 
-								int tileY_around = tileY + around;
-								int num_samples = 0;
-								for (int dy = -around; dy <= around; dy++) {
-									for (int dx = -around; dx <= around; dx++) {
-										int tile_around = (tileX_around + dx) + tilesX_around * (tileY_around + dy);
-										num_samples+=interp_list.get(tile_around).size();
-									}
-								}
-								int mindx = 0;
-								if (num_samples >= min_samples) {
-									// inter/extrapolate with regularization (for few samples) and weights
-									double [][][] mdata = new double [num_samples][3][];
-									for (int dy = -around; dy <= around; dy++) {
-										for (int dx = -around; dx <= around; dx++) {
-											int tile_around = (tileX_around + dx) + tilesX_around * (tileY_around + dy);
-											for (int nst:interp_list.get(tile_around)) {
-												double [] ref_xyd = reference_pXpYD[nst];
-												double [] scene_xyd = scene_pXpYD[nst];
-												double x = scene_xyd[0]-centerX;
-												double y = scene_xyd[1]-centerY;
-												double w = Math.exp (-(x*x + y*y)/sigma2);
-												mdata[mindx][0] = new double [2];
-												mdata[mindx][0][0] =  x;
-												mdata[mindx][0][1] =  y;
-												mdata[mindx][1] = new double [2]; // [3];
-												mdata[mindx][1][0] =  ref_xyd[0]; 
-												mdata[mindx][1][1] =  ref_xyd[1]; 
-//												mdata[mindx][1][2] =  ref_xyd[2];
-												mdata[mindx][2] = new double [1];
-												mdata[mindx][2][0] =  w;
-												mindx++;
-											}
-										}
-									}
-									double[][] approx2d = pa.quadraticApproximation(
-											mdata,
-											true,       // boolean forceLinear,  // use linear approximation
-											damping,    // double [] damping, null OK
-											-1);        // debug level
-									double px = approx2d[0][2];
-									double py = approx2d[1][2];
-									// interpolate disparity from disparity_ref
-									double tx = px/transform_size-0.5;
-									double ty = py/transform_size-0.5;
-									int tx0 = (int) Math.floor(tx);
-									int ty0 = (int) Math.floor(ty);
-									int tx1 = tx0+1;
-									int ty1 = ty0+1;
-									double kx = tx-tx0;
-									double ky = ty-ty0;
-									if (tx0 < 0) {
-										tx0 = 0;
-										tx1 = tx0+1;
-										kx = 0;
-									} else if (tx1 >= tilesX) {
-										tx1 = tilesX - 1;
-										tx0 = tx1-1;
-										kx = 1.0;
-									}
-									
-									if (ty0 < 0) {
-										ty0 = 0;
-										ty1 = ty0+1;
-										ky = 0;
-									} else if (ty1 >= tilesY) {
-										ty1 = tilesY - 1;
-										ty0 = ty1-1;
-										ky = 1.0;
-									}
-									double d = 
-											(1-kx)*(1-ky)*disparity_ref[tx0 + tilesX*ty0] +
-											(  kx)*(1-ky)*disparity_ref[tx1 + tilesX*ty0] +
-											(1-kx)*(  ky)*disparity_ref[tx0 + tilesX*ty1] +
-											(  kx)*(  ky)*disparity_ref[tx1 + tilesX*ty1];
-									reference_pXpYD_next[nTile] = new double[] {px, py, d};
-								}
-							}
-						}
-					};
-				}		      
-				ImageDtt.startAndJoin(threads);
-				System.arraycopy(
-						reference_pXpYD_next,
-						0,
-						reference_pXpYD,
-						0,
-						tiles);
-				// TODO: refine disparity from disparity_ref interpolation
-			} // if (nrefine == 0) else
-			if (debug) {
-				for (int nTile = 0; nTile < tiles; nTile++) if (reference_pXpYD[nTile] != null) {
-					debug_data[0][nrefine][nTile] = reference_pXpYD[nTile][0];
-					debug_data[1][nrefine][nTile] = reference_pXpYD[nTile][1];
-					debug_data[2][nrefine][nTile] = reference_pXpYD[nTile][2];
-				}
-			}
-		}
-		if (debug) {
-			ShowDoubleFloatArrays.showArraysHyperstack(
-					debug_data,         // double[][][] pixels, 
-					tilesX,             // int          width, 
-					debugTitle,         // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
-					debug_titles,       // String []    titles, // all slices*frames titles or just slice titles or null
-					debug_frame_titles, // String []    frame_titles, // frame titles or null
-					true);              // boolean      show)
-			
-		}		
-		
-		return reference_pXpYD; // ref_pXpYD_0; // 
-	}
+
 	
 	//TODO: refine inter-scene pose to accommodate refined disparity map
 	/**
@@ -5161,6 +4937,8 @@ public class OpticalFlow {
 				}
 			}
 		}
+		
+		
 		if (!build_ref_dsi && (quadCLTs[last_index] != null)) {
 			quadCLTs[last_index].restoreInterProperties(null, false, debugLevel); //null
 		}
@@ -5181,10 +4959,11 @@ public class OpticalFlow {
 					colorProcParameters, //
 					threadsMax,
 					debugLevel-2);
-			
-			
 		}
 		
+		//************* move loading center DSI here before building reference DSI. Or after "blue sky 
+
+		
 		while ((quadCLTs[last_index] == null) || !quadCLTs[last_index].hasBlueSky()) { // null
 			if (build_ref_dsi) {
 				TwoQuadCLT.copyJP4src( // actually there is no sense to process multiple image sets. Combine with other
@@ -5310,12 +5089,91 @@ public class OpticalFlow {
 
 		// Build initial orientations
 		boolean ims_use =      center_reference || clt_parameters.imp.ims_use; // center works with IMS only
-///		int center_index=ref_index;
 		int ref_index = last_index; // old versions
+		double [][] xyzatr_ims_center = null;
+		boolean use_cuas = clt_parameters.imp.cuas_rotation & ims_use; // lock_position & ims_use; // needs
+		QuadCLT center_CLT = null;
+		double [] atr_ims_center = null;
 		if (force_initial_orientations && !reuse_video) {
+			if (use_cuas) {
+				// try to restore if exists
+				center_CLT = QuadCLT.restoreCenterClt(
+						quadCLTs[last_index]); // ref_index]);
+				if (center_CLT != null) { // no sense to calculate center if CLT data is not available (for now will need to manually copy)
+					// to this vXX version
+					if (debugLevel > -3) {
+						System.out.println("Restored center CLT in CUAS mode");
+					}
+					// quadCLTs[last_index] should be known here
+//					for (int scene_index =  ref_index - 1; scene_index >= 0 ; scene_index--) if (quadCLTs[scene_index] == null){
+					for (int scene_index =  last_index - 1; scene_index >= 0 ; scene_index--) if (quadCLTs[scene_index] == null){
+						// to include ref scene photometric calibration
+						quadCLTs[scene_index] = quadCLTs[last_index].spawnNoModelQuadCLT(
+								set_channels[scene_index].set_name,
+								clt_parameters,
+								colorProcParameters, //
+								threadsMax,
+								debugLevel-2);
+					} // split cycles to remove output clutter
+					xyzatr_ims_center = Interscene.getXyzatrImsCenter( // use [0], [1] - radii
+							clt_parameters,          // final CLTParameters          clt_parameters,
+							use_ims_rotation,        // final boolean                compensate_ims_rotation,
+							inertial_only,           // final boolean                inertial_only, 			
+							quadCLTs[last_index],    // final QuadCLT                refClt,
+							quadCLTs,                //final QuadCLT[]              quadCLTs, //
+							debugLevel);             // final int                    debugLevel)
+					atr_ims_center = new double[] {xyzatr_ims_center[0][0],xyzatr_ims_center[0][1],xyzatr_ims_center[0][2]};
+					if (debugLevel > -3) {
+						System.out.println("xyzatr_ims_center[0]: A="+xyzatr_ims_center[0][0]+", T="+xyzatr_ims_center[0][1]+", R="+xyzatr_ims_center[0][2]);
+					}
+				} else {
+					if (debugLevel > -3) {
+						System.out.println("Center CLT for CUAS mode does not exist");
+					}
+				}
+				
+			}
+			
 			boolean OK = false;
 			int es1 = -1;
-			if (center_reference) {
+			/*
+	public static int setInitialOrientationsCuas(
+			final CLTParameters          clt_parameters,
+			final double [][]            center_xyzatr, // center_xyzatr[0]=ZERO3 - for this sequence
+			final QuadCLT                center_CLT,    // contains center CLT and DSI
+			final boolean                compensate_ims_rotation,
+			final boolean                inertial_only, 			
+			int                          min_num_scenes,
+			final ColorProcParameters    colorProcParameters,
+			final QuadCLT[]              quadCLTs, //
+			final QuadCLT.SetChannels [] set_channels,
+			final boolean                batch_mode,
+			int                          earliest_scene,
+			int []                       start_ref_pointers, // [0] - earliest valid scene, [1] ref_index. Set if non-null
+			final int                    threadsMax,  // int               threadsMax,
+			final boolean                updateStatus,
+			final int                    debugLevel) {
+
+			 */
+			if (center_CLT != null) {
+				double [][] center_xyzatr = {ZERO3,atr_ims_center};
+				// int [] start_ref_pointers = {earliest_scene, last_index};
+				es1 = Interscene.setInitialOrientationsCuas(
+						clt_parameters,      // final CLTParameters          clt_parameters,
+						center_xyzatr, // final double [][]            center_xyzatr, // center_xyzatr[0]=ZERO3 - for this sequence
+						center_CLT, // final QuadCLT                center_CLT,    // contains center CLT and DSI
+						use_ims_rotation,    // final boolean                compensate_ims_rotation,
+						inertial_only,       // final boolean                                inertial_only,
+						min_num_scenes,      // int                          min_num_scenes,
+						colorProcParameters, // final ColorProcParameters    colorProcParameters,
+						quadCLTs,            // final QuadCLT[]              quadCLTs, //
+						set_channels,        // final QuadCLT.SetChannels [] set_channels,
+						batch_mode,          // final boolean                batch_mode,
+						start_ref_pointers,  // int []                       start_ref_pointers, // [0] - earliest valid scene, [1] ref_index. Set if non-null
+						updateStatus,        // final boolean                updateStatus,
+						debugLevel);
+				OK = earliest_scene >=0;
+			} else if (center_reference) {
 				es1 = Interscene.setInitialOrientationsCenterIms(
 						clt_parameters,      // final CLTParameters          clt_parameters,
 						use_ims_rotation,    // final boolean                compensate_ims_rotation,
@@ -5902,7 +5760,7 @@ public class OpticalFlow {
 				
 			}
 		}
-		boolean extract_center_orientation = clt_parameters.imp.extract_center_orientation; // true; // false; // true; 
+		boolean extract_center_orientation = false; //  clt_parameters.imp.extract_center_orientation; // true; // false; // true; 
 		double [][] center_ATR = null; //  {{center_A, center_T, average_R},{radius_A, radius_T}}
 		double [] cuas_atr = ZERO3;
 		if (extract_center_orientation && clt_parameters.imp.lock_position) {
@@ -5915,12 +5773,79 @@ public class OpticalFlow {
 			cuas_atr = new double [] { center_ATR[0][0], center_ATR[0][1], center_ATR[0][2]};
 		}
 		
-		boolean combine_clt = (cuas_atr[0] != 0) || (cuas_atr[1] != 0) || (cuas_atr[2] != 0);
+		boolean combine_clt = clt_parameters.imp.cuas_rotation; // (cuas_atr[0] != 0) || (cuas_atr[1] != 0) || (cuas_atr[2] != 0);
 		if (combine_clt) {
-			boolean apply_clt = true; // set GPU with data
-			boolean show_clt =  true;
-			boolean merge_clt = true;
-			boolean save_clt = true;
+			// see if it exists
+			/*
+			QuadCLT center_CLT = QuadCLT.restoreCenterClt(
+					quadCLTs[ref_index]); // QuadCLT ref_clt);
+			*/
+			int sensor_mask_clt = -1; // all
+			if (center_CLT == null) {
+				boolean  condition_dsi = true;
+				if (combo_dsn_final == null) {
+					combo_dsn_final =quadCLTs[ref_index].restoreComboDSI(true); // also sets quadCLTs[ref_index].dsi and blue sky
+				}
+				center_CLT = Cuas.createCenterClt( // assuming cuas_rotation is true
+						clt_parameters,  // CLTParameters  clt_parameters,
+						quadCLTs,        // QuadCLT []     quadCLTs,
+						ref_index,       // int            ref_index,
+						null,            // int     []     range, // or null
+						combo_dsn_final, // double  [][]   ref_dsi, // DSI data for the reference scene (or null to read it from file)
+						condition_dsi,   // boolean        condition_dsi,
+						sensor_mask_clt, // int            sensor_mask, // -1 - all;
+						debugLevel);     // int            debugLevel)
+				if (debugLevel > -3) {
+					System.out.println("Created/saved center CLT "+center_CLT.getImagePath());
+				}
+			} else {
+				if (debugLevel > -3) {
+					System.out.println("Read existing CLT from "+center_CLT.getImagePath());
+				}
+			}
+			// just for verification
+			boolean apply_clt = true;
+			boolean show_clt = true;
+			int [] whc = new int[3];
+			if (apply_clt) { // set GPU with data
+				center_CLT.setQuadClt(); //2025 mark that GPU is used for center_CLT
+//				quadCLTs[ref_index].setComboToTD(
+				center_CLT.setComboToTD(
+						null, // new float [][] {center_CLT.getCenterClt()}, // ,combo_seq_clt,   // final float [][]  fclt,
+						true, // merge_clt,       // final boolean     merge_channels, // duplicate same data to all selected channels
+						sensor_mask_clt, // final int         sensor_mask, // only if merge_channels
+						whc,             // final int []      whc, // if int[2], will return width, height
+						false);          // final boolean     use_reference);
+				if (show_clt) {
+					String suffix="-virtual";
+//					ImagePlus imp_virtual = ref_clt.renderFromTD ( // do we need to update gpuQuad ?
+					ImagePlus imp_virtual = center_CLT.renderFromTD ( // do we need to update gpuQuad ?
+							sensor_mask_clt,         // final int         sensor_mask,
+			        		true, // merge_clt,               // boolean             merge_channels,
+			                clt_parameters,                                 // CLTParameters clt_parameters,
+//			                clt_parameters.getColorProcParameters(ref_clt.isAux()), //ColorProcParameters colorProcParameters,
+			                clt_parameters.getColorProcParameters(center_CLT.isAux()), //ColorProcParameters colorProcParameters,
+			                clt_parameters.getRGBParameters(),              //EyesisCorrectionParameters.RGBParameters rgbParameters,\
+			                whc,                                            // null, // int []  wh,
+			                false, // toRGB,                                // boolean toRGB,
+			                false, // use_reference,                        // boolean use_reference
+			                suffix);                                        // String  suffix)
+					imp_virtual.show();
+				}
+			}
+
+			
+			
+		}		
+		/*
+		if (combine_clt) {
+			boolean apply_clt =    true; // set GPU with data
+			boolean show_clt =     true;
+			boolean merge_clt =    true;
+			boolean save_clt =     true;
+			boolean save_dsi =     true; // translate each layer of *-INTER-INTRA-LMA.tiff to the virtual center view and save it
+			boolean save_in_ref =  true; // save in reference scene directory 
+			boolean save_in_last = true; // save in last scene directory
 			QuadCLT ref_clt = quadCLTs[ref_index];
 			int sensor_mask_clt = -1; // all
 			
@@ -5939,8 +5864,21 @@ public class OpticalFlow {
         			quadCLTs[ref_index], // QuadCLT        scene,
         			debugLevel);			
         	boolean save_weights = true;
+        	
+            final double [][] ref_pXpYD = getRefPxPyD( // Use to interpolate disparity layers
+            		clt_parameters, // CLTParameters  clt_parameters,
+            		combine_clt,    // boolean        mode_cuas,            		
+//            		save_weights,   // final boolean  save_weights, // output corresponding weights for each data
+            		null,           // Rectangle      fov_tiles,
+            		ZERO3,          // double []      stereo_xyz, // offset reference camera {x,y,z}
+            		cuas_atr,       // double []      stereo_atr_in, // offset reference orientation (cuas)
+            		ds_clt[0],      // double []      ref_disparity,			
+            		ref_clt,        // QuadCLT        refCLT, // should be the same instance if one of quadCLTs
+            		debugLevel);    // int            debugLevel)
+        	
 			float [][] combo_seq_clt_w = getTDComboSceneSequence(
 					clt_parameters,      // CLTParameters  clt_parameters,
+	        		ref_pXpYD,      // double [][]    ref_pXpYD,    		
 					save_weights,        // boolean        save_weights, // output corresponding weights for each data
 					merge_clt,           // boolean        merge_all,
 					sensor_mask_clt,     // int            sensor_mask,
@@ -5973,13 +5911,60 @@ public class OpticalFlow {
 				if (save_weights) {
 					suffix_clt+="_weights";
 				}
+				if (save_in_ref) {
 				ref_clt.saveFloatArrayInModelDirectory( // error
 						suffix_clt,     // String      suffix,
 						clt_titles,       // combo_dsn_titles_full, // null,          // String []   labels, // or null
 						combo_seq_clt_w,    // dbg_data,         // float [][] data,
 						width_clt,        // int         width,
 						height_clt);      // int         height)
+				}
+				if (save_in_last) {
+					quadCLTs[quadCLTs.length - 1].saveFloatArrayInModelDirectory( // error
+							suffix_clt,     // String      suffix,
+							clt_titles,       // combo_dsn_titles_full, // null,          // String []   labels, // or null
+							combo_seq_clt_w,    // dbg_data,         // float [][] data,
+							width_clt,        // int         width,
+							height_clt);      // int         height)
+				}
+			}
+			
+			if (save_dsi) {
+//				final double [][] combo_dsi_cli = ref_clt.restoreComboDSI (true);
+				TileProcessor tp =         quadCLTs[quadCLTs.length - 1].getTileProcessor();
+				final int transform_size = tp.getTileSize();
+				final int tilesX =         tp.getTilesX();
+				final int tilesY =         tp.getTilesY();
+				final double [][] combo_dsi_center = interpolateDSI(
+						ref_pXpYD,      // final double [][] ref_pXpYD,
+						combo_dsn_final, // combo_dsi_cli,  // final double [][] dsi_in,
+						transform_size, // final int         transform_size,
+						tilesX,         // final int         tilesX,
+						tilesY);        //final int         tilesY);
+				
+				String rslt_suffix = "-CENTER-INTER-INTRA";
+				// use quadCLTs[quadCLTs.length - 1].restoreComboDSI ("-CENTER",true); to restore
+				
+				rslt_suffix += (clt_parameters.correlate_lma?"-LMA":"-NOLMA");
+		    	if (save_in_ref) {
+		    		ref_clt.saveDoubleArrayInModelDirectory( // error
+		    				rslt_suffix,           // String      suffix,
+		    				OpticalFlow.COMBO_DSN_TITLES, // combo_dsn_titles_full, // null,          // String []   labels, // or null
+		    				combo_dsi_center,       // dbg_data,         // double [][] data,
+		    				tilesX,                // int         width,
+		    				tilesY);               // int         height)
+		    	}
+		    	if (save_in_last) {
+		    		quadCLTs[quadCLTs.length - 1].saveDoubleArrayInModelDirectory( // error
+		    				rslt_suffix,           // String      suffix,
+		    				OpticalFlow.COMBO_DSN_TITLES, // combo_dsn_titles_full, // null,          // String []   labels, // or null
+		    				combo_dsi_center,       // dbg_data,         // double [][] data,
+		    				tilesX,                // int         width,
+		    				tilesY);               // int         height)
+		    	}
 			}
+			
+			
 
 			int [] whc = new int[3];
 			if (apply_clt) { // set GPU with data
@@ -6005,7 +5990,7 @@ public class OpticalFlow {
 				}
 			}
 		}
-		
+		*/
 		
 		
 		if (generate_egomotion) {
@@ -7957,7 +7942,7 @@ public class OpticalFlow {
 			double      around_sigma =        4.0;
 			int         num_virtual_refines = 2;
 			String      debugSuffix=         null; // "virtual";
-			ref_pXpYD= transformFromVirtual(
+			ref_pXpYD= Cuas.transformFromVirtual(
 					ref_disparity, // double []         disparity_ref,
 					stereo_xyz, // final double []   scene_xyz,         // camera center in world (reference) coordinates
 					stereo_atr, // final double []   scene_atr,         // camera orientation relative to world (reference) frame
@@ -8199,217 +8184,6 @@ public class OpticalFlow {
     	return imp_scenes;
     }
     
-    public static float [][] getTDComboSceneSequence(
-    		CLTParameters  clt_parameters,
-    		final boolean  save_weights, // output corresponding weights for each data
-    		final boolean  merge_all,
-    		final int      sensor_mask,
-    		Rectangle      fov_tiles,
-    		double []      stereo_xyz, // offset reference camera {x,y,z}
-    		double []      stereo_atr_in, // offset reference orientation (cuas)
-    		double []      ref_disparity,			
-    		QuadCLT []     quadCLTs,
-    		QuadCLT        refCLT, // should be the same instance if one of quadCLTs
-//    		int            ref_index,
-    		int            debugLevel) {
-    	double [] stereo_atr =        (stereo_atr_in != null)? stereo_atr_in: ZERO3; // maybe later play with rotated camera
-    	boolean mode_cuas =           (stereo_atr[0] != 0) || (stereo_atr[1] != 0) || (stereo_atr[2] != 0);
-		boolean mb_en =               clt_parameters.imp.mb_en && (fov_tiles==null);
-		double  mb_tau =              clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
-		double  mb_max_gain =         clt_parameters.imp.mb_max_gain; // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
-		ErsCorrection ers_reference = refCLT.getErsCorrection();
-        int dbg_scene = -95;
-        double [][] ref_pXpYD;
-        double [][] ref_pXpYD_or_null = null; // debugging cuas mode keeping old
-        if (mode_cuas) { //  && (dbg_scene > 0)) {
-			int         around =              2;
-			double      around_sigma =        4.0;
-			int         num_virtual_refines = 2;
-			String      debugSuffix=         null; // "virtual";
-			ref_pXpYD= transformFromVirtual(
-					ref_disparity, // double []         disparity_ref,
-					stereo_xyz, // final double []   scene_xyz,         // camera center in world (reference) coordinates
-					stereo_atr, // final double []   scene_atr,         // camera orientation relative to world (reference) frame
-					refCLT, // quadCLTs[ref_index], // final QuadCLT     reference_QuadClt,
-					around, // final int         around,            // 2 search around for interpolation
-					around_sigma, // final double      sigma,
-					num_virtual_refines, // final int         num_refines,
-					debugSuffix); // final String      debugSuffix)
-			ref_pXpYD_or_null = ref_pXpYD;
-//			quadCLTs[ref_index].getErsCorrection().setupERS();
-			refCLT.getErsCorrection().setupERS();
-			System.out.println("Calculated virtual_PxPyD");
-        } else {
-        	ref_pXpYD = transformToScenePxPyD( // now should work with offset ref_scene
-        			fov_tiles,            // final Rectangle [] extra_woi,    // show larger than sensor WOI (or null)
-        			ref_disparity,        // final double []   disparity_ref, // invalid tiles - NaN in disparity
-        			ZERO3, // stereo_xyz, // ZERO3,                // final double []   scene_xyz, // camera center in world coordinates
-        			ZERO3, // stereo_atr, // ZERO3,                // final double []   scene_atr, // camera orientation relative to world frame
-        			refCLT, // quadCLTs[ref_index],  // final QuadCLT     scene_QuadClt,
-        			refCLT, // quadCLTs[ref_index],  // final QuadCLT     reference_QuadClt, // now - may be null - for testing if scene is rotated ref
-        			QuadCLT.THREADS_MAX);          // int               threadsMax)
-        }
-        //  quadCLTs[0].getCltLength(boolean use_ref)
-        int sc0 = -1;
-		for (int nscene =  0; nscene < quadCLTs.length ; nscene++) if (quadCLTs[nscene] != null){
-			sc0 = nscene;
-			break;
-		}
-		
-		final int num_slices = merge_all? 1 : quadCLTs[sc0].getNumSensors();
-        final float [][] sumFclt = new float [num_slices * (save_weights? 2:1)][quadCLTs[sc0].getCltSize(false)];
-        final int [][] numAcc = new int [num_slices][sumFclt[0].length];
-        // next two to improve multithreading performance
-        final int tile_size_td = 4 * GPUTileProcessor.DTT_SIZE * GPUTileProcessor.DTT_SIZE;
-        final int tiles_td = sumFclt[0].length/tile_size_td;
-        final int tiles_td_all = tiles_td * num_slices; // usually sumFclt.length==1
-        
-		final Thread[] threads = ImageDtt.newThreadArray(THREADS_MAX);
-		final AtomicInteger ai = new AtomicInteger(0);
-		for (int nscene =  0; nscene < quadCLTs.length ; nscene++) if (quadCLTs[nscene] != null){
-			if (nscene== dbg_scene) {
-				System.out.println("renderSceneSequence(): nscene = "+nscene);
-			}
-			String ts = quadCLTs[nscene].getImageName();
-			double []   scene_xyz = ZERO3;
-			double []   scene_atr = ZERO3;
-//			if (nscene != ref_index) { // Check even for raw, so video frames will match in all modes 
-			if (quadCLTs[nscene] != refCLT) { // Check even for raw, so video frames will match in all modes 
-				scene_xyz = ers_reference.getSceneXYZ(ts);
-				scene_atr = ers_reference.getSceneATR(ts);
-				if ((scene_atr==null) || (scene_xyz == null)) {
-					continue;
-				}
-				double []   scene_ers_xyz_dt = ers_reference.getSceneErsXYZ_dt(ts);
-				double []   scene_ers_atr_dt = ers_reference.getSceneErsATR_dt(ts);
-				quadCLTs[nscene].getErsCorrection().setErsDt(
-						scene_ers_xyz_dt, // double []    ers_xyz_dt,
-						scene_ers_atr_dt); // double []    ers_atr_dt)(ers_scene_original_xyz_dt);
-			}
-			if (!mode_cuas && (stereo_xyz != null)) { // offset all, including reference scene - now always, it is never null
-				double [][] combo_xyzatr = ErsCorrection.combineXYZATR(
-						stereo_xyz,  // double [] reference_xyz,
-						stereo_atr,  // double [] reference_atr, 
-						scene_xyz,   // double [] scene_xyz,
-						scene_atr);  // double [] scene_atr) 
-				scene_xyz = combo_xyzatr[0];
-				scene_atr = combo_xyzatr[1];
-			}
-			int sm = merge_all? -1: sensor_mask;
-			float [][] fclt = null;
-			double [][] dxyzatr_dt = null;
-			// should get velocities from HashMap at reference scene from timestamp , not re-calculate.
-			if (mb_en) {
-//				dxyzatr_dt = getVelocities(
-//						quadCLTs, // QuadCLT []     quadCLTs,
-//						nscene);  // int            nscene)
-				dxyzatr_dt = new double[][] { // for all, including ref
-					quadCLTs[nscene].getErsCorrection().getErsXYZ_dt(),
-					quadCLTs[nscene].getErsCorrection().getErsATR_dt()};				
-			}
-			
-			
-			if (mb_en && (dxyzatr_dt != null)) {
-				double [][] motion_blur = getMotionBlur(
-						refCLT, // quadCLTs[ref_index],   // QuadCLT        ref_scene,
-						quadCLTs[nscene],      // QuadCLT        scene,         // can be the same as ref_scene
-						ref_pXpYD,             // double [][]    ref_pXpYD,     // here it is scene, not reference!
-						scene_xyz,             // double []      camera_xyz,
-						scene_atr,             // double []      camera_atr,
-						dxyzatr_dt[0],         // double []      camera_xyz_dt,
-						dxyzatr_dt[1],         // double []      camera_atr_dt,
-						0,                     // int            shrink_gaps,  // will gaps, but not more that grow by this
-						debugLevel);           // int            debug_level)
-				fclt = QuadCLT.getTDCombo(
-						sm,                  // final int         sensor_mask,
-						merge_all,           // final boolean     merge_channels,
-						null,                // final Rectangle   full_woi_in,      // show larger than sensor WOI (or null)
-						clt_parameters,      // CLTParameters     clt_parameters,
-						ref_disparity,       // double []         disparity_ref,
-						ref_pXpYD_or_null,   // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
-						// motion blur compensation 
-						mb_tau,              // double            mb_tau,      // 0.008; // time constant, sec
-						mb_max_gain,         // double            mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
-						motion_blur,         // double [][]       mb_vectors,  //
-						scene_xyz,           // final double []   scene_xyz, // camera center in world coordinates
-						scene_atr,           // final double []   scene_atr, // camera orientation relative to world frame
-						quadCLTs[nscene],    // final QuadCLT     scene,
-						refCLT, // quadCLTs[ref_index], // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
-						true,                // final boolean     show_nan,
-						QuadCLT.THREADS_MAX, // int               threadsMax,
-						debugLevel);         // final int         debugLevel) 
-
-			} else {
-				fclt = QuadCLT.getTDCombo(
-						sm,                  // final int         sensor_mask,
-						merge_all,           // final boolean     merge_channels,
-						null,                // final Rectangle   full_woi_in,      // show larger than sensor WOI (or null)
-						clt_parameters,      // CLTParameters     clt_parameters,
-						ref_disparity,       // double []         disparity_ref,
-						ref_pXpYD_or_null,   // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
-						// motion blur compensation 
-						0,                   // double            mb_tau,      // 0.008; // time constant, sec
-						0,                   // double            mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
-						null,                // double [][]       mb_vectors,  //
-						scene_xyz,           // final double []   scene_xyz, // camera center in world coordinates
-						scene_atr,           // final double []   scene_atr, // camera orientation relative to world frame
-						quadCLTs[nscene],    // final QuadCLT     scene,
-						refCLT, // quadCLTs[ref_index], // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
-						true,                // final boolean     show_nan,
-						QuadCLT.THREADS_MAX, // int               threadsMax,
-						debugLevel);         // final int         debugLevel) 
-			}
-			final float [][] ffclt = fclt;
-			ai.set(0);
-			for (int ithread = 0; ithread < threads.length; ithread++) {
-				threads[ithread] = new Thread() {
-					public void run() {
-						for (int nTileAll = ai.getAndIncrement(); nTileAll < tiles_td_all; nTileAll = ai.getAndIncrement()) {
-							int ntile = nTileAll % tiles_td;
-							int nsens = nTileAll / tiles_td;
-							int indx0 = ntile * tile_size_td;
-							int indx1 = indx0 + tile_size_td;
-							for (int indx = indx0; indx < indx1; indx++) {
-								float d =ffclt[nsens][indx];
-								if (!Float.isNaN(d)){
-									sumFclt[nsens][indx] += d;
-									numAcc[nsens][indx]++ ;
-								}
-							}
-						}
-					}
-				};
-			}		      
-			ImageDtt.startAndJoin(threads);
-		}
-		// calculate averages
-		ai.set(0);
-		for (int ithread = 0; ithread < threads.length; ithread++) {
-			threads[ithread] = new Thread() {
-				public void run() {
-					for (int nTileAll = ai.getAndIncrement(); nTileAll < tiles_td_all; nTileAll = ai.getAndIncrement()) {
-						int ntile = nTileAll % tiles_td;
-						int nsens = nTileAll / tiles_td;
-						int indx0 = ntile * tile_size_td;
-						int indx1 = indx0 + tile_size_td;
-						for (int indx = indx0; indx < indx1; indx++) {
-							if (numAcc[nsens][indx] > 0) {
-								sumFclt[nsens][indx]/=numAcc[nsens][indx];
-							} else {
-								sumFclt[nsens][indx] = Float.NaN;
-							}
-							if (save_weights) {
-								sumFclt[nsens + num_slices][indx] = numAcc[nsens][indx];
-							}
-						}
-					}
-				}
-			};
-		}		      
-		ImageDtt.startAndJoin(threads);
-    	return sumFclt;
-    }
-    
     
     public static double [][] getSceneSZXY(
     		QuadCLT   scene,
@@ -10578,7 +10352,7 @@ public class OpticalFlow {
 // And NaN in strength cause (at least) NaN for double avg_z = quadCLTs[ref_index].getAverageZ(true);		
 //		public static int COMBO_DSN_INDX_STRENGTH =    1; // strength, FG
 //		public static int COMBO_DSN_INDX_STRENGTH_BG = 6; // background strength
-		for (int slice:new int[] {COMBO_DSN_INDX_STRENGTH,COMBO_DSN_INDX_STRENGTH_BG}) {
+		for (int slice: COMBO_DSN_NONNAN) { // new int[] {COMBO_DSN_INDX_STRENGTH,COMBO_DSN_INDX_STRENGTH_BG}) {
 			if (combo_dsn_final[slice] != null) {
 				for (int i = 0; i <combo_dsn_final[slice].length; i++) {
 					if (Double.isNaN(combo_dsn_final[slice][i])) {
@@ -11504,7 +11278,6 @@ public class OpticalFlow {
 					"-results-nonoise" + (read_nonoise_lma?"-lma":"-nolma"), // String      suffix,
 					combo_dsn_titles.length - 1, // 4
 					null); // int []      wh);
-
 		}
 		
 		
@@ -13878,6 +13651,7 @@ public class OpticalFlow {
 					QuadCLT.THREADS_MAX);                  // final int                 threadsMax)  // maximal number of threads to launch
 			ref_scene.saveQuadClt(); // to re-load new set of Bayer images to the GPU (do nothing for CPU) and Geometry
 			image_dtt.setReferenceTD( // tp_tasks_ref will be updated
+					null, 			            // final float [][]          fclt, 
 					erase_clt,
 					null,                       // final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 					clt_parameters.img_dtt,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others

src/main/java/com/elphel/imagej/tileprocessor/QuadCLT.java

@@ -1464,6 +1464,7 @@ public class QuadCLT extends QuadCLTCPU {
 	    			debugLevel);                // final int                 globalDebugLevel);
 	    } else {
 	    	image_dtt.setReferenceTD( // change to main?
+					null, 			            // final float [][]          fclt, 
 	    			erase_clt, //final int                 erase_clt,
 	    			wh, // null,                       // final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 	    			clt_parameters.img_dtt,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
@@ -1608,7 +1609,7 @@ public class QuadCLT extends QuadCLTCPU {
 			final boolean     merge_channels,
 			final Rectangle   full_woi_in,      // show larger than sensor WOI in tiles (or null)
 			CLTParameters     clt_parameters,
-			double []         disparity_ref,
+			double []         disparity_ref, // may be null if ref_pXpYD!=null
 			double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
   			// motion blur compensation 
 			double            mb_tau,      // 0.008; // time constant, sec
@@ -1667,8 +1668,8 @@ public class QuadCLT extends QuadCLTCPU {
 			int dbg_width = rendered_width/GPUTileProcessor.DTT_SIZE;
 			int dbg_height = pXpYD.length/dbg_width;
 			String [] dbg_titles = (mb_vectors!=null)?
-					(new String[] {"pX","pY","Disparity","mb_X","mb_Y","disparity_ref"}):
-						(new String[] {"pX","pY","Disparity","disparity_ref"});
+					(new String[] {"pX","pY","Disparity","mb_X","mb_Y"}):
+						(new String[] {"pX","pY","Disparity"});
 			double [][] dbg_img = new double [dbg_titles.length][pXpYD.length]; // 3 + ((mb_vectors!=null)? 2:0)][pXpYD.length];
 			for (int i = 0; i < dbg_img.length; i++) {
 				Arrays.fill(dbg_img[i], Double.NaN);
@@ -1685,7 +1686,7 @@ public class QuadCLT extends QuadCLTCPU {
 					}
 				}
 			}
-			dbg_img[dbg_img.length-1] = disparity_ref;
+//			dbg_img[dbg_img.length-1] = disparity_ref;
 			ShowDoubleFloatArrays.showArrays( // out of boundary 15
 					dbg_img,
 					dbg_width,
@@ -1756,7 +1757,8 @@ public class QuadCLT extends QuadCLTCPU {
 	    			debugLevel);                // final int                 globalDebugLevel);
 	    } else {
 	    	image_dtt.setReferenceTD( // change to main?
-	    			erase_clt, //final int                 erase_clt,
+					null, 			            // final float [][]          fclt, 
+					erase_clt, //final int                 erase_clt,
 	    			wh, // null,                       // final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 	    			clt_parameters.img_dtt,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
 	    			use_reference, // true, // final boolean             use_reference_buffer,
@@ -1923,6 +1925,7 @@ public class QuadCLT extends QuadCLTCPU {
 	    			debugLevel);                // final int                 globalDebugLevel);
 	    } else {
 	    	image_dtt.setReferenceTD( // change to main?
+					null, 			            // final float [][]          fclt, 
 	    			erase_clt, //final int                 erase_clt,
 	    			wh, // null,                       // final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 	    			clt_parameters.img_dtt,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
@@ -2117,6 +2120,7 @@ public class QuadCLT extends QuadCLTCPU {
 					debugLevel);                // final int                 globalDebugLevel);
 		} else {
 			image_dtt.setReferenceTD( // change to main?
+					null, 			            // final float [][]          fclt, 
 					erase_clt, //final int                 erase_clt,
 					null, // wh, // null,                       // final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 					clt_parameters.img_dtt,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
@@ -2408,6 +2412,7 @@ public class QuadCLT extends QuadCLTCPU {
 					debugLevel);                // final int                 globalDebugLevel);
 		} else {
 			image_dtt.setReferenceTD( // change to main?
+					null, 			            // final float [][]          fclt, 
 					erase_clt, //final int                 erase_clt,
 					null, // wh, // null,                       // final int []              wh,               // null (use sensor dimensions) or pair {width, height} in pixels
 					clt_parameters.img_dtt,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
@@ -2744,11 +2749,15 @@ public class QuadCLT extends QuadCLTCPU {
 	 * @param use_reference set the reference buffers (false - set main buffers)
 	 */
 	void setComboToTD(
-			final float [][]  fclt,
-			final boolean     merge_channels, // duplicate same data to all selected channels
-			final int         sensor_mask, // only if merge_channels
-			final int []      whc, // if int[2], will return width, height
-			final boolean     use_reference){
+			float [][]  fclt,
+			boolean     merge_channels, // duplicate same data to all selected channels
+			int         sensor_mask, // only if merge_channels
+			int []      whc, // if int[2], will return width, height
+			boolean     use_reference){
+		if (fclt == null) {
+			fclt = new float[][] {getCenterClt()};
+			merge_channels = true;
+		}
 		final int [] width_height = gpuQuad.getWH(use_reference);
 		final int num_colors = gpuQuad.getNumColors();
 		if (whc != null) {

src/main/java/com/elphel/imagej/tileprocessor/QuadCLTCPU.java

@@ -132,6 +132,9 @@ public class QuadCLTCPU {
 	public static final boolean   USE_PRE_2021 = false; // temporary
 	public static final int       THREADS_MAX =       100;
 	public static final double [] ZERO3 = {0.0,0.0,0.0};
+	
+	public static final String CENTER_DIR_SUFFIX = "-CENTER";
+	public static final String CENTER_CLT_SUFFIX = "-CLT";
 
 //	public GPUTileProcessor.GpuQuad gpuQuad =      null;
 
@@ -206,6 +209,107 @@ public class QuadCLTCPU {
     public String      timestamp_first =     null;
     public String      timestamp_last =      null;
     
+    // combined clt for center view, only used in cuas mode
+    public float []    center_clt =           null; // clt data (single accumulated channel)
+    public int   []    center_clt_num =       null; // same size as combo_clt - number of accumulated scenes for this point
+    public int         center_clt_width =     0;    // image width to save combo_clt (and combo_clt_num)  as image
+    
+    public String getImagePath() {
+    	return image_path;
+    }
+    public void setImagePath(String path) {
+    	image_path = path;
+    }
+    
+    public boolean hasCenterClt() {
+    	return (center_clt != null) && (center_clt_num != null) && (center_clt_width != 0);
+    }
+    
+    public void setCenterClt(
+    		float [] clt,
+    		int []   clt_num,
+    		int      clt_width) {
+    	center_clt =       clt;
+    	center_clt_num =   clt_num;
+    	center_clt_width = clt_width;
+    }
+    public float [] getCenterClt() {
+    	return center_clt;
+    }
+    public int [] getCenterCltNum() {
+    	return center_clt_num;
+    }
+    public int [] getCenterCltWH() {
+    	return new int [] {center_clt_width, center_clt.length / center_clt_width};
+    }
+    
+	public boolean saveCenterClt() {
+		if (!hasCenterClt()) {
+			System.out.println("saveCenterCli(): center CLT data is not set");
+			return false;
+		}
+		float [] fcenter_clt_num = new float [center_clt_num.length];
+		for (int i = 0; i <center_clt_num.length; i++) {
+			fcenter_clt_num[i] = center_clt_num[i];
+		}
+		String [] clt_titles= {"CLT","num-acc"};
+		float [][] fdata = {center_clt, fcenter_clt_num};
+		saveFloatArrayInModelDirectory( // error
+				CENTER_CLT_SUFFIX,    // String      suffix,
+				clt_titles,           // String []   labels, // or null
+				fdata,                // dbg_data,         // float [][] data,
+				getCenterCltWH()[0],  // int         width,
+				getCenterCltWH()[1]); // int         height)
+		return true;
+	}
+    
+	public static String getCenterDirName(String ref_name) {
+		return ref_name + CENTER_DIR_SUFFIX;
+	}
+	
+	public static QuadCLT restoreCenterClt(QuadCLT ref_clt) {
+		String center_name =  QuadCLT.getCenterDirName(ref_clt.getImageName());
+		String ref_dir_path = ref_clt.getX3dDirectory(center_name);
+		File cdir = new File(ref_dir_path);
+		if (!cdir.exists() || !cdir.isDirectory()) {
+			System.out.println("restoreCenterClt(): directory does not exist: "+ref_dir_path);
+			return null;
+		}
+		QuadCLT center_CLT = new QuadCLT(ref_clt, center_name);
+		center_CLT.setImagePath(cdir.getPath());		
+		int [] wh = new int[2];
+		float [][] fclt_w = center_CLT.readFloatArrayFromModelDirectory(
+				CENTER_CLT_SUFFIX,  // String      suffix,
+				0,                  // int         num_slices, // (0 - all)
+				wh);                // int []      wh)
+		if (fclt_w != null) {
+			if (fclt_w.length != 2) {
+				System.out.println("restoreCenterClt(): expected a 2-slice data, got "+fclt_w.length+" slices.");
+				return null;
+			}
+			int [] num = new int [fclt_w[1].length];
+			for (int j = 0; j < num.length; j++) {
+				num[j] = (int) Math.round(fclt_w[1][j]);
+			}
+			center_CLT.setCenterClt(
+					fclt_w[0], // float [] clt,
+					num,       // int []   clt_num,
+					wh[0]);    // int      clt_width);
+		}
+		// TwoQuadCLT.DSI_COMBO_SUFFIX
+		/*
+		center_CLT.restoreDSI(
+				TwoQuadCLT.DSI_COMBO_SUFFIX, // String suffix,
+	    		true); // boolean silent)
+	    		*/
+		center_CLT.restoreComboDSI(true);
+		
+		
+		return center_CLT;
+	}
+	
+	
+	
     public String getReferenceTimestamp() {
     	return timestamp_reference;
     }
@@ -509,7 +613,8 @@ public class QuadCLTCPU {
 				ref_index,      // final int           ref_index,
 				null,           // double [][][]       dxyzatr,
 				earliest_scene, // final int           early_index,
-				last_index      //(quadCLTs.length -1) // int           last_index,
+				last_index,     //(quadCLTs.length -1) // int           last_index,
+				debugLevel      // final int           debugLevel
 	    		);
 		double [][][] xyzatr = new double [quadCLTs.length][][];
 		ErsCorrection ers_ref = quadCLTs[ref_index].getErsCorrection();
@@ -848,7 +953,9 @@ public class QuadCLTCPU {
     				ref_index,      // final int           ref_index,
     				null,           // double [][][]       dxyzatr,
     				early_index,    // final int           early_index,
-    				last_index);    //(quadCLTs.length -1) // int           last_index,
+    				last_index,     //(quadCLTs.length -1) // int           last_index,
+    				debugLevel);    // final int           debugLevel
+
     				
     	}
     	double [][][] xyzatr_pull = new double [xyzatr.length][][];
@@ -1202,7 +1309,8 @@ public class QuadCLTCPU {
     		final int           ref_index,
     		double [][][]       dxyzatr,
     		final int           early_index,
-    		final int           last_index){
+    		final int           last_index,
+    		final int           debugLevel){
     	boolean renormalize = true;
     	double [] timestamps = new double [quadCLTs.length];
     	for (int nscene = early_index; nscene <= last_index; nscene++) {
@@ -1268,7 +1376,20 @@ public class QuadCLTCPU {
     			prev_scene = nscene;
     		}
     	}
-    	return            xyzatr;
+    	if (debugLevel>0) {
+    		System.out.println(String.format(
+    				"%3s\t%9s\t%9s\t%9s\t%9s\t%9s\t%9s",
+    				"N","X","Y","Z","A","T","R"));
+        	for (int nscene = 0; nscene < quadCLTs.length; nscene++) if (quadCLTs[nscene] != null) {
+           		System.out.println(String.format(
+        				"%3d\t%9.5f\t%9.5f\t%9.5f\t%9.5f\t%9.5f\t%9.5f",
+        				nscene,
+        				xyzatr[nscene][0][0],xyzatr[nscene][0][1],xyzatr[nscene][0][2],
+        				xyzatr[nscene][1][0],xyzatr[nscene][1][1],xyzatr[nscene][1][2]
+        				));
+        	}    		
+    	}    	
+    	return xyzatr;
     }
 
     /**
@@ -1292,8 +1413,10 @@ public class QuadCLTCPU {
     	double [] ims_mount_atr = clt_parameters.imp.getImsMountATR(); // converts to radians
     	double [][][] scenes_xyzatr = new double [quadCLTs.length][][];
     	//		scenes_xyzatr[ref_index] =         new double[2][3]; // all zeros
-
     	Did_ins_2   d2_ref = did_ins_2;
+    	if (d2_ref == null) {
+    		d2_ref=quadCLTs[quadCLTs.length-1].did_ins_2;
+    	}
     	// apply correction to orientation - better not to prevent extra loops.
     	// And it only applies to rotations
     	double [] cam_quat_ref_enu =Imx5.quaternionImsToCam(d2_ref.getQEnu() ,
@@ -1336,7 +1459,7 @@ public class QuadCLTCPU {
     	}
     	if (debugLevel>0) {
     		System.out.println(String.format(
-    				"%3s\t%9s\t%9s\t%9s\t%9s\t%9s",
+    				"%3s\t%9s\t%9s\t%9s\t%9s\t%9s\t%9s",
     				"N","X","Y","Z","A","T","R"));
         	for (int nscene = 0; nscene < quadCLTs.length; nscene++) if (quadCLTs[nscene] != null) {
            		System.out.println(String.format(
@@ -3575,6 +3698,10 @@ public class QuadCLTCPU {
     }
     public boolean[] getBooleanBlueSky() {
     	double [] double_blue_sky = getDoubleBlueSky();
+    	return getBooleanBlueSky(double_blue_sky);
+    }
+
+    public static boolean[] getBooleanBlueSky(double [] double_blue_sky) {
     	if (double_blue_sky == null) {
     		return null;
     	}
@@ -3584,7 +3711,9 @@ public class QuadCLTCPU {
     	}
     	return blue_sky;
     }
-
+    
+    
+    
     public boolean hasBlueSky() {
     	return (getDoubleBlueSky() != null);
 
@@ -3835,7 +3964,7 @@ public class QuadCLTCPU {
      * @return combo DSI if read, null if failed to read. Result has full length
      *               (OpticalFlow.COMBO_DSN_TITLES.length), missing slices are null
      */
-    public double [][] restoreComboDSI (boolean silent) {
+    /*public double [][] restoreComboDSI (boolean silent) {
     	double [][] combo_dsi = new double [OpticalFlow.COMBO_DSN_TITLES.length][];
     	for (int indx: new int[]{INDEX_INTER_LMA, INDEX_INTER}) {
     		int num_slices = restoreDSI(
@@ -3850,7 +3979,32 @@ public class QuadCLTCPU {
     	}
     	return null;
     }
-
+    */
+    public double [][] restoreComboDSI (boolean silent) {
+    	return restoreComboDSI (null, silent);
+    }
+    
+    public double [][] restoreComboDSI (
+    		String extra_suffix,
+    		boolean silent) {
+    	if (extra_suffix == null) {
+    		extra_suffix = "";
+    	}
+    	double [][] combo_dsi = new double [OpticalFlow.COMBO_DSN_TITLES.length][];
+    	for (int indx: new int[]{INDEX_INTER_LMA, INDEX_INTER}) {
+    		int num_slices = restoreDSI(
+    				extra_suffix+DSI_SUFFIXES[indx], // String suffix, // "-DSI_COMBO", "-DSI_MAIN" (DSI_COMBO_SUFFIX, DSI_MAIN_SUFFIX)
+    				combo_dsi,          // double [][] dsi, // if null - just check file exists
+    				silent);            // boolean silent);
+    		if (num_slices >= 0) {
+    			System.out.println ("restoreComboDSI(): used "+getX3dDirectory()+ Prefs.getFileSeparator() + image_name + extra_suffix + DSI_SUFFIXES[indx] + ".tiff");
+    			setDSIFromCombo(combo_dsi); // reformat
+    			return combo_dsi;
+    		}
+    	}
+    	return null;
+    }
+    
     public double [][] readComboDSI (boolean silent) {
     	double [][] combo_dsi = new double [OpticalFlow.COMBO_DSN_TITLES.length][];
     	for (int indx: new int[]{INDEX_INTER_LMA, INDEX_INTER}) {
@@ -3883,7 +4037,7 @@ public class QuadCLTCPU {
     			silent);
     }
     //min_ref_frac	
-    public boolean [] getReliableTiles(
+    public boolean [] getReliableTiles_old(
     		boolean   use_combo,
     		double    min_strength,
     		double    min_ref_frac,
@@ -3907,8 +4061,6 @@ public class QuadCLTCPU {
     		needs_lma = needs_lma_combo;
     	} else {
     		main_dsi =  readDsiMain();
-//    		this.dsi =  readDsiMain(); // this.dsiis used in getBooleanBlueSky()
-//    		main_dsi =  this.dsi;
     	}
     	if (main_dsi == null) {
     		return null;
@@ -3917,6 +4069,9 @@ public class QuadCLTCPU {
     	double [] disparity_lma = main_dsi[isAux()?TwoQuadCLT.DSI_DISPARITY_AUX_LMA:TwoQuadCLT.DSI_DISPARITY_MAIN_LMA];
     	double [] strength = main_dsi[isAux()?TwoQuadCLT.DSI_STRENGTH_AUX:TwoQuadCLT.DSI_STRENGTH_MAIN];
     	double [] sfm_gain = main_dsi[isAux()?TwoQuadCLT.DSI_SFM_GAIN_AUX:TwoQuadCLT.DSI_SFM_GAIN_MAIN];
+    	
+    	
+    	
     	if ((strength == null) || (needs_lma && (disparity_lma == null) )) {
     		return null;
     	}
@@ -4007,7 +4162,161 @@ public class QuadCLTCPU {
     	}
     	return reliable;
     }
+    public boolean [] getReliableTiles(
+    		boolean   use_combo,
+    		double    min_strength,
+    		double    min_ref_frac,
+    		boolean   needs_lma,
+    		boolean   needs_lma_combo,
+    		boolean   sfm_filter, 
+    		double    sfm_minmax, 
+    		double    sfm_fracmax,
+    		double    sfm_fracall,    		
+    		double [] reduced_strength, // if not null will return >0 if had to reduce strength (no change if did not reduce)
+    		int       debugLevel
+    		) {
+    	// 10.15.2023 - was error here,	readComboDSI (silent) returns combo_dsi, not converted to this.dsi format;	
+    	//		double [][] main_dsi = use_combo? readComboDSI (silent): readDsiMain();
+    	double [][] main_dsi = null;
+    	boolean silent = false;
+    	if (use_combo) {
+    		readComboDSI (silent);
+    		main_dsi = this.dsi;
+    		needs_lma = needs_lma_combo;
+    	} else {
+    		main_dsi =  readDsiMain();
+    	}
+    	if (main_dsi == null) {
+    		return null;
+    	}
+    	boolean [] blue_sky = getBooleanBlueSky();
+    	double [] disparity_lma = main_dsi[isAux()?TwoQuadCLT.DSI_DISPARITY_AUX_LMA:TwoQuadCLT.DSI_DISPARITY_MAIN_LMA];
+    	double [] strength = main_dsi[isAux()?TwoQuadCLT.DSI_STRENGTH_AUX:TwoQuadCLT.DSI_STRENGTH_MAIN];
+    	double [] sfm_gain = main_dsi[isAux()?TwoQuadCLT.DSI_SFM_GAIN_AUX:TwoQuadCLT.DSI_SFM_GAIN_MAIN];
+    	return getReliableTiles(
+    			min_strength,     // double     min_strength,
+    			min_ref_frac,     // double     min_ref_frac,
+    			needs_lma,        // boolean    needs_lma,
+    			sfm_filter,       // boolean    sfm_filter, 
+    			sfm_minmax,       // double     sfm_minmax, 
+    			sfm_fracmax,      // double     sfm_fracmax,
+    			sfm_fracall,      // double     sfm_fracall,    		
+    			reduced_strength, // double []  reduced_strength, // if not null will return >0 if had to reduce strength (no change if did not reduce)
+    			disparity_lma,    // double []  disparity_lma,
+    			strength,         // double []  strength,
+    			sfm_gain,         // double []  sfm_gain,
+    			blue_sky,         // boolean [] blue_sky,
+    			debugLevel);      // int        debugLevel)
+    }
+    
+    public static boolean [] getReliableTiles(
+    		double     min_strength,
+    		double     min_ref_frac,
+    		boolean    needs_lma,
+    		boolean    sfm_filter, 
+    		double     sfm_minmax, 
+    		double     sfm_fracmax,
+    		double     sfm_fracall,    		
+    		double []  reduced_strength, // if not null will return >0 if had to reduce strength (no change if did not reduce)
+        	double []  disparity_lma,
+        	double []  strength,
+        	double []  sfm_gain,
+        	boolean [] blue_sky,
+       		int        debugLevel) {
+    	int NUM_BINS = 1024;
+    	if ((strength == null) || (needs_lma && (disparity_lma == null) )) {
+    		return null;
+    	}
+    	int	min_reliable = (int) Math.round (strength.length * min_ref_frac);
+    	strength = strength.clone();
+    	boolean [] reliable = new boolean [strength.length];
+    	for (int i = 0; i < reliable.length; i++) {
+    		reliable[i] = (strength[i] >= min_strength) &&
+    				(!needs_lma || !Double.isNaN(disparity_lma[i]));
+    	}
+//    	boolean [] blue_sky = getBooleanBlueSky();
+    	if (blue_sky != null) {
+    		for (int i = 0; i < reliable.length; i++) {
+    			if (blue_sky[i]){
+    				reliable[i] = false;
+    				strength[i] = 0.0;
+    			}
+    		}
+    	}
+    	if (sfm_filter && (sfm_gain != null)) {
+    		double sfm_max = 0.0;
+    		for (int i = 0; i < reliable.length; i++) {
+    			if (sfm_gain[i] > sfm_max) sfm_max = sfm_gain[i];
+    		}
+    		if (sfm_max > sfm_minmax) {
+    			double sfm_thresh = sfm_max * sfm_fracmax;
+    			boolean [] sfm_good = new boolean [reliable.length];
+    			int num_sfm_gain = 0;
+    			for (int i = 0; i < reliable.length; i++) {
+    				if (sfm_gain[i] >= sfm_thresh) {
+    					sfm_good[i] = true;
+    					num_sfm_gain++;
+    				}
+    			}
+    			if (num_sfm_gain > (sfm_fracall * reliable.length)) {
+    				if (debugLevel > -3) {
+    					System.out.println("getReliableTiles(): Using SfM filter");
+    				}
+    				for (int i = 0; i < reliable.length; i++) {
+    					if (!sfm_good[i]) { // sfm_gain[i] <= 0){
+    						reliable[i] = false;
+    						strength[i] = 0.0;
+    					}
+    				}
+    			}
+    		}
+    	}
+    	int num_reliable = 0;
+    	for (boolean b: reliable) if (b) num_reliable++;
+		if (debugLevel > -3) {
+			System.out.println("getReliableTiles(): num_reliable="+num_reliable);
+		}
 
+    	if (num_reliable < min_reliable) { // not enough, select best tiles, ignoring LMA
+    		double max_str = 0;
+    		for (double s:strength) if ((s > max_str) )max_str= s; // NaN OK
+    		if (max_str == 0) return null;
+    		int [] hist = new int[NUM_BINS];
+    		int num_gt0 = 0;
+    		for (double s:strength) if (s > 0) { // verify enough > 0
+    			int bin = (int) Math.floor(NUM_BINS * s /max_str);
+    			if (bin >= NUM_BINS) bin = NUM_BINS - 1;
+    			hist[bin]++;
+    			num_gt0++;
+    		}
+    		if (num_gt0 < min_reliable) {
+    			return null; // not enough non-zero values
+    		}
+    		num_reliable = 0;
+    		int num_prev = 0;
+    		int bin = NUM_BINS - 1;
+    		for (; num_reliable < min_reliable; bin--) {
+    			num_prev =      num_reliable;
+    			num_reliable +=	hist[bin];
+    		}
+    		double threshold = (bin + 1.0 * (num_reliable - min_reliable) / (num_reliable - num_prev)) * max_str / NUM_BINS;
+    		num_reliable = 0;
+    		for (int i = 0; i < reliable.length; i++) {
+    			reliable[i] = (strength[i] > threshold);
+    			if (reliable[i]) num_reliable++;
+    		}
+    		if (reduced_strength != null) {
+    			reduced_strength[0] = threshold;
+    		}
+    		if (debugLevel > -3) {
+    			System.out.println("getReliableTiles(): updated num_reliable="+num_reliable);
+    		}
+    	}
+    	return reliable;
+    }
+    
+    
+    
     public double [][] readDsiMain(){
     	double [][] main_dsi = new double [TwoQuadCLT.DSI_SLICES.length][];
     	int slices = restoreDSI(
@@ -5245,7 +5554,21 @@ public class QuadCLTCPU {
 			int         num_slices, // (0 - all)
 			int []      wh)
 	{
-//		final int [] image_wh = geometryCorrection.getSensorWH();
+		float [][] fdata = readFloatArrayFromModelDirectory(
+				suffix,     // String      suffix,
+				num_slices, // int         num_slices, // (0 - all)
+				wh);        // int []      wh)
+		if (fdata == null) {
+			return null;
+		}
+		double [][] result = new double [fdata.length][]; 
+		for (int n = 0; n < fdata.length; n++) if (fdata[n]!=null) {
+			result[n] = new double [fdata[n].length];
+			for (int i = 0; i < fdata[n].length; i++) {
+				result[n][i] = fdata[n][i];
+			}
+		}
+		/*
 		String x3d_path = getX3dDirectory();
 		String file_name = image_name + suffix;
 		String file_path = x3d_path + Prefs.getFileSeparator() + file_name + ".tiff";
@@ -5277,9 +5600,42 @@ public class QuadCLTCPU {
 			wh[0] = imp.getWidth();
 			wh[1] = imp.getHeight();
 		}
+		*/
 		return result;
 	}
 		
+	public float [][] readFloatArrayFromModelDirectory(
+			String      suffix,
+			int         num_slices, // (0 - all)
+			int []      wh)
+	{
+		String x3d_path = getX3dDirectory();
+		String file_name = image_name + suffix;
+		String file_path = x3d_path + Prefs.getFileSeparator() + file_name + ".tiff";
+		ImagePlus imp = null;
+		try {
+			imp = new ImagePlus(file_path);
+		} catch (Exception e) {
+			System.out.println ("Failed to open "+file_path+", will generate it");
+		}
+		if ((imp == null) || (imp.getTitle() == null) || (imp.getTitle().equals(""))) {
+			return null;
+		}
+		ImageStack imageStack = imp.getStack();
+		int nChn=imageStack.getSize();
+		if ((num_slices > 0) && (num_slices < nChn)) {
+			nChn = num_slices;
+		}
+		float [][] result = new float [nChn][]; 
+		for (int n = 0; n < nChn; n++) {
+			result[n] = (float[]) imageStack.getPixels(n + 1);
+		}
+		if (wh != null) {
+			wh[0] = imp.getWidth();
+			wh[1] = imp.getHeight();
+		}
+		return result;
+	}
 	
 	
 

src/main/java/com/elphel/imagej/tileprocessor/sfm/StructureFromMotion.java

@@ -176,6 +176,7 @@ public class StructureFromMotion {
 				ref_scene,                 // QuadCLT            ref_scene,
 				mb_ref_disparity,          // double []          ref_disparity, // null or alternative reference disparity
 				ref_pXpYD,                 // double [][]        ref_pXpYD,
+				null, 			    // final float [][]          fclt, 
 				reliable_ref,              // final boolean []   selection, // may be null, if not null do not  process unselected tiles
 				margin,                    // final int          margin,
 				// motion blur compensation 
@@ -2103,6 +2104,7 @@ public class StructureFromMotion {
 				scenes[0],           // ref_scene,                 // QuadCLT            ref_scene,
 				null,                // ref_disparity,             // double []          ref_disparity, // null or alternative reference disparity
 				scenes_pXpYD[0],     // ref_pXpYD,                 // double [][]        ref_pXpYD,
+				null, 			    // final float [][]          fclt, 
 				null,                // reliable_ref,              // final boolean []   selection, // may be null, if not null do not  process unselected tiles
 				margin,              // final int          margin,
 				// motion blur compensation