package com.elphel.imagej.cuas;

import java.text.SimpleDateFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Calendar;
import java.util.concurrent.atomic.AtomicInteger;

import com.elphel.imagej.cameras.CLTParameters;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.gpu.GPUTileProcessor;
import com.elphel.imagej.gpu.GpuQuad;
import com.elphel.imagej.gpu.TpTask;
import com.elphel.imagej.ims.UasLogReader;
import com.elphel.imagej.tileprocessor.Correlation2d;
import com.elphel.imagej.tileprocessor.ErsCorrection;
import com.elphel.imagej.tileprocessor.GeometryCorrection;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.ImageDttParameters;
import com.elphel.imagej.tileprocessor.OpticalFlow;
import com.elphel.imagej.tileprocessor.QuadCLT;

import ij.ImagePlus;
import ij.Prefs;

public class CuasRanging {
	public static final String   TARGET_RANGING_LOGS_SUFFIX =      "-TARGET_RANGING.log";
	public static final String   TARGET_RANGING_LOGS_FULL_SUFFIX = "-TARGET_RANGING_FULL.log";
	public static final String   TARGET_DISPARITIES_SUFFIX =       "-TARGET_DISPARITIES";
	public static final String   TARGET_GLOBALS_SUFFIX =           "-TARGET_GLOBALS";
	public static final String   TARGET_STATS_SUFFIX =             "-TARGETS"; // *.csv

	
	final QuadCLT           center_CLT;
	final QuadCLT []        scenes;
	final double [][][]     img_um;
	final CLTParameters     clt_parameters;
	CuasMotion              cuasMotion = null;
	double [][][]           ranging_verify = null;            // just for debugging, length - total number of scenes (some will be null)
	double []               ranging_verify_weight = null;     // accumulated scene weight (divide ranging_verify pixels)
	public int              debugLevel = 0;
	
	public CuasRanging 	(
			CLTParameters     clt_parameters,
			QuadCLT           center_CLT,
			QuadCLT []        scenes,
			int               debugLevel) {
		this.clt_parameters = clt_parameters;
		this.center_CLT = center_CLT;
		this.scenes = scenes;
		this.img_um =         new double [scenes.length][][];
	}
	
	public CuasMotion detectTargets(
			UasLogReader    uasLogReader,
			boolean         batch_mode){
		boolean save_linear_cuas = true;
		boolean save_um_cuas = true;
		double  um_sigma =           clt_parameters.imp.um_sigma;
		double  um_weight =          clt_parameters.imp.um_weight;
		boolean mono_fixed =         clt_parameters.imp.mono_fixed;
		double  mono_range =         clt_parameters.imp.mono_range;
		double rng_um_sigma =        clt_parameters.imp.cuas_rng_um_sigma;
///		boolean rng_um_twice =       clt_parameters.imp.cuas_rng_um2;
		double rng_radius0 =         clt_parameters.imp.cuas_rng_radius0;
		double rng_radius =          clt_parameters.imp.cuas_rng_radius;
		double rng_blur =            clt_parameters.imp.cuas_rng_blur;
		
		boolean rng_img =            clt_parameters.imp.cuas_rng_img;
		boolean rng_glob =           clt_parameters.imp.cuas_rng_glob;
		boolean rng_disp =           clt_parameters.imp.cuas_rng_disp;
		boolean recalc_target_id = true; // temporalily bug fixing, may be removed later
		
		boolean rng_vfy =            clt_parameters.imp.cuas_rng_vfy; // Generate/save ranging verification images (per-sensor and combined rendering from the same data)
		boolean reset_disparity =    clt_parameters.imp.cuas_reset_disparity;
		
		double cuas_initial_disparity=clt_parameters.imp.cuas_initial_disparity; // 1.0;  // Start correlation with this disparity (in addition to infinity) after reset
		double  cuas_infinity =      clt_parameters.imp.cuas_infinity ; // 0.63;  // disparity at infinity for targets
		boolean smooth_omegas =      clt_parameters.imp.cuas_smooth_omegas;
		boolean log_ranging =        clt_parameters.imp.cuas_log_ranging; //  true;  // Log ranging iterations
		boolean reuse_disparity =    clt_parameters.imp.cuas_reuse_disparity;
		boolean radar_mode =         clt_parameters.imp.cuas_reuse_globals;
///		double  rng_limit =          clt_parameters.imp.cuas_rng_limit;
		
		boolean generate_output =    clt_parameters.imp.cuas_generate;      // generate and save targets Tiff and/or video files
		boolean generate_csv =       clt_parameters.imp.cuas_generate_csv;  // generate and save targets as csv file
		boolean clean_video =        clt_parameters.imp.cuas_clean_video;   //true;// save video without any debug information for targets, output in TIFF files. False - same output for video and TIFFs
		
    			
		double [][]combo_dsi = center_CLT.comboFromMain();
    	double [][] dls = { 
    			combo_dsi[OpticalFlow.COMBO_DSN_INDX_DISP], // **** null on second scene sequence
    			combo_dsi[OpticalFlow.COMBO_DSN_INDX_LMA],
    			combo_dsi[OpticalFlow.COMBO_DSN_INDX_STRENGTH]
    	};
    	double [][] ds = OpticalFlow.conditionInitialDS(
    			true,                // boolean        use_conf,       // use configuration parameters, false - use following  
    			clt_parameters,      // CLTParameters  clt_parameters,
    			dls,                 // double [][]    dls
    			center_CLT,          // quadCLTs[ref_index], // QuadCLT        scene,
    			debugLevel-1);			
    	double [] disparity_fg = ds[0]; // combo_dsn_final[COMBO_DSN_INDX_DISP_FG];
    	double [] strength_fg =  ds[1];
		if (strength_fg != null) { // for FG/BG only, fixing for transformCameraVew()
			for (int i = 0; i < disparity_fg.length; i++) {
				if (!Double.isNaN(disparity_fg[i]) && (strength_fg[i] == 0)) strength_fg[i] = 0.01; // transformCameraVew ignores strength= 0
			}
		}
    	
    	double [] xyz_offset= {0,0,0};
		double [][] ds_vantage = new double[][] {disparity_fg, strength_fg};

		boolean debug_vantage = false;
    	double [][] dbg_vantage = debug_vantage ? (new double[7][]): null;
		if (dbg_vantage != null) {
			for (int i = 0; i < 3; i++) {
				dbg_vantage[i] = dls[i].clone();
			}
			for (int i = 0; i < 2; i++) {
				dbg_vantage[i+3] = ds_vantage[i].clone();
			}
		}
		ds_vantage = OpticalFlow.transformCameraVew(
				null,              // (debug_ds_fg_virt?"transformCameraVew":null),     // final String    title,
				ds_vantage,        // final double [][] dsrbg_camera_in,
				xyz_offset,        // xyz_offset, // _inverse[0], // final double [] scene_xyz, // camera center in world coordinates
				OpticalFlow.ZERO3, // _inverse[1], // final double [] scene_atr, // camera orientation relative to world frame
				center_CLT,        // quadCLTs[ref_index],      // final QuadCLT   scene_QuadClt,
				center_CLT,        // quadCLTs[ref_index],      // final QuadCLT   reference_QuadClt,
				8); // iscale);                  // final int       iscale);
		if (dbg_vantage != null) {
			for (int i = 0; i < 2; i++) {
				dbg_vantage[i+5] = ds_vantage[i].clone();
			}
    		ShowDoubleFloatArrays.showArrays(
    				dbg_vantage,
    				center_CLT.getTileProcessor().getTilesX(),
    				center_CLT.getTileProcessor().getTilesY(),
    				true,
    				center_CLT.getImageName()+"-ds_vantage", // "-corr2d"+"-"+frame0+"-"+frame1+"-"+corr_pairs,
    				new String[] {"disp0","lma0", "str0", "disp","str","virt_disp", "virt_str"});
		}

		float [] average_pixels =    (float []) center_CLT.getCenterAverage().getProcessor().getPixels();
		float [][] average_channels = new float [][] {average_pixels}; // for future color images

		double [] cuas_atr = OpticalFlow.ZERO3;

		String scenes_suffix = center_CLT.getImageName()+"-CUAS"; // "1747829900_781803-SEQ-FG-MONO-FPN";
		ImagePlus imp_targets= OpticalFlow.renderSceneSequence(
				clt_parameters,     // CLTParameters clt_parameters,
				true,  // center_CLT.hasCenterClt(), // boolean        mode_cuas,
				false, // clt_parameters.imp.um_mono, // boolean        um_mono,
				clt_parameters.imp.calculate_average,      // boolean        insert_average, // then add new parameter, keep add average
				average_channels,     //  average_slice,
				clt_parameters.imp.subtract_average, // boolean    subtract_average,
				clt_parameters.imp.running_average,  // int            running_average,
				null, // fov_tiles,   // Rectangle     fov_tiles,
				1, // mode3d,         // int           mode3d,
				false, // toRGB,      // boolean       toRGB,
				xyz_offset,           // double []     stereo_offset, // offset reference camera {x,y,z}
				cuas_atr,             // double []      stereo_atr, // offset reference orientation (cuas)
				1, // sensor_mask,    // int           sensor_mask,
				scenes_suffix,        // String        suffix,
				ds_vantage[0],        // selected_disparity, // double []     ref_disparity,			
				scenes,               // QuadCLT []    quadCLTs,
				center_CLT,           // ref_index,          // int           ref_index,
				ImageDtt.THREADS_MAX, // threadsMax,         // int           threadsMax,
				debugLevel);          // int           debugLevel);
        if (save_linear_cuas) {	   
        	center_CLT.saveImagePlusInModelDirectory(
                    null, // "GPU-SHIFTED-D"+clt_parameters.disparity, // String      suffix,
                    imp_targets); // imp_scenes); // ImagePlus   imp)
        }
        	// always generating UM, even if not saving - needed for targets
        String cuas_title = imp_targets.getTitle();
        String um_suffix = "";
        if (mono_fixed) {
        	um_suffix = String.format("-UM%.1f_%.3f_%.0f",um_sigma,um_weight,mono_range);
        } else {
        	um_suffix = String.format("-UM%.1f_%.3f_A",um_sigma,um_weight);
        }
        imp_targets = OpticalFlow.applyUM ( // apply UM
        		cuas_title + um_suffix , // final String title, // should include -UM...
        		imp_targets,                  // final ImagePlus imp,
        		um_sigma,                  // final double um_sigma,
        		um_weight);                // final double um_weight)
        if (save_um_cuas) {
        	center_CLT.saveImagePlusInModelDirectory(
        			null,      // "GPU-SHIFTED-D"+clt_parameters.disparity, // String      suffix,
        			imp_targets); // imp_scenes); // ImagePlus   imp)
        }
		
        boolean insert_average = (center_CLT.getCenterAverage() != null) || clt_parameters.imp.calculate_average;
        System.out.println("Will generate targets images/videos");
        int first_corr = insert_average? 1:0; // skip average
        int num_scenes = imp_targets.getStack().getSize()- first_corr; // includes average
        float [][] fpixels = new float[num_scenes][];
        String [] scene_titles = new String [num_scenes];
        for (int nscene = 0; nscene < fpixels.length; nscene++) {
            fpixels[nscene] = (float[]) imp_targets.getStack().getPixels(nscene+first_corr+1);
            String s = imp_targets.getStack().getSliceLabel(nscene+first_corr+1);
            if (s.indexOf("-0") >=0) {
                s=s.substring(0, s.indexOf("-0"));
            }
            scene_titles[nscene] = s;// imp_targets.getStack().getSliceLabel(nscene+first_corr+1);
        }
        boolean skip_targets = false;// true;
        if (skip_targets) {
        	return null;
        }
		this.cuasMotion = new CuasMotion (
				clt_parameters, // CLTParameters     clt_parameters,
				scene_titles,   // String []         scene_titles,
				center_CLT,     // QuadCLT           parentCLT,
				uasLogReader,   // UasLogReader      uasLogReader,
				debugLevel);    // int               debugLevel)
		double[][][] targets; //  = null;

		/*        
		cuasMotion.processMovingTargetsMulti( // will remove rendering 
				batch_mode,     // final boolean         batch_mode,
				fpixels,        // final float [][]      fpixels,
				debugLevel);    // final int             debugLevel) {
		double[][][] targets = cuasMotion.getTargets();
*/
//		String model_prefix = center_CLT.getImageName()+CuasMotion.getParametersSuffixRanging(clt_parameters,null);
//		double [][] glob_stats = null;
        if (radar_mode) {
    		double [][][] disparity_targets = null;
			String disparity_path = center_CLT.getImageName()+CuasMotion.getParametersSuffixRanging(clt_parameters,TARGET_GLOBALS_SUFFIX); // may later be directly provided (e.g. select specific version)
			if (disparity_path.indexOf(Prefs.getFileSeparator()) < 0) { // now always
				String x3d_path = center_CLT.getX3dDirectory();
				disparity_path = x3d_path + Prefs.getFileSeparator() + disparity_path+".tiff";
			}
			disparity_targets = CuasMotion.getTargetsFromHyperAugment(disparity_path);
			if (disparity_targets == null) {
				System.out.println("processMovingTargetsMulti(): FAILED TO READ TARGET DATA FROM "+disparity_path);
				return null;
			} else {
				targets = disparity_targets;
				System.out.println("processMovingTargetsMulti(): re-using target disparities from "+disparity_path+", generating radar mode images");
			}
			cuasMotion.setTargets(targets);
        } else {
    		cuasMotion.processMovingTargetsMulti( // will remove rendering 
    				batch_mode,     // final boolean         batch_mode,
    				fpixels,        // final float [][]      fpixels,
    				debugLevel);    // final int             debugLevel) {
//    		double[][][] 
    		targets = cuasMotion.getTargets();

        	if (rng_disp || rng_glob) {
        		double [][][] disparity_targets = null;
        		if (reuse_disparity) {
        			//				String disparity_path = model_prefix+TARGET_DISPARITIES_SUFFIX; // may later be directly provided (e.g. select specific version)
        			String disparity_path = center_CLT.getImageName()+CuasMotion.getParametersSuffixRanging(clt_parameters,TARGET_DISPARITIES_SUFFIX); // may later be directly provided (e.g. select specific version)
        			if (disparity_path.indexOf(Prefs.getFileSeparator()) < 0) { // now always
        				String x3d_path = center_CLT.getX3dDirectory();
        				disparity_path = x3d_path + Prefs.getFileSeparator() + disparity_path+".tiff";
        			}
        			disparity_targets = CuasMotion.getTargetsFromHyperAugment(disparity_path);
        			if (disparity_targets == null) {
        				System.out.println("processMovingTargetsMulti(): FAILED TO READ TARGET DATA FROM "+disparity_path);
        			} else {
        				targets = disparity_targets;
        				System.out.println("processMovingTargetsMulti(): re-using target disparities from "+disparity_path);
        			}
        		}
        		//			if ((disparity_targets == null) || rng_glob) { // for now - always recalculate if 
        		if (disparity_targets == null) { 
        			if (debugLevel > -4) {
        				System.out.println("detectTargets(): Generating target disparities");
        			}
        			if (reset_disparity) {
        				for (int nseq = 0; nseq < targets.length; nseq++) {
        					for (int ntile = 0; ntile < targets[nseq].length; ntile++) if (targets[nseq][ntile] != null) {
        						targets[nseq][ntile][CuasMotionLMA.RSLT_DISPARITY] = cuas_initial_disparity + cuas_infinity; //  0;
        						targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_STR] = 0;
        					}
        				}
        			}
        			if (smooth_omegas) { // may be empty
        				if (debugLevel > -4) {
        					System.out.println ("Recalculating omegas to better fit target positions in consecutive key frames.");
        				}
        				cuasMotion.recalcOmegas(
        						true, // boolean       recalc,
        						targets,     // final double [][][] target_single_in,
        						batch_mode,  // boolean       batch_mode,
        						debugLevel); // final int             debugLevel){

        			} else {
        				cuasMotion.recalcOmegas(
        						false, // boolean       recalc,
        						targets,     // final double [][][] target_single_in,
        						batch_mode,  // boolean       batch_mode,
        						debugLevel); // final int             debugLevel){
        				if (debugLevel > -4) {
        					System.out.println ("Skipping omegas recalculation omegas, still calculating local target index (called GLOBAL_INDEX).");
        				}
        			}
        			//	String model_prefix = center_CLT.getImageName()+CuasMotion.getParametersSuffix(clt_parameters,null);
        			if (rng_glob) {
        				// glob_stats =  
        				rangeSingleTargets(
        						targets,          //final double [][][] targets,   // centers
        						log_ranging,      // final boolean       log_ranging,
        						rng_vfy,          // final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data)
        						debugLevel);      // final int           debugLevel)
        			}

        			if (disparity_targets == null) { // recalculate if was not read 
        				rangeTargets(
        						targets,          //final double [][][] targets,   // centers
        						//					rng_combine,      // final int           rng_combine,
        						log_ranging,      // final boolean       log_ranging,
        						rng_vfy,          // final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data)
        						debugLevel);      // final int           debugLevel)
        			}
        			if (!generate_csv) { // if csv is generated, data will be saved later
        				ImagePlus imp_with_range = CuasMotion.showTargetSequence(
        						targets,          // double [][][] vector_fields_sequence,
        						cuasMotion.getSliceTitles(),         // String []     titles, // all slices*frames titles or just slice titles or null
        						center_CLT.getImageName()+CuasMotion.getParametersSuffixRanging(clt_parameters,TARGET_DISPARITIES_SUFFIX),
        						!batch_mode,          // boolean       show,
        						cuasMotion.getTilesX());   //  int           tilesX) {
        				center_CLT.saveImagePlusInModelDirectory(imp_with_range);  // ImagePlus   imp)
        			}
        		} else if (recalc_target_id){// if (disparity_targets == null) { // at com.elphel.imagej.cuas.CuasMotion$46.run(CuasMotion.java:7269)
    				cuasMotion.recalcOmegas( // does not work - 
    						false, // boolean       recalc,
    						targets,     // final double [][][] target_single_in,
    						batch_mode,  // boolean       batch_mode,
    						debugLevel); // final int             debugLevel){
    				if (debugLevel > -4) {
    					System.out.println ("(Re)calculating local target index (called GLOBAL_INDEX). May be remolved later");
    				}
        		}
        	}
        	// calculate range from disparity
        	{ // always save after calculating ranges and adding UAS data. Later may be moved to after rangeTargets() above, only when calculated
        		getRangeFromDisparity(targets); // adds RSLT_INFINITY field
        		// generate results video (move from earlier)
        		addUasData(targets); // add flight log data to the nearest tiles, either existing or new
        		// re-saving data with flight log (ground truth) additions
        		ImagePlus imp_with_range = CuasMotion.showTargetSequence(
        				targets,          // double [][][] vector_fields_sequence,
        				cuasMotion.getSliceTitles(),         // String []     titles, // all slices*frames titles or just slice titles or null
        				center_CLT.getImageName()+CuasMotion.getParametersSuffixRanging(clt_parameters,TARGET_DISPARITIES_SUFFIX),
        				!batch_mode,          // boolean       show,
        				cuasMotion.getTilesX());   //  int           tilesX) {
        		center_CLT.saveImagePlusInModelDirectory(imp_with_range);  // ImagePlus   imp)
        	}
        	if (generate_csv) {
        		saveTargetStats(targets); // final double [][][] targets_single) {
        	}
        }
		
		
		if (generate_output) {
			cuasMotion.generateExtractFilterMovingTargets( // move parameters to clt_parameters
					radar_mode,              // final boolean       radar_mode,
					false,                   // final boolean       video_pass, // if clt_parameters.cuas_clean_video=true, video_pass=0 - output TIFFS, but no video. If video_pass==1 - only video with no debug
					batch_mode,              // final boolean       batch_mode,
					fpixels,                 // final float [][]    fpixels,
					targets, // _nonconflict,     // final double [][][] vf_sequence,  // center tiles (not extended), null /non-null only
					debugLevel);             // final int            debugLevel)
			if (clean_video) {
				cuasMotion.generateExtractFilterMovingTargets( // move parameters to clt_parameters
						radar_mode,              // final boolean       radar_mode,
						true, // final boolean       video_pass, // if clt_parameters.cuas_clean_video=true, video_pass=0 - output TIFFS, but no video. If video_pass==1 - only video with no debug
						batch_mode,              // final boolean       batch_mode,
						fpixels,                 // final float [][]    fpixels,
						targets, // _nonconflict,     // final double [][][] vf_sequence,  // center tiles (not extended), null /non-null only
						debugLevel);             // final int            debugLevel)
			}

		} else {
			System.out.println("Output images/videos with detected targets is disabled by \"Generate and save detected targets\" parameter, skipping it.");
			
		}

		
		
		if (rng_img) {
			double [][][] extended_targets = CuasMotion.extendMotionScan(
					targets,                 // final double [][][] motion_scan,
					null,                    // filter5,      // final boolean [][]  filtered, // centers, should be non-overlapped
					cuasMotion.getTilesX(),  // final int           tilesX)
					2,                       // final int           range, // 1 or 2
					null); // remain);              // final int     []    remain)
			if (debugLevel > -4) {
				System.out.println("detectTargets(): Generating and saving per-sensor target images");
			}
			double [][][] rendered_keyframes =  renderKeyFrames(
					targets,          //final double [][][] targets,   // centers
					extended_targets, // final double [][][] targets5x5,
					debugLevel);      // final int           debugLevel)
			boolean show = false; // true;
			int num_sens =   center_CLT.getNumSensors();
			String [] sens_titles = new String[num_sens+1];
			for (int i = 0; i < num_sens; i++) {
				sens_titles[i] = "SENS-"+i;
			}
			sens_titles[num_sens] = "COMBO";
			ImagePlus imp_rendered = ShowDoubleFloatArrays.showArraysHyperstack(
					rendered_keyframes,       // double[][][] pixels, 
					center_CLT.getWidth(),         // int          width, 
					center_CLT.getImageName()+"-RENDERED_TARGETS-R"+rng_radius0+":"+rng_radius+
					"-B"+rng_blur+"-UM"+rng_um_sigma,          // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					sens_titles,         // String []    titles, // all slices*frames titles or just slice titles or null
					cuasMotion.getSliceTitles(), // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
					show);          // boolean      show)
			center_CLT.saveImagePlusInModelDirectory(imp_rendered);            // ImagePlus   imp)
		}
		return cuasMotion;
	}
	
	public void getRangeFromDisparity(
			final double [][][] targets){
		//clt_parameters.imp.
		final double   max_disp_diff = clt_parameters.imp.cuas_max_disp_diff; //    0.05;  // Maximal disparity difference during last change to consider disparity valid
		final double   min_disp_str =  clt_parameters.imp.cuas_min_disp_str;   //    0.4;   // Minimal disparity strength to consider disparity valid
		final double   rng_limit =     clt_parameters.imp.cuas_rng_limit;         // 5000;     // Maximal displayed distance to target
		final double   cuas_infinity = clt_parameters.imp.cuas_infinity;
		final GeometryCorrection gc = center_CLT.getGeometryCorrection();
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nSeq = ai.getAndIncrement(); nSeq < targets.length; nSeq = ai.getAndIncrement()) {
						for (int ntile = 0; ntile < targets[nSeq].length; ntile++) if (targets[nSeq][ntile] != null){
							targets[nSeq][ntile][CuasMotionLMA.RSLT_INFINITY] = cuas_infinity;
							if (!Double.isNaN(targets[nSeq][ntile][CuasMotionLMA.RSLT_DISPARITY]) && (targets[nSeq][ntile][CuasMotionLMA.RSLT_DISPARITY] != 0.0)) {
								double disparity =      targets[nSeq][ntile][CuasMotionLMA.RSLT_DISPARITY] - cuas_infinity; // corrected disparity
								double disparity_diff = targets[nSeq][ntile][CuasMotionLMA.RSLT_DISP_DIFF];
								double strength =       targets[nSeq][ntile][CuasMotionLMA.RSLT_DISP_STR];
								if ((disparity_diff > max_disp_diff) || (strength < min_disp_str)) {
									targets[nSeq][ntile][CuasMotionLMA.RSLT_RANGE] = Double.NaN;
								} else if (disparity <= 0){
									targets[nSeq][ntile][CuasMotionLMA.RSLT_RANGE] = Double.POSITIVE_INFINITY;
								} else {
									double z = gc.getZFromDisparity(disparity); // or should we use sqrt(x^2 + y^2 + z^2)?
									if (z > rng_limit) {
										targets[nSeq][ntile][CuasMotionLMA.RSLT_RANGE] = Double.POSITIVE_INFINITY;
									} else {
										targets[nSeq][ntile][CuasMotionLMA.RSLT_RANGE] = Math.abs(z); // was it positive already?
									}
								}
							} else {
								targets[nSeq][ntile][CuasMotionLMA.RSLT_RANGE] = Double.NaN;
							}
							
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return;
	}
	
	
	
	public double [][][] renderKeyFrames(
			final double [][][] targets,   // centers
			final double [][][] targets5x5,
			final int           debugLevel) {
		int num_seq = cuasMotion.getNumCorrSamples();
		double [][][] render_targets = new double [num_seq][][];
		for (int nseq = 0; nseq < num_seq; nseq++) {
			render_targets[nseq] = renderKeyFrame(
					targets,     // final double [][][] targets,   // centers
					targets5x5,  // final double [][][] targets5x5,
					nseq);       // final int           nseq)
		}
		return render_targets;
	}
	
	/**
	 * Preparation for ranging - rendering per-pixel, around targets (now 5x5 or 3x3 tiles), separate for each sensor
	 * @param targets5x5 targets, expanded around each detected target (will not be needed for ranging)
	 * @param radius radius of the mask around the target for each sensor individually
	 * @param nseq   number of sequence (keyframe)
	 * @return rendered images for each sensor separately [sensor][pixel]
	 */
	public double [][] renderKeyFrame(
			final double [][][] targets,   // centers
			final double [][][] targets5x5,
			final int           nseq) {
		final boolean um_en =          clt_parameters.imp.cuas_rng_um;
		final boolean um_all =         clt_parameters.imp.cuas_rng_um_all;
		final double um_sigma =        clt_parameters.imp.cuas_rng_um_sigma;
		final boolean um_twice =       clt_parameters.imp.cuas_rng_um2;
		final double um_weight =       clt_parameters.imp.cuas_rng_um_weight;
		final double radius0 =         clt_parameters.imp.cuas_rng_radius0;
		final boolean mb_en =          clt_parameters.imp.mb_en; //  && (fov_tiles==null) && (mode3d > 0);
		final double  mb_tau =         clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
		final 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
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);

		final int frame_center =     cuasMotion.getFrameCenter(nseq);
		final int half_accum_range = cuasMotion.getSeqLength()/2;
		//		final boolean   unsharped 

		final double [] window_full= cuasMotion.getSegmentWindow(
				clt_parameters.imp.cuas_rng_coswnd, // boolean smooth)
				true); // boolean normalize)
		//radius_blur
		final int start_scene = frame_center - half_accum_range;
		final int num_scenes = 2*half_accum_range + 1;
		boolean   clean_up = true;  // remove unneded data
		final double [][][] img_um_seq = getUMSequence(
				start_scene, // final int       start_scene,
				num_scenes,   // final int       num_scenes,
				clean_up,     // final boolean   clean_up, // remove unneded data
				false,        // final boolean   erase_source_images, // erase source images (will not clean up this.img_um)
				debugLevel);  // final int debugLevel)
		/*
		final double [][][] img_um_seq = QuadCLT.unsharpMaskSourceMono(
				scenes,      // final QuadCLT[] scenes,
				start_scene, // final int       start_scene,
				num_scenes,  // final int       num_scenes,
				um_en,       // final boolean   um_en,
				um_all,      // final boolean   unsharped,
				um_sigma,    // final double    um_sigma,
				um_twice,    // final boolean   um_twice,
				um_weight,   // final double    um_weight,
				debugLevel); // final int       debugLevel)
				*/
		double [][][] pXpYD5x5s = cuasMotion.targetPxPyD(
				targets5x5[nseq]); // final double [][] targets)
		double [][][] pXpYDs = cuasMotion.targetPxPyD(
				targets[nseq]); // final double [][] targets)

		///   	boolean batch_run =clt_parameters.batch_run; // may be modified for debug
		///		boolean cuas_debug =  clt_parameters.imp.cuas_debug;  // save debug images (and show them if not in batch mode)
		final int tilesX =     center_CLT.getTilesX();
		final int tilesY =     center_CLT.getTilesY();
		final int tileSize =   center_CLT.getTileSize();
		final int sensor_mask = -1; // all sensors
		final int num_sens =   center_CLT.getNumSensors();
		final int width = tilesX * tileSize;
		final int height = tilesY * tileSize;
		ErsCorrection ers_reference = center_CLT.getErsCorrection();
		int num_used_sens = 0;
		for (int i = 0; i < num_sens; i++) if (((sensor_mask >> i) & 1) != 0) num_used_sens++;
		int [] channels = new int [num_used_sens];
		int nch = 0;
		for (int i = 0; i < num_sens; i++) if (((sensor_mask >> i) & 1) != 0) channels[nch++] = i;
		double [][] rendered = new double [num_sens + 1][width*height];
		for (int dseq = half_accum_range; dseq <= half_accum_range; dseq++) {
			final double [][] ref_pXpYD5x5 = pXpYD5x5s[dseq + half_accum_range];
			final double [][] ref_pXpYD =    pXpYDs[dseq + half_accum_range]; // center target
			
			final int nscene = frame_center + dseq;
			final int sm = -1; // merge_all? -1: sensor_mask;
			double [][] drender;
			final double [][] dxyzatr_dt = (mb_en ? new double[][] { // for all, including ref
				scenes[nscene].getErsCorrection().getErsXYZ_dt(),
				scenes[nscene].getErsCorrection().getErsATR_dt()} : null);
			final String ts = scenes[nscene].getImageName();
			final double [] scene_xyz = ers_reference.getSceneXYZ(ts);
			final double [] scene_atr = ers_reference.getSceneATR(ts);
			if ((scene_atr==null) || (scene_xyz == null)) {  // should not happen
				System.out.println("renderKeyFrame() BUG1");
				continue;
			}
		    ImageDtt image_dtt = new ImageDtt(
		    		center_CLT.getNumSensors(),
		    		clt_parameters.transform_size,
		    		clt_parameters.img_dtt,
		    		center_CLT.isAux(),
		    		center_CLT.isMonochrome(),
		    		center_CLT.isLwir(),
		    		clt_parameters.getScaleStrength(center_CLT.isAux()),
		    		center_CLT.getGPU());
			prepareAltImages(
					image_dtt,        // final ImageDtt      image_dtt,
					ref_pXpYD,        // final double [][]   ref_pXpYD,
					img_um_seq,       // final double [][][] img_um_seq,
					scene_xyz,        // final double []     scene_xyz,       
					scene_atr,        // final double []     scene_atr,
			        radius0,          // final double        radius,
					nseq,             // final int           nseq,
					-half_accum_range,// final int           first_offs,
					dseq,             // final int           dseq,
					debugLevel);      // final int           debugLevel)

			
			if (mb_en && (dxyzatr_dt != null)) {
				double [][] motion_blur = OpticalFlow.getMotionBlur(
						center_CLT, // quadCLTs[ref_index],   // QuadCLT        ref_scene,
						scenes[nscene],        // QuadCLT        scene,         // can be the same as ref_scene
						ref_pXpYD5x5,             // 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)
				// TODO: Add target movement (recalculate velocities) 

				drender = QuadCLT.renderDoubleGPUFromDSI(
						sm,                  // final int         sensor_mask,
						false, // merge_all,           // final boolean     merge_channels,
						null,                // final Rectangle   full_woi_in,      // show larger than sensor WOI (or null)
						clt_parameters,      // CLTParameters     clt_parameters,
						null,                // ref_disparity,       // double []         disparity_ref,
						ref_pXpYD5x5,           // 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
						scenes[nscene],    // final QuadCLT     scene,
						center_CLT, // quadCLTs[ref_index], // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
						false,               // final boolean     toRGB,
						clt_parameters.imp.show_mono_nan,
						debugLevel);         // int         debugLevel)
			} else {
				drender = QuadCLT.renderDoubleGPUFromDSI(
						sm,                  // final int         sensor_mask,
						false, // merge_all,           // final boolean     merge_channels,
						null, // fov_tiles,           // testr, // null,                // final Rectangle   full_woi_in,      // show larger than sensor WOI (or null)
						clt_parameters,      // CLTParameters     clt_parameters,
						null, // ref_disparity,       // double []         disparity_ref,
						ref_pXpYD5x5,   // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
						0.0,              // double            mb_tau,      // 0.008; // time constant, sec
						0.0,         // double            mb_max_gain, // 5.0;   // motion blur maximal gain (if more - move second point more than a pixel
						null,         // double [][]       mb_vectors,  //
						// not used, just as null/not null now
						// null means uniform grid, no view transform. even with 0 rot ERS was changing results
						scene_xyz,           // final double []   scene_xyz, // camera center in world coordinates
						scene_atr,           // final double []   scene_atr, // camera orientation relative to world frame
						scenes[nscene],    // final QuadCLT     scene,
						center_CLT, // quadCLTs[ref_index], // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
						false,               // final boolean     toRGB,
						clt_parameters.imp.show_mono_nan,
						debugLevel);         // int         debugLevel)
			}
			ai.set(0);
			final double scale = window_full[half_accum_range + dseq];
			for (int ithread = 0; ithread < threads.length; ithread++) {
				threads[ithread] = new Thread() {
					public void run() {
						for (int nPixY = ai.getAndIncrement(); nPixY < height; nPixY = ai.getAndIncrement()) {
							for (int npixX = 0; npixX < width; npixX++) {
								int npix = npixX + width* nPixY;
								if (!Double.isNaN(rendered[0][npix])) {
									for (int nchn = 0; nchn < num_sens; nchn++) {
										rendered[nchn][npix] += scale * drender[nchn][npix];
									}
								}
							}
						}
					}
				};
			}		      
			ImageDtt.startAndJoin(threads);
		}
		Arrays.fill(rendered[num_sens], Double.NaN);
		ai.set(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nPixY = ai.getAndIncrement(); nPixY < height; nPixY = ai.getAndIncrement()) {
						for (int npixX = 0; npixX < width; npixX++) {
							int npix = npixX + width* nPixY;
							if (!Double.isNaN(rendered[0][npix])) {
								double d = 0;
								for (int nchn = 0; nchn < num_sens; nchn++) {
									d += rendered[nchn][npix];
								}
								rendered[num_sens][npix] = d/num_sens;
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return rendered;
	}
	
	public void rangeSingleTargets(
			final double [][][] targets,   // centers
			final boolean       log_ranging,
			final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data) NOT implemented
			final int           debugLevel) {
		int [][][] linked_targets = CuasMotion.getLinkedTargets(targets);
///		double [][] target_glob_stats = new double [linked_targets.length][];
		int    dbg_rng_tgt = clt_parameters.imp.cuas_dbg_rng_tgt;
		int    show_mode =   clt_parameters.imp.cuas_dbg_show_mode; //  2; // 7; // 0; //  4; // 0; //3; +1 - disparity, +2 - corr
		
		float [][][]  accum_2d_corr = null;
		if (log_ranging) {
			String log_head = getRangingLogParameters();
			log_head +="target\tnseq\tframe_center\tnrefine\ttileX\ttileY\ttile\tdisparity\tdisp_diff\thalf\tstrength\tlast_str\tmax_diff\n";
			center_CLT.appendStringInModelDirectory(log_head, TARGET_RANGING_LOGS_SUFFIX);
			center_CLT.appendStringInModelDirectory(log_head, TARGET_RANGING_LOGS_FULL_SUFFIX);
		}
		int target_id0 = 1;
		int target_id1 = linked_targets.length;
		if (dbg_rng_tgt >0) {
			target_id0  = dbg_rng_tgt;
			target_id1  = dbg_rng_tgt;
		}
		System.out.println ("rangeSingleTargets(): Processing target_id from "+target_id0+" to "+target_id1+" (inclusive)"); // for debugging
		for (int target_id = target_id0; target_id <= target_id1; target_id++) {
			boolean rng_vfy_this = (target_id == target_id0) && rng_vfy; // only for the first target in range
///			double [] target_stats = 
			rangeSingleTarget( // return {disparity, strength, last disparity difference}
					targets,        // final double [][][] targets,   // centers
					linked_targets, // final int [][][]    linked_targets,
					target_id,      // final int           target_id, // target number, starting with 1
					log_ranging,    // final boolean       log_ranging,
					accum_2d_corr,  // final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
					show_mode,      // final int           show_mode, //+1 - disparity, +2 - correlations +4 - debug single run
					rng_vfy_this,   // final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data) NOT implemented
					debugLevel) ;   // final int           debugLevel) ;
///			target_glob_stats[target_id-1] = target_stats;
		}
		return; //  target_glob_stats;
	}
	
	
	
	public void rangeTargets(
			final double [][][] targets,   // centers
//			final int           rng_combine,
			final boolean       log_ranging,
			final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data)
			final int           debugLevel) {
		int num_seq =      cuasMotion.getNumCorrSamples();
		int show_mode =    clt_parameters.imp.cuas_dbg_show_mode; //  2; // 7; // 0; //  4; // 0; //3; +1 - disparity, +2 - corr
		int rng_combine =  clt_parameters.imp.cuas_rng_combine; // // combine multiple scenes before intrascene correlation
		int dbg_rng_seq =  clt_parameters.imp.cuas_dbg_rng_seq;
		float [][][]  accum_2d_corr = null;
		int nseq_first = 0;
		int nseq_last = num_seq-1;
		if (dbg_rng_seq >= 0) {
			nseq_first = dbg_rng_seq;
			nseq_last =  dbg_rng_seq;
		}
		System.out.println("rangeTargets(): debug/change nseq_first="+nseq_first+", nseq_last="+nseq_last);
		// for rng_vfy debugging
		final int half_accum_range = cuasMotion.getSeqLength()/2;
		final int num_sens =   center_CLT.getNumSensors();
		final int num_scenes =   2*half_accum_range+1;
		final int num_scene_grp = (num_scenes + (rng_combine - 1)) / rng_combine;
//		final double [][][] ranging_verify = (rng_vfy)  ? (new double [num_scenes + 1][num_sens+1][]) : null;
		final double [][][] ranging_verify = (rng_vfy)  ? (new double [num_scene_grp + 1][num_sens + 1][]) : null;
		
		if (log_ranging) {
			String log_head = getRangingLogParameters();
			log_head +="target\tnseq\tframe_center\tnrefine\ttileX\ttileY\ttile\tdisparity\tdisp_diff\thalf\tstrength\tlast_str\tmax_diff\n";
			center_CLT.appendStringInModelDirectory(log_head, TARGET_RANGING_LOGS_SUFFIX);
			center_CLT.appendStringInModelDirectory(log_head, TARGET_RANGING_LOGS_FULL_SUFFIX);
		}
		
		for (int nseq = nseq_first; nseq <= nseq_last; nseq++) {
			int show_mode1 = show_mode & ((nseq == nseq_first)? 7 : 3);
			double [][][] ranging_verify_use = (nseq == nseq_first) ? ranging_verify : null ;
			rangeTargets(
					 targets,            // final double [][][] targets,   // centers
					 log_ranging,        // final boolean       log_ranging,
					 accum_2d_corr,      // final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
					 nseq,               // final int           nseq,
					 show_mode1,         // final int           show_mode, //+1 - disparity, +2 - correlations
					 ranging_verify_use, // final double [][][] ranging_verify,
					 debugLevel);        // final int           debugLevel)
		}
		return;
	}

	String getRangingLogParameters() {
    	StringBuffer sb = new StringBuffer();
    	sb.append("\n------------------------\n");
    	sb.append(new SimpleDateFormat("yyyy/MM/dd HH:mm:ss").format(Calendar.getInstance().getTime())+"\n");
    	sb.append("cuas_smooth_omegas = "   + clt_parameters.imp.cuas_smooth_omegas	  + "\n"); // Recalculate omegas from continuing targets positions 
    	sb.append("cuas_rng_um = "    	    + clt_parameters.imp.cuas_rng_um          + "\n"); // apply UM to images for ranging (assumed true)
    	sb.append("cuas_rng_um_all = "      + clt_parameters.imp.cuas_rng_um_all      + "\n"); // Unsharp mask all images (after FPN and Row/Col), unchecked - when needed. Checked faster, but original data is lost. 
    	sb.append("cuas_rng_um_sigma = "    + clt_parameters.imp.cuas_rng_um_sigma    + "\n"); // unsharp mask sigma for ranging images
    	sb.append("cuas_rng_um2 = "         + clt_parameters.imp.cuas_rng_um2         + "\n"); //  apply UM to images twice
    	sb.append("cuas_rng_um_weight = "   + clt_parameters.imp.cuas_rng_um_weight   + "\n"); // unsharp mask weight for ranging images
    	sb.append("cuas_rng_coswnd = "      + clt_parameters.imp.cuas_rng_coswnd      + "\n"); // Use cosine scenes window (false - rectangular)
    	sb.append("cuas_rng_combine = "     + clt_parameters.imp.cuas_rng_combine     + "\n"); // combine multiple scenes before intrascene correlation
    	sb.append("cuas_rng_radius0 = "     + clt_parameters.imp.cuas_rng_radius0     + "\n"); // mask out data outside peripheral areas, keep 0.5 at the radius - first iteration 
    	sb.append("cuas_rng_radius = "      + clt_parameters.imp.cuas_rng_radius      + "\n"); // mask out data outside peripheral areas, keep 0.5 at the radius  - next iteration
    	sb.append("cuas_rng_blur = "        + clt_parameters.imp.cuas_rng_blur        + "\n"); // relative transition radius range (for 0.5 the mask inside (1-0.5) radius will be 1.0, outside (1+0.5) radius will be 0
    	sb.append("cuas_rng_fz0 = "         + clt_parameters.imp.cuas_rng_fz0         + "\n"); // Fat zero in target ranging mode (first iteration)
    	sb.append("cuas_rng_fz = "          + clt_parameters.imp.cuas_rng_fz          + "\n"); // Fat zero in target ranging mode (next iterations)
    	sb.append("cuas_rng_scale = "       + clt_parameters.imp.cuas_rng_scale       + "\n"); // Scale alt_data to use same strength as for normal ranging because targets
    	sb.append("cuas_lma_bypass = "      + clt_parameters.imp.cuas_lma_bypass      + "\n"); //  Bypass solution tests, OK weak 
    	sb.append("cuas_mcorr_sel = "       + clt_parameters.imp.cuas_mcorr_sel       + "\n"); // all pairs and diagonals (as was for non-cuas ranging)
    	sb.append("cuas_mcorr_sel_lma = "   + clt_parameters.imp.cuas_mcorr_sel_lma   + "\n"); // all pairs (will use all selected by cuas_mcorr_sel)
    	sb.append("cuas_max_rel_rms = "     + clt_parameters.imp.cuas_max_rel_rms     + "\n"); // maximal relative (to average max/min amplitude LMA RMS) // May be up to 0.3)
    	sb.append("cuas_min_strength = "    + clt_parameters.imp.cuas_min_strength    + "\n"); // minimal composite strength (sqrt(average amp squared over absolute RMS)
    	sb.append("cuas_ac_offset = "       + clt_parameters.imp.cuas_ac_offset       + "\n"); // add to a,c coefficients for near-lines where A,C could become negative because of window 
    	sb.append("cuas_min_max_ac = "      + clt_parameters.imp.cuas_min_max_ac      + "\n"); // 0.175; // 0.14;  // LWIR16: 0.01  maximal of a and C coefficients minimum (measures sharpest point/line)
    	sb.append("cuas_min_min_ac = "      + clt_parameters.imp.cuas_min_min_ac      + "\n"); // LWIR16: 0.007 minimal of a and C coefficients minimum (measures sharpest point)
    	sb.append("cuas_reset_disparity = " + clt_parameters.imp.cuas_reset_disparity + "\n"); // reset target disparities from infinity
    	sb.append("cuas_initial_disparity= "+ clt_parameters.imp.cuas_initial_disparity+ "\n");// Start correlation with this disparity (in addition to infinity) after reset
    	sb.append("cuas_infinity = "        + clt_parameters.imp.cuas_infinity        + "\n"); // disparity at infinity for targets
    	sb.append("cuas_rng_niterate = "    + clt_parameters.imp.cuas_rng_niterate    + "\n"); // number of disparity iterations
    	sb.append("cuas_rng_diff = "        + clt_parameters.imp.cuas_rng_diff        + "\n"); // exit when disparity difference is smaller
    	
    	sb.append("cuas_max_disp_diff = "   + clt_parameters.imp.cuas_max_disp_diff   + "\n"); // Maximal disparity difference during last change to consider disparity valid
    	sb.append("cuas_min_disp_str = "    + clt_parameters.imp.cuas_min_disp_str    + "\n"); // Minimal disparity strength to consider disparity valid
    	sb.append("cuas_rng_limit = "       + clt_parameters.imp.cuas_rng_limit       + "\n"); // Maximal displayed distance to target
    	sb.append("cuas_infinity = "        + clt_parameters.imp.cuas_infinity        + "\n"); // Disparity at infinity
    	
    	sb.append("cuas_flw_levels = "      + clt_parameters.imp.cuas_flw_levels      + "\n"); // 1 - all, 2 - all and two halves, 3 - all, two halves and 4 quaters
    	sb.append("cuas_flw_len = "         + clt_parameters.imp.cuas_flw_len         + "\n"); // Minimal number of key frames in each segment after division;
    	sb.append("cuas_flw_diff = "        + clt_parameters.imp.cuas_flw_diff        + "\n"); // (pix) Minimal difference between halves to use "tilted" non-constant disparity 

    	return sb.toString();
	}
	
	public void rangeSingleTarget( // return {disparity, strength, last disparity difference}
			final double [][][] targets,   // centers
			final int [][][]    linked_targets,
			final int           target_id, // target number, starting with 1
			final boolean       log_ranging,
			final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
			final int           show_mode, //+1 - disparity, +2 - correlations +4 - debug single run
			final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data) NOT implemented
			final int           debugLevel) {
		int    flw_levels =    clt_parameters.imp.cuas_flw_levels;
		int    flw_len =       clt_parameters.imp.cuas_flw_len;
///		double flw_diff =      clt_parameters.imp.cuas_flw_diff;
		double max_disp_diff = clt_parameters.imp.cuas_max_disp_diff; //    0.05;  // Maximal disparity difference during last change to consider disparity valid
		double min_disp_str =  clt_parameters.imp.cuas_min_disp_str;   //    0.4;   // Minimal disparity strength to consider disparity valid
		double rng_limit =     clt_parameters.imp.cuas_rng_limit;         // 5000;     // Maximal displayed distance to target
		double cuas_infinity = clt_parameters.imp.cuas_infinity;
		
		// extract target from the first to the last defined sequence number
		int very_first_seq = -1;
		for (int i = 0; i < targets.length; i++) {
			if (linked_targets[target_id-1][i] != null) {
				very_first_seq = i;
				break;
			}
		}
		if (very_first_seq < 0) {
			System.out.println("rangeSingleTarget(): target "+target_id+" is undefined");
			return; //  null;
		}
		int iseq;
		for (iseq = very_first_seq; (iseq < targets.length) && (linked_targets[target_id-1][iseq] != null); iseq++ );
		int [] target_tiles_full = new int [iseq - very_first_seq];
		for (int i = 0; i < target_tiles_full.length; i++) target_tiles_full[i] = linked_targets[target_id-1][very_first_seq+i][0];
		final GeometryCorrection gc = center_CLT.getGeometryCorrection();
		process_segments: {
			double [] seg_rslt = null;
			for (int seg_lev = 0; seg_lev < flw_levels; seg_lev++) {
				int num_segm = (1 << seg_lev);
				for (int segm_num = 0; segm_num < num_segm; segm_num++) {
					int seg_start = (target_tiles_full.length * segm_num) /num_segm;
					int seg_end = ((target_tiles_full.length * (segm_num + 1)) /num_segm) - 1;
					int [] target_tiles = new int [seg_end-seg_start+1];
					System.arraycopy(target_tiles_full, seg_start, target_tiles, 0, target_tiles.length);
					seg_rslt = rangeSingleTarget( // return {disparity, strength, last disparity difference}
							targets,       // final double [][][] targets,   // centers
							target_tiles,  // final int []        target_tiles,
							seg_start+very_first_seq,     // final int           first_seq, // 
							target_id,     // final int           target_id, // target number, starting with 1
							log_ranging,   // final boolean       log_ranging,
							accum_2d_corr, // final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
							show_mode,     // final int           show_mode, //+1 - disparity, +2 - correlations +4 - debug single run
							rng_vfy,       // final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data) NOT implemented
							debugLevel);   // final int           debugLevel)
					if (seg_rslt == null) {
						System.out.println("rangeSingleTarget():failed to calculate disparity for target #"+target_id+
								", segment level="+seg_lev+", segment number="+segm_num+", first_seq="+very_first_seq+", num_seq="+target_tiles.length);
					} else {
						int seg_length = seg_end - seg_start +1;
						int middle_index = seg_start + seg_length/2; // middle for odd, next after middle for even
						int ntile = target_tiles[middle_index - seg_start];
						double [] target = targets[middle_index + very_first_seq][ntile];
						double disparity =      seg_rslt[0] - cuas_infinity;
						double strength =       seg_rslt[1];
						double disparity_diff = seg_rslt[2];
						double range =          Double.NaN;
						if ((disparity_diff > max_disp_diff) || (strength < min_disp_str)) {
							disparity = Double.NaN;
						} else if (disparity <= 0){
							range = Double.POSITIVE_INFINITY;
						} else {
							double z = gc.getZFromDisparity(disparity); // or should we use sqrt(x^2 + y^2 + z^2)?
							if (z > rng_limit) {
								range = Double.POSITIVE_INFINITY;
							} else {
								range = z; // was it positive already? Yes
							}
						}
						target[CuasMotionLMA.RSLT_GLENGTH] =    seg_length; // null pointer
						target[CuasMotionLMA.RSLT_GDISPARITY] = seg_rslt[0]; // disparity; here - raw disparity, not corrected
						target[CuasMotionLMA.RSLT_GDISP_DIFF] = disparity_diff;
						target[CuasMotionLMA.RSLT_GDISP_STR] =  strength;
						target[CuasMotionLMA.RSLT_GRANGE] =     range;
					}
				}
				// verify next leve's lengths
				if ((target_tiles_full.length >> (seg_lev + 1)) < flw_len) {
					break process_segments;
				}
			}
		}
        return; // null;
	}
	
	
	
	public double[] rangeSingleTarget( // return {disparity, strength, last disparity difference}
			final double [][][] targets,   // centers
			final int []        target_tiles,
			final int           first_seq, // 
			final int           target_id, // target number, starting with 1
			final boolean       log_ranging,
			final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
			final int           show_mode, //+1 - disparity, +2 - correlations +4 - debug single run
			final boolean       rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data) NOT implemented
			final int           debugLevel) {
		int refine_half_alt = 2;    // when refine >=  refine_half_alt and diff changed sign - multiply it by 0.5
		final int end_mode = clt_parameters.imp.cuas_glob_ends; // 1; // 0 - same as internal, 1 - cosine extended ends, 2 - rectangular extended ends // make a parameter
		final boolean first_last_extra = (end_mode > 0); 
		final boolean smooth_ends = (end_mode == 1);
		/*
		// extract target from the first to the last defined sequence number
		int first_seq = -1;
		for (int i = 0; i < targets.length; i++) {
			if (linked_targets[target_id-1][i] != null) {
				first_seq = i;
				break;
			}
		}
		if (first_seq < 0) {
			System.out.println("rangeSingleTarget(): target "+target_id+" is undefined");
			return null;
		}
		int iseq;
		for (iseq = first_seq; (iseq < targets.length) && (linked_targets[target_id-1][iseq] != null); iseq++ );
		int [] target_tiles = new int [iseq - first_seq];
		for (int i = 0; i < target_tiles.length; i++) target_tiles[i] = linked_targets[target_id-1][first_seq+i][0];
		*/
		
		boolean show_disparity=     (show_mode & 1) != 0;
		boolean show_corr=          (show_mode & 2) != 0;
		boolean show_debug_single = (show_mode & 4) != 0;
		
		final int rng_niterate =  clt_parameters.imp.cuas_rng_niterate;
		final double rng_diff =   clt_parameters.imp.cuas_rng_diff;
		double cuas_initial_disparity=clt_parameters.imp.cuas_initial_disparity; // 1.0;  // Start correlation with this disparity (in addition to infinity) after reset
		double cuas_infinity =      clt_parameters.imp.cuas_infinity ; // 0.63;  // disparity at infinity for targets
		
		final int tilesX =     center_CLT.getTilesX();
//		final int tilesY =     center_CLT.getTilesY();
//		final int frame_center =     cuasMotion.getFrameCenter(nseq); // for debug only
		final boolean use_non_lma = false;
		
		// TODO: move to parameters!
		final double worsen_str = 0.3; // 9; // do not update if new strength is below worsen_str * old_str
		// Show disparity_map;
		StringBuffer sb =      null;
		StringBuffer sb_full = null;
		StringBuffer sb_last = null;
    	if (log_ranging) {
    		sb =     new StringBuffer();
    		sb_full = new StringBuffer();
    	}
		int refine_vfy = 0; // on which refine step to generate verify images
		// debugging
		final boolean center_only = debugLevel < 1000; // > 1000; // < 1000;
		boolean show_dbg = false; // true;
		
		double disparity =  cuas_initial_disparity + cuas_infinity; //  0;
		double strength =   0;
		double disparity_diff = Double.NaN;
		int ref_tile = target_tiles[0]; // moving targets will have accumulated correlations in this tile
		int ref_tileX = ref_tile % tilesX; 
		int ref_tileY = ref_tile / tilesX; 
		for (int nrefine = 0; nrefine < rng_niterate; nrefine++) {
			boolean rng_vfy_this = rng_vfy && (nrefine == refine_vfy);
			double [][] disparity_map = refineSingleTargetDisparity(
					targets,          // final double [][][] targets,   // centers
					accum_2d_corr,    //final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
					first_seq,        // final int    first_seq,
					target_tiles,     // final int [] target_tiles,
					disparity,        // final double disparity0,
					nrefine,          // final int    nrefine, // number of the refine to select radius (larger first time)
					first_last_extra, // final boolean first_last_extra,
					smooth_ends,      // final boolean smooth_ends,
					show_corr,        // final boolean show_corr,
					rng_vfy_this,     // final boolean rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data)
					debugLevel);      // final int           debugLevel)
			// render this.ranging_verify here, regardless of the disparity_map
			if (rng_vfy_this) {
				ImagePlus imp_vfy = renderRangingVerifySingle(
						target_id,  // final int           target_id,
						nrefine,    // final int           nrefine,
						false);     // final boolean       show)
				center_CLT.saveImagePlusInModelDirectory(imp_vfy);
			}
			
			if (disparity_map == null) {
				System.out.println("rangeTargets()): disparity_map==null on nrefine="+nrefine);
				break;
			}
			double disp_diff = disparity_map[ImageDtt.DISPARITY_INDEX_POLY][ref_tile]; // null
			double str = disparity_map[ImageDtt.DISPARITY_INDEX_POLY+1][ref_tile];
			if (Double.isNaN(disp_diff)) {
				if (use_non_lma) {
					disp_diff = disparity_map[ImageDtt.DISPARITY_INDEX_CM][ref_tile];
					str = disparity_map[ImageDtt.DISPARITY_INDEX_CM+1][ref_tile];
				}
				if (Double.isNaN(disp_diff)) {
					System.out.println ("rangeSingleTarget(): Failed to range target #"+target_id+" for sequence "+ first_seq+" ("+ref_tileX+":"+ref_tileY+").");
					return null;
				}
			}
			
			String shalf = " ";
			if ((nrefine >= refine_half_alt) && ((disparity_diff * disp_diff) < 0)) {  // disp_diff changed sign
				disp_diff *= 0.5;
				shalf = "*";
			}
			
			boolean improved = false;
			if ((str > 0) && !((worsen_str * strength) > str)) {
				disparity += disp_diff;
				disparity_diff =  disp_diff;
				strength = str;
				improved = true;
			}
			double abs_diff = Math.abs(disp_diff);
			if ((sb != null) || (sb_full != null)) {
				sb_last = new StringBuffer();
			}
			if (sb_last != null) {
				
				sb_last.append(target_id+"\t"+first_seq+"\t"+target_tiles.length+"\t"+nrefine+"\t"+ref_tileX+"\t"+ref_tileY+"\t"+ref_tile+"\t"+
						disparity+"\t"+
						disp_diff+"\t"+shalf+"\t"+
						strength+"\t"+
						str+"\t"+
						abs_diff+"\n");
			}
			
			if (debugLevel > -4) {
				System.out.println(String.format("rangeSingleTarget(), target=%2d, first_seq=%3d (%3d), nrefine=%2d, %3d:%3d (%4d), disparity=%7.4f (%7.4f %s), strength=%7.5f (%7.5f) abs_diff=%7.5f",
						target_id, first_seq,target_tiles.length,nrefine,ref_tileX,ref_tileY,ref_tile,
						disparity,
						disp_diff, shalf,
						strength,str,abs_diff));
			}
			
			if (sb_full != null) {
				sb_full.append(sb_last);
//				sb_full.append("\n");
			}
			if (!improved) {
				System.out.println ("rangeSingleTarget(): No improvement for target #"+target_id+" for sequence "+ first_seq+" ("+ref_tileX+":"+ref_tileY+").");
				break;
			}
			if (abs_diff <= rng_diff) {
				break;
			}
		} //		for (int nrefine = 0; nrefine < rng_niterate; nrefine++) {
		// if never was set - make it NaN
		if (strength <= 0) { // never was set
			disparity = Double.NaN;
			strength = 0;
		}
		if (sb != null) {
			sb.append(sb_last);
			center_CLT.appendStringInModelDirectory(sb.toString(), TARGET_RANGING_LOGS_SUFFIX);
		}
		if (sb_full != null) {
			center_CLT.appendStringInModelDirectory(sb_full.toString(), TARGET_RANGING_LOGS_FULL_SUFFIX);
		}
		return new double [] {disparity, strength, disparity_diff};
	}
	
	

	public void rangeTargets(
			final double [][][] targets,   // centers
			final boolean       log_ranging,
			final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
			final int           nseq,
			final int           show_mode, //+1 - disparity, +2 - correlations +4 - debug single run
			final double [][][] ranging_verify,
			final int           debugLevel) {
		boolean show_disparity=     (show_mode & 1) != 0;
		boolean show_corr=          (show_mode & 2) != 0;
		boolean show_debug_single = (show_mode & 4) != 0;
		int refine_half_alt = 2;    // when refine >=  refine_half_alt and diff changed sign - multiply it by 0.5
		final int rng_niterate =  clt_parameters.imp.cuas_rng_niterate;
		final double rng_diff =   clt_parameters.imp.cuas_rng_diff;
		double cuas_initial_disparity=clt_parameters.imp.cuas_initial_disparity; // 1.0;  // Start correlation with this disparity (in addition to infinity) after reset
		double cuas_infinity =      clt_parameters.imp.cuas_infinity ; // 0.63;  // disparity at infinity for targets
		
		final int tilesX =     center_CLT.getTilesX();
		final int tilesY =     center_CLT.getTilesY();
		final int frame_center =     cuasMotion.getFrameCenter(nseq); // for debug only
		final boolean use_non_lma = false;
		
		// TODO: move to parameters!
		final double worsen_str = 0.3; // 9; // do not update if new strength is below worsen_str * old_str
		// Show disparity_map;
		StringBuffer sb =      null;
		StringBuffer sb_full = null;
		StringBuffer sb_last = null;
    	if (log_ranging) {
    		sb =     new StringBuffer();
    		sb_full = new StringBuffer();
    	}
		int refine_vfy = 0; // on which refine step to generate verify images
		// debugging
		final boolean center_only = debugLevel < 1000; // > 1000; // < 1000;
		boolean show_dbg = false; // true;
		for (int nrefine = 0; nrefine < rng_niterate; nrefine++) {
			boolean reset_disp = debugLevel > 1000;
			if (reset_disp) {
				for (int dbg_nseq = 0; dbg_nseq < targets.length; dbg_nseq++) {
					for (int dbg_ntile = 0; dbg_ntile < targets[nseq].length; dbg_ntile++) if (targets[dbg_nseq][dbg_ntile] != null) {
						targets[dbg_nseq][dbg_ntile][CuasMotionLMA.RSLT_DISPARITY] = cuas_initial_disparity + cuas_infinity; //  0;;
						targets[dbg_nseq][dbg_ntile][CuasMotionLMA.RSLT_DISP_STR] = 0;
					}
				}
				String dbg_title = center_CLT.getImageName()+
						"TARGETS-RESET-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma +  
						"-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
						"-F"+ (cuasMotion.getFrameCenter(nseq))+":"+nrefine+"-M"+clt_parameters.imp.cuas_mcorr_sel+
						"-"+clt_parameters.imp.cuas_mcorr_sel_lma;
				CuasMotion.showTargetSequence(
						targets,        // double [][][] vector_fields_sequence,
						cuasMotion.getSliceTitles(),         // String []     titles, // all slices*frames titles or just slice titles or null
						dbg_title,// String        title,
						true,          // boolean       show,
						cuasMotion.getTilesX());   //  int           tilesX) {
			}
			
			
			boolean vfy_now = true; // (nrefine == refine_vfy);
			double [][][] ranging_verify_this = (vfy_now? ranging_verify : null);
			
//			double [][][] extended_targets = null;
			double [][][] var_targets = targets;
			
			if (!center_only && (ranging_verify_this != null)) {
				var_targets = CuasMotion.extendMotionScan(
						targets,                 // final double [][][] motion_scan,
						null,                    // filter5,      // final boolean [][]  filtered, // centers, should be non-overlapped
						cuasMotion.getTilesX(),  // final int           tilesX)
						2,                       // final int           range, // 1 or 2
						null); // remain);              // final int     []    remain)
			}
			if (show_dbg) {
				String dbg_title = center_CLT.getImageName()+
						"VAR_TARGETS-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma +  
						"-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
						"-F"+ (cuasMotion.getFrameCenter(nseq))+":"+nrefine+"-M"+clt_parameters.imp.cuas_mcorr_sel+
						"-"+clt_parameters.imp.cuas_mcorr_sel_lma;
				CuasMotion.showTargetSequence(
						var_targets,        // double [][][] vector_fields_sequence,
						cuasMotion.getSliceTitles(),         // String []     titles, // all slices*frames titles or just slice titles or null
						dbg_title,// String        title,
						true,          // boolean       show,
						cuasMotion.getTilesX());   //  int           tilesX) {
			}
			

			double [][] disparity_map = refineTargetDisparity( // returns disparity_map
					var_targets,         // targets,             // final double [][][] targets,   // centers
					accum_2d_corr,       //final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
					ranging_verify_this, // final double [][][] ranging_verify,
					nseq,                // final int           nseq,
					nrefine,             // final int           nrefine, // number of the refine (just for debug)
					show_corr,           // final boolean       show_corr,
					show_debug_single && (nrefine == 0), // final boolean       show_debug_single,
					debugLevel);         // final int           debugLevel)
			if (disparity_map == null) {
				System.out.println("rangeTargets()): disparity_map==null on nrefine="+nrefine);
				break;
			}
			if (ranging_verify_this != null) {
				int half_accum_range = cuasMotion.getSeqLength()/2;
				ImagePlus imp_vfy = renderRangingVerify(
						ranging_verify_this,
						frame_center-half_accum_range, // final int nscene0);
						nrefine,                       // final int           nrefine,
						false);                        // final boolean       show)
				center_CLT.saveImagePlusInModelDirectory(imp_vfy);
			}
			
			if (show_disparity) {
				ShowDoubleFloatArrays.showArrays(
						disparity_map,
						tilesX,
						tilesY,
						true,
						center_CLT.getImageName()+"-FZ"+clt_parameters.imp.cuas_rng_fz0+":"+clt_parameters.imp.cuas_rng_fz+
						"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma +  
						"-accumulated_disparity_map_R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
						"-F"+ (cuasMotion.getFrameCenter(nseq))+":"+nrefine+"-M"+clt_parameters.imp.cuas_mcorr_sel+
						"-"+clt_parameters.imp.cuas_mcorr_sel_lma,
						ImageDtt.getDisparityTitles(center_CLT.getNumSensors(),center_CLT.isMonochrome()) // ImageDtt.DISPARITY_TITLES
						);
			}
			if ((sb != null) || (sb_full != null)) {
				sb_last = new StringBuffer();
			}

			double max_diff = 0.0;
			for (int ntile = 0; ntile < targets[nseq].length; ntile++) if ((targets[nseq][ntile] != null) && 
//					(disparity_map != null) &&
					(disparity_map[ImageDtt.DISPARITY_INDEX_POLY] != null)) {
				double disp_diff = disparity_map[ImageDtt.DISPARITY_INDEX_POLY][ntile]; // null
				double str = disparity_map[ImageDtt.DISPARITY_INDEX_POLY+1][ntile];
				int itarget = (int) targets[nseq][ntile][CuasMotionLMA.RSLT_GLOBAL];

				if (Double.isNaN(disp_diff)) {
					if (!use_non_lma) {
						continue;
					}
					disp_diff = disparity_map[ImageDtt.DISPARITY_INDEX_CM][ntile]; // /0.85;
					str = disparity_map[ImageDtt.DISPARITY_INDEX_CM+1][ntile]; // /0.85;
					if (Double.isNaN(disp_diff)) {
						disp_diff = 0;
						str=0;
						continue;
					}
				}
				if (Double.isNaN(targets[nseq][ntile][CuasMotionLMA.RSLT_DISPARITY])) {
					targets[nseq][ntile][CuasMotionLMA.RSLT_DISPARITY] = 0.0;
				}
				String shalf = " ";
				if ((nrefine >= refine_half_alt) && ((targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_DIFF] * disp_diff) < 0)) { // disp_diff changed sign
					disp_diff *= 0.5;
					shalf = "*";
				}
				if ((str > 0) && !((worsen_str * targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_STR]) > str)) {
					targets[nseq][ntile][CuasMotionLMA.RSLT_DISPARITY] += disp_diff;
					targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_DIFF] =  disp_diff;
					targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_STR] = str;
				} else {
					disp_diff = 0;
				}
				max_diff = Math.max(max_diff, Math.abs(disp_diff));
				if (sb_last != null) {
					int dbg_tilex = ntile % tilesX;
					int dbg_tiley = ntile / tilesX;
					
					sb_last.append(itarget+"\t"+nseq+"\t"+frame_center+"\t"+nrefine+"\t"+dbg_tilex+"\t"+dbg_tiley+"\t"+ntile+"\t"+
							targets[nseq][ntile][CuasMotionLMA.RSLT_DISPARITY]+"\t"+
							disp_diff+"\t"+
							shalf+"\t"+
							targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_STR]+"\t"+
							str+"\t"+
							max_diff+"\n");
				}
				
				if (debugLevel > -4) {
					int dbg_tilex = ntile % tilesX;
					int dbg_tiley = ntile / tilesX;
					System.out.println(String.format("rangeTargets(), target=%2d, nseq=%3d (%3d), nrefine=%2d, %3d:%3d (%4d), disparity=%7.4f (%7.4f %s), strength=%7.5f (%7.5f) max_diff=%7.5f",
							itarget, nseq,frame_center,nrefine,dbg_tilex,dbg_tiley,ntile,
							targets[nseq][ntile][CuasMotionLMA.RSLT_DISPARITY],
							disp_diff,shalf,
							targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_STR],str,max_diff));
					
				}
				
				if (sb_full != null) {
					sb_full.append(sb_last);
//					sb_full.append("\n");
				}
			}
			if (max_diff <= rng_diff) {
				break;
			}
		}
		// if never was set - make it NaN
		for (int ntile = 0; ntile < targets[nseq].length; ntile++) if ((targets[nseq][ntile] != null) && !(targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_STR] > 0)) {
			targets[nseq][ntile][CuasMotionLMA.RSLT_DISP_STR] = 0;
			targets[nseq][ntile][CuasMotionLMA.RSLT_DISPARITY] = Double.NaN;
		}

		
		if (sb != null) {
			sb.append(sb_last);
			center_CLT.appendStringInModelDirectory(sb.toString(), TARGET_RANGING_LOGS_SUFFIX);
		}
		if (sb_full != null) {
			center_CLT.appendStringInModelDirectory(sb_full.toString(), TARGET_RANGING_LOGS_FULL_SUFFIX);
		}
		return;
	}
	
	public ImagePlus renderRangingVerify(
			final double [][][] ranging_verify, // normalized, each scene group value is divided by number of averaged scenes
			final int           nscene0,
			final int           nrefine,
			final boolean       show) {
		if (ranging_verify == null) {
			return null;
		}

		final int half_accum_range = cuasMotion.getSeqLength()/2;
		final int num_scenes =   2*half_accum_range+1;
		int rng_combine =            clt_parameters.imp.cuas_rng_combine; // // combine multiple scenes before intrascene correlation
		final int num_scene_grp = (num_scenes + (rng_combine - 1)) / rng_combine; //==ranging_verify.length - 1;

		//		final int num_scenes =       ranging_verify.length-1;
		final int num_sens =         ranging_verify[0].length-1;
		final int num_pix =          ranging_verify[0][0].length;
//		final int half_accum_range = (num_scenes-1)/2;

		String [] titles_scenes = new String [num_scene_grp + 1];
		for (int igrp = 0; igrp < (titles_scenes.length-1); igrp++) {
			titles_scenes[igrp] = cuasMotion.getSceneTitles()[nscene0+igrp * rng_combine];
		}
		titles_scenes[titles_scenes.length-1]="COMBO";
		String [] titles_sensors = new String [num_sens+1];
		for (int i = 0; i < (titles_sensors.length-1); i++) {
			titles_sensors[i] = "SENS_"+i;
		}
		titles_sensors[titles_sensors.length-1]="ALL";
		// calculate averages: per scene (all sensors)
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nSceneGrp = ai.getAndIncrement(); nSceneGrp < num_scene_grp; nSceneGrp = ai.getAndIncrement()) {
//						int iscene =  nSceneGrp * rng_combine;
//						int grp_len = Math.min(rng_combine, num_scenes-iscene);
						ranging_verify[nSceneGrp][num_sens] = new double [num_pix];
						for (int npix = 0; npix < num_pix; npix++) {
							double d = 0;
							for (int nsens = 0; nsens < num_sens; nsens++) {
								d += ranging_verify[nSceneGrp][nsens][npix];
							}
							ranging_verify[nSceneGrp][num_sens][npix] = d/num_sens;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		ai.set(0);
		// calculate averages: per senSors (including "all")
		final double [] window_full= cuasMotion.getSegmentWindow( // normalized, sum == 1.0
				clt_parameters.imp.cuas_rng_coswnd, // boolean smooth)
				false); // boolean normalize)
		final double [] window_grp = new double [num_scene_grp];
		double sw = 0;
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			double w = window_full[nscene];
			window_grp[nscene / rng_combine] += w;
			sw += w;
		}
		final double sum_window = sw;

		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nSens = ai.getAndIncrement(); nSens < (num_sens+1); nSens = ai.getAndIncrement()) {
						ranging_verify[num_scene_grp][nSens] = new double [num_pix];
						for (int npix = 0; npix < num_pix; npix++) {
							double d = 0;
							for (int nscene_grp = 0; nscene_grp < num_scene_grp; nscene_grp++) {
								d += ranging_verify[nscene_grp][nSens][npix]; // sum before normalizing
								ranging_verify[nscene_grp][nSens][npix] /= window_grp[nscene_grp]; // normalize each group image 
							}
							ranging_verify[num_scene_grp][nSens][npix] = d / sum_window;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		
		ImagePlus imp_vfy = ShowDoubleFloatArrays.showArraysHyperstack(
				ranging_verify,       // double[][][] pixels, 
				center_CLT.getWidth(),         // int          width, 
				center_CLT.getImageName()+"-RANGING_VERIFY-R"+"-FZ"+clt_parameters.imp.cuas_rng_fz0+":"+clt_parameters.imp.cuas_rng_fz+
				"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma +  
	            "-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
	            "-F"+ (nscene0 + half_accum_range)+":"+nrefine+"-M"+clt_parameters.imp.cuas_mcorr_sel+
	            "-"+clt_parameters.imp.cuas_mcorr_sel_lma,
	            titles_sensors,         // String []    titles, // all slices*frames titles or just slice titles or null
	            titles_scenes, // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
				show);          // boolean      show)
		return imp_vfy;
	}

	public ImagePlus renderRangingVerifySingle(
			final int           target_id,
			final int           nrefine,
			final boolean       show) {
		if (ranging_verify == null) {
			return null;
		}
		String [] scene_titles = cuasMotion.getSceneTitles();
		final int num_scenes = scene_titles.length;
		ArrayList<Integer> indices_list = new ArrayList<Integer>();
		for (int i = 0; i < num_scenes; i++) if (this.ranging_verify[i] != null) {
			indices_list.add(i);
		}
		final int[] indices = indices_list.stream().mapToInt(Integer::intValue).toArray();
		final int num_used_scenes = indices.length;
		String [] titles_scenes =      new String [num_used_scenes + 1];
		final double [][][] ranging_verify = new double [indices.length + 1][][];
		final double []     ranging_verify_weight = new double [indices.length + 1];
		int sw = 0;
		for (int i = 0; i < indices.length; i++) {
			ranging_verify[i] =        this.ranging_verify[indices[i]];
			ranging_verify_weight[i] = this.ranging_verify_weight[indices[i]];
			sw += ranging_verify_weight[i];
			titles_scenes[i] = scene_titles[indices[i]];
		}
		titles_scenes[num_used_scenes]="COMBO";
		final int num_sens =                     ranging_verify[0].length-1;
		ranging_verify[num_used_scenes] =        new double [num_sens+1][]; // ranging_verify[0][0].length];
		ranging_verify_weight[num_used_scenes] = sw;

		String [] titles_sensors = new String [num_sens+1];
		for (int i = 0; i < (titles_sensors.length-1); i++) {
			titles_sensors[i] = "SENS_"+i;
		}
		titles_sensors[titles_sensors.length-1]="ALL";
///		final int half_accum_range = cuasMotion.getSeqLength()/2;
///		final int num_scenes =   2*half_accum_range+1;
		int rng_combine =            clt_parameters.imp.cuas_rng_combine; // // combine multiple scenes before intrascene correlation
///		final int num_scene_grp = (num_scenes + (rng_combine - 1)) / rng_combine; //==ranging_verify.length - 1;

///		final int num_scenes =       ranging_verify.length-1;
		final int num_pix =          ranging_verify[0][0].length;
///		final int half_accum_range = (num_scenes-1)/2;

		// calculate averages: per scene (all sensors)
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int iScene = ai.getAndIncrement(); iScene < num_used_scenes; iScene = ai.getAndIncrement()) {
						ranging_verify[iScene][num_sens] = new double [num_pix];
						for (int npix = 0; npix < num_pix; npix++) {
							double d = 0;
							for (int nsens = 0; nsens < num_sens; nsens++) {
								d += ranging_verify[iScene][nsens][npix];
							}
							ranging_verify[iScene][num_sens][npix] = d/num_sens;
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		ai.set(0);
		// calculate averages: per sensors (including "all")

		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nSens = ai.getAndIncrement(); nSens < (num_sens+1); nSens = ai.getAndIncrement()) {
						ranging_verify[num_used_scenes][nSens] = new double [num_pix];
						for (int npix = 0; npix < num_pix; npix++) {
							double d = 0;
							for (int iscene = 0; iscene < num_used_scenes; iscene++) {
								d += ranging_verify[iscene][nSens][npix]; // sum before normalizing
								ranging_verify[iscene][nSens][npix] /= ranging_verify_weight[iscene]; // normalize each group image 
							}
							ranging_verify[num_used_scenes][nSens][npix] = d / ranging_verify_weight[num_used_scenes];
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		
		ImagePlus imp_vfy = ShowDoubleFloatArrays.showArraysHyperstack(
				ranging_verify,       // double[][][] pixels, 
				center_CLT.getWidth(),         // int          width, 
				center_CLT.getImageName()+"-RANGING_VERIFY-R"+"-FZ"+clt_parameters.imp.cuas_rng_fz0+":"+clt_parameters.imp.cuas_rng_fz+
				"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma +  
	            "-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
	            "-ID"+ target_id+":"+nrefine+"-M"+clt_parameters.imp.cuas_mcorr_sel+
	            "-"+clt_parameters.imp.cuas_mcorr_sel_lma,
	            titles_sensors,         // String []    titles, // all slices*frames titles or just slice titles or null
	            titles_scenes, // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
				show);          // boolean      show)
		return imp_vfy;
	}
	
	
	
	
	public double [][][] getUMSequence(
			final int       start_scene,
			final int       num_scenes,
			final boolean   clean_up, // remove unneded data
			final boolean   erase_source, // erase source images (will not clean up this.img_um, will erase used source images instead)
			final int debugLevel) {
		final boolean um_en =       clt_parameters.imp.cuas_rng_um;
		final double  um_sigma =    clt_parameters.imp.cuas_rng_um_sigma;
		final boolean um_twice =    clt_parameters.imp.cuas_rng_um2;
		final double  um_weight =   clt_parameters.imp.cuas_rng_um_weight;
		if (clean_up) {
			for (int nscene = 0; nscene < start_scene; nscene++) {
				if (!erase_source) { // keep cache if source images are deleted
					this.img_um[nscene] = null;
				}
				scenes[nscene].setImageDataAlt(null);
			}
			for (int nscene = start_scene+num_scenes; nscene < this.img_um.length; nscene++) {
				if (!erase_source) { // keep cache if source images are deleted
					this.img_um[nscene] = null;
				}
				scenes[nscene].setImageDataAlt(null);
			}
			Runtime runtime = Runtime.getRuntime();
		    runtime.gc();
			if (debugLevel > -3) System.out.println("getUMSequence(): Free memory="+runtime.freeMemory()+" (of "+runtime.totalMemory()+")");
		}
		
		final double [][][] img_um_seq = 
				QuadCLT.unsharpMaskSourceMono(
						scenes,      // final QuadCLT[] scenes,
						start_scene, // final int       start_scene,
						num_scenes,  // final int       num_scenes,
			    		um_en,       // final boolean   um_en,
			    		img_um,      // final double [][][] img_um, // cache for unsharped images, same length as scenes. May be null
			    		false,       // final boolean   unsharped, // ==false, not used 
			    		false,       // final boolean   update_source, // ==false, not used if update_source will not generate output, return null
			    		erase_source,// final boolean   erase_source,
			    		um_sigma,    // final double    um_sigma,
			    		um_twice,    // final boolean   um_twice,
			    		um_weight,   // final double    um_weight,
			    		debugLevel); // final int       debugLevel)
		return img_um_seq;
	}
	
	
	
	/**
	 * Calculate disparity for a sequence of keyframes following a single target that may change tiles along the way
	 * @param targets      targets array with a single target per tile and marked target numbers in CuasMotionLMA.RSLT_GLOBAL field
	 * @param first_seq    first sequence (key frame) number (0<= first_seq < targets.length)
	 * @param target_tiles array of target indices for consecutive key frames
	 * @param disparity0   initial disparity to be refined 
	 * @param nrefine      number of refine cycle (0-th uses larger radius) 
	 * @param first_last_extra add longer scene sequence before the first and after the last keyframe (same as fo per-keyframe ranging, apply windowing)  
	 * @param debugLevel   debug level
	 * @return disparity   map (only one tile will be used - the first one with the selected target
	 */
	public double [][] refineSingleTargetDisparity(
			final double [][][] targets,   // centers
			final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
			final int    first_seq,
			final int [] target_tiles,
			final double disparity0,
			final int    nrefine, // number of the refine to select radius (larger first time)
			final boolean first_last_extra,
			final boolean smooth_ends,
			final boolean show_corr,
			final boolean rng_vfy, // Generate/save ranging verification images (per-sensor and combined rendering from the same data)
			final int    debugLevel){
		final boolean show_sensors_vfy = true; // generate (unmasked) sensor images when rng_vfy is true;
		final boolean show_alt_vfy = true; // generate (unmasked) sensor images when rng_vfy is true;
//		final boolean show_corr = true;
//		final boolean first_last_extra = true; // use longer sequences for the first/last keyframe
//		final int num_tiles = targets[0].length;
//		final double  um_sigma =    clt_parameters.imp.cuas_rng_um_sigma; // only for file names
		final double  rng_fz =      (nrefine == 0)?clt_parameters.imp.cuas_rng_fz0 : clt_parameters.imp.cuas_rng_fz;
		final boolean mb_en =       clt_parameters.imp.mb_en; //  && (fov_tiles==null) && (mode3d > 0);
		final double  mb_tau =      clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
		final 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
		final double gpu_sigma_corr =     clt_parameters.getGpuCorrSigma(center_CLT.isMonochrome());
		final double gpu_sigma_rb_corr =  center_CLT.isMonochrome()? 1.0 : clt_parameters.gpu_sigma_rb_corr;
		final double gpu_sigma_log_corr = clt_parameters.getGpuCorrLoGSigma(center_CLT.isMonochrome());
		final double radius =       (nrefine == 0)? clt_parameters.imp.cuas_rng_radius0: clt_parameters.imp.cuas_rng_radius;
		final int rng_combine = clt_parameters.imp.cuas_rng_combine; // // combine multiple scenes before intrascene correlation
		final int half_accum_range = cuasMotion.getSeqLength()/2;
		boolean   clean_up = true;  // remove unneeded data
		boolean   erase_source = true; // check everything works
		final int frame_center0 =     cuasMotion.getFrameCenter(first_seq);
		final int frame_step =        cuasMotion.getCorrInc(); // 20
		final int hframe_step = frame_step/2;
		final int start_scene = frame_center0 - (first_last_extra ? half_accum_range: hframe_step);
		final int num_seq = target_tiles.length;
//		final int num_scenes = first_last_extra ? ((num_seq -1) * frame_step + 2* half_accum_range + 1) : (num_seq * frame_step);
		final int num_scenes = ((num_seq-1) * frame_step) + 2 * (first_last_extra ? half_accum_range : (frame_step/ 2)) + 1;

		
		final double [] window_nnorm = new double [num_scenes];
		Arrays.fill(window_nnorm, 1.0);
		if (rng_vfy) { // create for each scene, unused will be nulls
			ranging_verify = new double [cuasMotion.getSceneTitles().length + 1][][];
			ranging_verify_weight = new double [ranging_verify.length];
		}
		if (first_last_extra ) {
			final double [] window1= cuasMotion.getSegmentWindow( // normalized, sum == 1.0
					smooth_ends, // boolean smooth)
					false);      // boolean normalize) // not normalized, values up to 1.0
			for (int i = 0; i <= half_accum_range; i++) {
				window_nnorm[i] = window1[i];
				window_nnorm[num_scenes - 1 - i] = window1[i];
			}
		}
		double sum_wnd = 0;
		for (int i = 0; i < window_nnorm.length; i++) {
			sum_wnd+=window_nnorm[i];
		}
		
		
		
		final double [][][] img_um_seq = getUMSequence(
				start_scene,  // final int       start_scene,
				num_scenes,   // final int       num_scenes,
				clean_up,     // final boolean   clean_up, // remove unneded data
				erase_source, // final boolean   erase_source, // erase source images (will not clean up this.img_um)
				debugLevel);  // final int debugLevel)
		if (rng_vfy && show_sensors_vfy) {
			String [] titles_sensors = new String[img_um_seq[0].length];
			for (int i = 0; i < titles_sensors.length; i++) titles_sensors[i] = "SENS-"+i;
			String [] titles_scenes = new String[img_um_seq.length];
			for (int i = 0; i < titles_scenes.length; i++) {
				titles_scenes[i] = cuasMotion.getSceneTitles()[start_scene+i];
			}
			boolean   show = false; // true;
			ImagePlus imp_sensors_um =ShowDoubleFloatArrays.showArraysHyperstack(
					img_um_seq,       // double[][][] pixels, 
					center_CLT.getWidth(),         // int          width, 
					center_CLT.getImageName()+"-SENSOR_FULL_IMAGES"+
					//"-R"+"-FZ"+clt_parameters.imp.cuas_rng_fz+
					//"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+
					"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma  
		            //"-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
		            //"-ID"+ target_id+":"+nrefine+
		            //"-M"+clt_parameters.imp.cuas_mcorr_sel+
		            //"-"+clt_parameters.imp.cuas_mcorr_sel_lma
		            ,
		            titles_sensors,         // String []    titles, // all slices*frames titles or just slice titles or null
		            titles_scenes, // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
					show);          // boolean      show)
        	center_CLT.saveImagePlusInModelDirectory(
        			null,            // "GPU-SHIFTED-D"+clt_parameters.disparity, // String      suffix,
        			imp_sensors_um); // imp_scenes); // ImagePlus   imp)
		}
		
		
/*
		double [][][] pXpYDs = cuasMotion.targetPxPyD(
				targets[nseq]); // final double [][] targets)
*/
		///   	boolean batch_run =clt_parameters.batch_run; // may be modified for debug
		///		boolean cuas_debug =  clt_parameters.imp.cuas_debug;  // save debug images (and show them if not in batch mode)
		final int tilesX =     center_CLT.getTilesX();
		final int tilesY =     center_CLT.getTilesY();
		final int sensor_mask = -1; // all sensors
		final int num_sens =   center_CLT.getNumSensors();
		
		ImageDttParameters imgdtt_params = null;
		try {
			imgdtt_params = clt_parameters.img_dtt.clone();
		} catch (CloneNotSupportedException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		imgdtt_params.updateFromCuas(clt_parameters.imp, num_sens);
		
		final int mcorr_sel = ImageDttParameters.corrSelEncode(imgdtt_params,num_sens);
		ErsCorrection ers_reference = center_CLT.getErsCorrection();
		int num_used_sens = 0;
		for (int i = 0; i < num_sens; i++) if (((sensor_mask >> i) & 1) != 0) num_used_sens++;
		int [] channels = new int [num_used_sens];
		int nch = 0;
		for (int i = 0; i < num_sens; i++) if (((sensor_mask >> i) & 1) != 0) channels[nch++] = i;

		// from OF13477: 	public static double[][] correlateInterscene(
		final int num_pairs = Correlation2d.getNumPairs(num_sens);
		final float  [][][][]   fcorr_td_acc  = new float [tilesY][tilesX][][];
		final float  [][][]     accum_weights = new float [tilesY][tilesX][num_pairs];
		final boolean show_2d_corr = false;
		boolean show_accumulated_correlations = show_2d_corr || debugLevel > -5;
		//		boolean show_reference_correlations =  show_2d_corr || debugLevel > -5;
		//		final float  [][][]       fclt_corr = ((accum_2d_corr != null) || show_accumulated_correlations || show_reference_correlations) ?
		final float  [][][]       fclt_corr = ((accum_2d_corr != null) || show_accumulated_correlations) ?
				(new float [tilesX * tilesY][][]) : null;
		final boolean        use_rms = true; // DISPARITY_STRENGTH_INDEX means LMA RMS (18/04/2023) (from OF11168:		boolean        use_rms = true; )
		final boolean        no_map = false;

		ImageDtt image_dtt = new ImageDtt(
				center_CLT.getNumSensors(),
				clt_parameters.transform_size,
				imgdtt_params, // clt_parameters.img_dtt,
				center_CLT.isAux(),
				center_CLT.isMonochrome(),
				center_CLT.isLwir(),
				clt_parameters.getScaleStrength(center_CLT.isAux()),
				center_CLT.getGPU());
		image_dtt.getCorrelation2d(); // initiate image_dtt.correlation2d, needed if disparity_map != null  
		final double[][] disparity_map = no_map ? null : new double [image_dtt.getDisparityTitles().length][];

		TpTask[][] tp_tasks = null; // will contain last tp_tasks
		TpTask[][] tp_tasks_ref = null; // copy from the first segment dseq==0
		int ntile0 = target_tiles[0];
		int tileX0 = ntile0 % tilesX;
		int tileY0 = ntile0 / tilesX;
		
		// initialize fclt and normalize window
		for (int rnseq = 0; rnseq < num_seq; rnseq++) { // relative sequence number, starting from 0
			int nseq = first_seq + rnseq; // absolute sequence number
			int frame_center = cuasMotion.getFrameCenter(nseq);
//			int first_offs = ((rnseq==0) && first_last_extra)?             -half_accum_range  : -hframe_step;
//			int last_offs =  ((rnseq==(num_seq - 1)) && first_last_extra)?  half_accum_range  : (frame_step - hframe_step);
			int first_offs = -hframe_step;
			if ((rnseq==0) && first_last_extra) {
				first_offs = -half_accum_range;
			}
			int last_offs = frame_step - hframe_step - 1; // -1 for inclusive
			if (rnseq==(num_seq - 1)) {
				if (first_last_extra) {
					last_offs = half_accum_range;
				} else {
					last_offs = frame_step - hframe_step; // 1 longer
				}
			}
			
//			int rscene0 = frame_center + first_offs; // negative offset 
//			int rscene1 = frame_center + last_offs;  // positive offset
			int ntile = target_tiles[rnseq];
			int tileX = ntile % tilesX;
			int tileY = ntile / tilesX;
			double [][][] pXpYDs = cuasMotion.targetPxPyDSingleTile(
					targets[nseq], // final double [][] targets,
					ntile,          // final int         ntile,
					disparity0,     // final double      disparity,
					first_offs,     // final int         first_offs,
					last_offs);     // final int         last_offs) {
			// split segment by integration groups
			double wgrp = 0;
			for (int dseq = first_offs; dseq <= last_offs; dseq++) {
//				final int iscene = dseq + half_accum_range; // 0, 1,..., 2* half_accum_range.
				final int iscene = dseq - first_offs; // starting from 0;
				final double [][] ref_pXpYD =    pXpYDs[iscene]; // center target
				final int nscene = frame_center + dseq;  // absolute scene number
				final int iscene_grp = iscene / rng_combine;
				final boolean first_in_grp = (iscene % rng_combine) == 0;
				final boolean last_in_grp =  (((iscene + 1) % rng_combine) == 0) || (dseq == last_offs);
				if (first_in_grp) {
					wgrp = 0;
				}
				wgrp += window_nnorm[iscene];

				final int sm = -1; // merge_all? -1: sensor_mask;
				final double [][] dxyzatr_dt = (mb_en ? new double[][] { // for all, including ref
					scenes[nscene].getErsCorrection().getErsXYZ_dt(),
					scenes[nscene].getErsCorrection().getErsATR_dt()} : null);
				final String ts = scenes[nscene].getImageName();
				final double [] scene_xyz = ers_reference.getSceneXYZ(ts);
				final double [] scene_atr = ers_reference.getSceneATR(ts);
				if ((scene_atr==null) || (scene_xyz == null)) {  // should not happen
					System.out.println("renderKeyFrame() BUG1");
					continue;
				}
				int seq_offs = first_offs - rnseq * frame_step; // start of he img_um_seq[] relative to the center of he interval 
				prepareAltImages(
						image_dtt,   //  final ImageDtt      image_dtt,
						ref_pXpYD,   // final double [][]   ref_pXpYD,
						img_um_seq,  // final double [][][] img_um_seq,
						scene_xyz,   // final double []     scene_xyz,       
						scene_atr,   // final double []     scene_atr,
						radius,      // final double        radius,
						seq_offs, // first_offs,  // final int           first_offs,
						nseq,        // final int           nseq,
						dseq,        // final int           dseq,
						debugLevel); // final int           debugLevel)
				tp_tasks = new TpTask[(mb_en && (dxyzatr_dt != null)) ? 2 : 1][];
				if (mb_en && (dxyzatr_dt != null)) {
					double [][] motion_blur = OpticalFlow.getMotionBlur(
							center_CLT, // quadCLTs[ref_index],   // QuadCLT        ref_scene,
							scenes[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)
					// TODO: Add target movement (recalculate velocities) 

					QuadCLT.preRenderGPUFromDSI(
							!first_in_grp,    // final boolean     accumulate,
							window_nnorm[iscene],// final double      global_scale, // <=1.0
							sm,               // final int         sensor_mask,
							false,            // final boolean     merge_channels,
							0,                // final int         discard_border,
							0,                // final double      max_fold,
							0,                // final int         min_in_row_col,   // Minimal number of defined tiles in a row/column
							null,             // final Rectangle   full_woi_in,      // show larger than sensor WOI in tiles (or null)
							clt_parameters,   // CLTParameters     clt_parameters,
							null,             // double []         disparity_ref,
							// only center tiles here!
							ref_pXpYD,        // 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,  // now [2][ntiles];
							scene_xyz,        // double []         scene_xyz, // camera center in world coordinates. If null - no shift, no ers
							scene_atr,        // double []         scene_atr, // camera orientation relative to world frame
							scenes[nscene],   //final QuadCLT     scene,
							center_CLT,       // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
							tp_tasks, 			  // TpTask [][]       tasks,
							false, // clt_parameters.imp.show_mono_nan,         // final boolean     show_nan,
							ImageDtt.THREADS_MAX, // threadsMax,       // int               threadsMax,
							debugLevel);      // final int         debugLevel)
				} else {
					QuadCLT.preRenderGPUFromDSI(
							!first_in_grp,    // final  boolean     accumulate,
							window_nnorm[iscene], // final double      global_scale, // <=1.0
							sm,               // final int         sensor_mask,
							false,            // final boolean     merge_channels,
							0,                // final int         discard_border,
							0,                // final double      max_fold,
							0,                // final int         min_in_row_col,   // Minimal number of defined tiles in a row/column
							null,             // final Rectangle   full_woi_in,      // show larger than sensor WOI in tiles (or null)
							clt_parameters,   // CLTParameters     clt_parameters,
							null,             // double []         disparity_ref,
							// only center tiles here!
							ref_pXpYD,        // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
							// motion blur compensation 
							0.0,              // double            mb_tau,      // 0.008; // time constant, sec
							0.0,              // 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];
							scene_xyz,        // double []         scene_xyz, // camera center in world coordinates. If null - no shift, no ers
							scene_atr,        // double []         scene_atr, // camera orientation relative to world frame
							scenes[nscene],   //final QuadCLT     scene,
							center_CLT,       // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
							tp_tasks, 		  // TpTask [][]       tasks,
							false, // clt_parameters.imp.show_mono_nan,         // final boolean     show_nan,
							ImageDtt.THREADS_MAX, // threadsMax,       // int               threadsMax,
							debugLevel);      // final int         debugLevel)
				}
				// copy center tasks from the first series to the reference
				if ((rnseq == 0) && (dseq == 0)) {
					tp_tasks_ref = tp_tasks; // for the middle frame
					for (int i = 0; i < tp_tasks.length; i++) { // here just one
						tp_tasks_ref[i] = tp_tasks[i].clone();
					}
				}
				// correlate only after the last in group, others - just accumulate in TD
				if (last_in_grp) {
					if (rng_vfy) { // this should not be in a thread
						if (ranging_verify == null) { // save to this.ranging_verify sparce array
							ranging_verify = new double [cuasMotion.getSceneTitles().length + 1][][];
							ranging_verify_weight = new double [ranging_verify.length];
						}
						if (ranging_verify[nscene] == null) {
							ranging_verify[nscene]= new double [num_sens+1][];
						}
						double [][] scene_render = center_CLT.renderDoubleFromTDMono ( // [sensor][pixel]
								-1,     // int                 sensor_mask,
								null,   // int []              wh, // may be null, or {width, height}
								false); // boolean             use_reference
						for (int nsens = 0; nsens < scene_render.length; nsens++) {
							ranging_verify[nscene][nsens] = scene_render[nsens];
						}
						ranging_verify_weight[nscene] = wgrp;
					}

					float  [][][][]     fcorr_td =       new float[tilesY][tilesX][][];
					image_dtt.quadCorrTD(
							imgdtt_params, // clt_parameters.img_dtt,            // final ImageDttParameters imgdtt_params,    // Now just extra correlation parameters, later will include, most others
							tp_tasks[0], // *** will be updated inside from GPU-calculated geometry
							fcorr_td, // fcorrs_td[nscene],                 // [tilesY][tilesX][pair][4*64] transform domain representation of 6 corr pairs
							clt_parameters.gpu_sigma_r,        // 0.9, 1.1
							clt_parameters.gpu_sigma_b,        // 0.9, 1.1
							clt_parameters.gpu_sigma_g,        // 0.6, 0.7
							clt_parameters.gpu_sigma_m,        //  =       0.4; // 0.7;
							gpu_sigma_rb_corr,                 // final double              gpu_sigma_rb_corr, //  = 0.5; // apply LPF after accumulating R and B correlation before G, monochrome ? 1.0 : gpu_sigma_rb_corr;
							gpu_sigma_corr,                    //  =    0.9;gpu_sigma_corr_m
							gpu_sigma_log_corr,                // final double              gpu_sigma_log_corr,   // hpf to reduce dynamic range for correlations
							clt_parameters.corr_red,           // +used
							clt_parameters.corr_blue,          // +used
							mcorr_sel,                         // final int                 mcorr_sel,    // Which pairs to correlate // +1 - all, +2 - dia, +4 - sq, +8 - neibs, +16 - hor + 32 - vert
							ImageDtt.THREADS_MAX,       // maximal number of threads to launch
							debugLevel);
					if (image_dtt.getGPU().getGpu_debug_level() > -1) {
						System.out.println("==ooo=after image_dtt.quadCorrTD()");
					}
					// Verify tasks are now updated
					/*
					OpticalFlow.accumulateCorrelations(
							//						window_full[iscene] * (2 * half_accum_range + 1),        // final double         weight, // keep same sum as before
							wgrp, // window_grp[iscene_grp],        // final double         weight, // keep same sum as before
							accum_weights, // final int [][][]     num_acc,          // number of accumulated tiles [tilesY][tilesX][pair]
							fcorr_td,      // final float [][][][] fcorr_td,         // [tilesY][tilesX][pair][256] sparse transform domain representation of corr pairs 
							fcorr_td_acc); // final float [][][][] fcorr_td_acc      // [tilesY][tilesX][pair][256] sparse transform domain representation of corr pairs
					*/
					if (fcorr_td[tileY][tileX] != null) {
						if (fcorr_td_acc[tileY0][tileX0] == null) {
							fcorr_td_acc[tileY0][tileX0] = new float [fcorr_td[tileY][tileX].length][];
						}
						OpticalFlow.accumulateCorrelations(
								wgrp, // double     weight,
								accum_weights[tileY0][tileX0], // float []   num_acc,          // number of accumulated tiles [pair]
								fcorr_td[tileY][tileX], // float [][] fcorr_td,         // [pair][256] sparse transform domain representation of corr pairs 
								fcorr_td_acc[tileY0][tileX0]); // float [][] fcorr_td_acc      // [pair][256] sparse transform domain representation of corr pairs // should not be null, same length as fcorr_td
					}
					
					if (image_dtt.getGPU().getGpu_debug_level() > -1) {
						System.out.println("==ooo=accumulateCorrelations()");
					}
				} // if (last_in_grp) {

			} // for (int dseq = first_offs; dseq <= last_offs; dseq++) {
			
		} // for (int nseq = 0; nseq < num_seq; nseq++) {
		
		// Normalize accumulated correlations
		OpticalFlow.accumulateCorrelationsAcOnly( 
				accum_weights,       // final float [][][]     num_acc,          // number of accumulated tiles [tilesY][tilesX][pair]
				fcorr_td_acc); // final float [][][][] fcorr_td_acc      // [tilesY][tilesX][pair][256] sparse transform domain representation of corr pairs 
		double [][][]     dcorr_tiles = (fclt_corr != null)? (new double [tp_tasks_ref[0].length][][]):null;
		image_dtt.clt_process_tl_correlations( // convert to pixel domain and process correlations already prepared in fcorr_td and/or fcorr_combo_td
				imgdtt_params, // clt_parameters.img_dtt,		   // final ImageDttParameters  imgdtt_params,   // Now just extra correlation parameters, later will include, most others
				fcorr_td_acc,		 	       // final float  [][][][]     fcorr_td,        // [tilesY][tilesX][pair][4*64] transform domain representation of all selected corr pairs
				accum_weights,                 // float [][][]                num_acc,         // number of accumulated tiles [tilesY][tilesX][pair] (or null)       
				null, // dcorr_weight,         // double []                 dcorr_weight,    // alternative to num_acc, compatible with CPU processing (only one non-zero enough)
				clt_parameters.gpu_corr_scale, //  final double              gpu_corr_scale,  //  0.75; // reduce GPU-generated correlation values
				//				clt_parameters.getGpuFatZero(center_CLT.isMonochrome()),   // final double     gpu_fat_zero,    // clt_parameters.getGpuFatZero(is_mono);absolute == 30.0
				rng_fz,                        // final double     gpu_fat_zero,    // clt_parameters.getGpuFatZero(is_mono);absolute == 30.0
				image_dtt.transform_size - 1,  // final int                 gpu_corr_rad,    // = transform_size - 1 ?
				// The tp_tasks data should be decoded from GPU to get coordinates
				tp_tasks_ref[0],               // final TpTask []           tp_tasks,        // data from the reference frame - will be applied to LMW for the integrated correlations
				null,                          // final double [][]         far_fgbg,        // null, or [nTile]{disp(fg)-disp(bg), str(fg)-str(bg)} hints for LMA FG/BG split 
				center_CLT.getErsCorrection().getRXY(false), // final double [][]         rXY,             // from geometryCorrection
				// next both can be nulls
				null,                          // final double [][][][]     clt_corr_out,   // sparse (by the first index) [type][tilesY][tilesX][(2*transform_size-1)*(2*transform_size-1)] or null
				// combo will be added as extra pair if mcorr_comb_width > 0 and clt_corr_out has a slot for it
				// to be converted to float
				dcorr_tiles,                   // final double  [][][]      dcorr_tiles,     // [tile][pair][(2*transform_size-1)*(2*transform_size-1)] // if null - will not calculate
				// When clt_mismatch is non-zero, no far objects extraction will be attempted
				use_rms,                       // final boolean             use_rms, // DISPARITY_STRENGTH_INDEX means LMA RMS (18/04/2023)
				//optional, may be null
				disparity_map,                 // final double [][]         disparity_map,   // [8][tilesY][tilesX], only [6][] is needed on input or null - do not calculate
				null,                          // final double [][]         ddnd,            // data for LY. SHould be either null or [num_sensors][]
				clt_parameters.correlate_lma,  // final boolean             run_lma,         // calculate LMA, false - CM only
				// define combining of all 2D correlation pairs for CM (LMA does not use them)
				imgdtt_params.mcorr_comb_width, //final int                 mcorr_comb_width,  // combined correlation tile width (set <=0 to skip combined correlations)
				imgdtt_params.mcorr_comb_height,//final int                 mcorr_comb_height, // combined correlation tile full height
				imgdtt_params.mcorr_comb_offset,//final int                 mcorr_comb_offset, // combined correlation tile height offset: 0 - centered (-height/2 to height/2), height/2 - only positive (0 to height)
				imgdtt_params.mcorr_comb_disp,	 //final double              mcorr_comb_disp,   // Combined tile per-pixel disparity for baseline == side of a square
				clt_parameters.clt_window,     // final int                 window_type,     // GPU: will not be used
				clt_parameters.tileX,          // final int                 debug_tileX,
				clt_parameters.tileY,          // final int                 debug_tileY,
				ImageDtt.THREADS_MAX,                    // final int                 threadsMax,      // maximal number of threads to launch
				"accumulated",                 // final String              debug_suffix,
				debugLevel + 0); // 2); // -1 );     // final int                 globalDebugLevel)
		ImageDtt.convertFcltCorr(
				dcorr_tiles, // double [][][] dcorr_tiles,// [tile][sparse, correlation pair][(2*transform_size-1)*(2*transform_size-1)] // if null - will not calculate
				fclt_corr);  // float  [][][] fclt_corr) //  new float [tilesX * tilesY][][] or null

		if (show_corr || (accum_2d_corr != null)){ // -1
			float [][] accum_2d_img = ImageDtt.corr_partial_dbg( // not used in lwir
					fclt_corr, // final float  [][][]     fcorr_data,       // [tile][pair][(2*transform_size-1)*(2*transform_size-1)] // if null - will not calculate
					tp_tasks_ref[0], // final TpTask []         tp_tasks,        //
					tilesX,    //final int               tilesX,
					tilesY,    //final int               tilesX,
					2*image_dtt.transform_size - 1,	// final int               corr_size,
					1000, // will be limited by available layersfinal int               layers0,
					clt_parameters.corr_border_contrast, // final double            border_contrast,
					ImageDtt.THREADS_MAX, // final int               threadsMax,     // maximal number of threads to launch
					debugLevel); // final int               globalDebugLevel)
			String [] titles = new String [accum_2d_img.length]; // dcorr_tiles[0].length];
			int ind_length = image_dtt.getCorrelation2d().getCorrTitles().length;
			System.arraycopy(image_dtt.getCorrelation2d().getCorrTitles(), 0, titles, 0, ind_length);
			String [] combos= {"combo-all", "combo-diameters"};
			for (int i = ind_length; i < titles.length; i++) {
				titles[i] = combos[i-ind_length]; //  "combo-"+(i - ind_length);
			}
			if ((accum_2d_corr != null)) {
				accum_2d_corr[0] = accum_2d_img;
			}
			if (show_accumulated_correlations) {
				ShowDoubleFloatArrays.showArrays( // out of boundary 15
						accum_2d_img,
						tilesX*(2*image_dtt.transform_size),
						tilesY*(2*image_dtt.transform_size),
						true,
						center_CLT.getImageName()+"-FZ"+clt_parameters.imp.cuas_rng_fz+
						"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma +  
						"-CORR-ACCUM_R"+clt_parameters.imp.cuas_rng_radius+
						"-T"+ntile0+":"+nrefine+"-M"+mcorr_sel+"-"+clt_parameters.imp.cuas_mcorr_sel_lma,
						titles);      // image_dtt.getCorrelation2d().getCorrTitles()); //CORR_TITLES);
			}
		}
		if (rng_vfy && show_alt_vfy) {
			String [] titles_sensors = new String[img_um_seq[0].length];
			for (int i = 0; i < titles_sensors.length; i++) titles_sensors[i] = "SENS-"+i;
			String [] titles_scenes = new String[img_um_seq.length];
			for (int i = 0; i < titles_scenes.length; i++) {
				titles_scenes[i] = cuasMotion.getSceneTitles()[start_scene+i];
			}
			boolean   show = false; // true;
			double [][][] alt_images = new double [titles_scenes.length][][];
			for (int i = 0; i < alt_images.length; i++) {
				double [][][] alt_sens_col = scenes[start_scene+i].getAltImageData();
				alt_images[i] = new double[alt_sens_col.length][];
				for (int j = 0; j < alt_sens_col.length; j++) {
					alt_images[i][j] = alt_sens_col[j][0];
				}
			}
			//img_um_seq
			ImagePlus imp_alt =ShowDoubleFloatArrays.showArraysHyperstack(
					alt_images,       // double[][][] pixels, 
					center_CLT.getWidth(),         // int          width, 
					center_CLT.getImageName()+"-SENSOR_ALT_IMAGES"+
					//"-R"+"-FZ"+clt_parameters.imp.cuas_rng_fz+
					//"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+
					"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma  
		            //"-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
		            //"-ID"+ target_id+":"+nrefine+
		            //"-M"+clt_parameters.imp.cuas_mcorr_sel+
		            //"-"+clt_parameters.imp.cuas_mcorr_sel_lma
		            ,
		            titles_sensors,         // String []    titles, // all slices*frames titles or just slice titles or null
		            titles_scenes, // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
					show);          // boolean      show)
        	center_CLT.saveImagePlusInModelDirectory(
        			null,            // "GPU-SHIFTED-D"+clt_parameters.disparity, // String      suffix,
        			imp_alt); // imp_scenes); // ImagePlus   imp)
		}
		/*
		if (rng_vfy && show_sensors_vfy) { // verify img_um_seq did not change - yes, same as before
			String [] titles_sensors = new String[img_um_seq[0].length];
			for (int i = 0; i < titles_sensors.length; i++) titles_sensors[i] = "SENS-"+i;
			String [] titles_scenes = new String[img_um_seq.length];
			for (int i = 0; i < titles_scenes.length; i++) {
				titles_scenes[i] = cuasMotion.getSceneTitles()[start_scene+i];
			}
			boolean   show = false; // true;
			ImagePlus imp_sensors_um =ShowDoubleFloatArrays.showArraysHyperstack(
					img_um_seq,       // double[][][] pixels, 
					center_CLT.getWidth(),         // int          width, 
					center_CLT.getImageName()+"-SENSOR_FULL_IMAGES2"+
					//"-R"+"-FZ"+clt_parameters.imp.cuas_rng_fz+
					//"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+
					"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma  
		            //"-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
		            //"-ID"+ target_id+":"+nrefine+
		            //"-M"+clt_parameters.imp.cuas_mcorr_sel+
		            //"-"+clt_parameters.imp.cuas_mcorr_sel_lma
		            ,
		            titles_sensors,         // String []    titles, // all slices*frames titles or just slice titles or null
		            titles_scenes, // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
					show);          // boolean      show)
        	center_CLT.saveImagePlusInModelDirectory(
        			null,            // "GPU-SHIFTED-D"+clt_parameters.disparity, // String      suffix,
        			imp_sensors_um); // imp_scenes); // ImagePlus   imp)
		}
		*/
		return disparity_map; // disparity_map
	}
	

	public double [][] refineTargetDisparity( // returns disparity_map
			final double [][][] targets,   // centers
			final float [][][]  accum_2d_corr, // if [1][][] - return accumulated 2d correlations (all pairs)
			final double [][][] ranging_verify, // Now per group (rng_combine) [2*half_accum_range +1 +1][num_sens+1][] // [num_pix]
			final int           nseq,
			final int           nrefine, // number of the refine to select radius (larger first time)
			final boolean       show_corr,
			final boolean       show_debug_single,
			final int           debugLevel) {

		final double  um_sigma =    clt_parameters.imp.cuas_rng_um_sigma;
		final double  rng_fz =      (nrefine == 0)?clt_parameters.imp.cuas_rng_fz0 : clt_parameters.imp.cuas_rng_fz;
		final boolean mb_en =       clt_parameters.imp.mb_en; //  && (fov_tiles==null) && (mode3d > 0);
		final double  mb_tau =      clt_parameters.imp.mb_tau;      // 0.008; // time constant, sec
		final 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

		final double gpu_sigma_corr =     clt_parameters.getGpuCorrSigma(center_CLT.isMonochrome());
		final double gpu_sigma_rb_corr =  center_CLT.isMonochrome()? 1.0 : clt_parameters.gpu_sigma_rb_corr;
		final double gpu_sigma_log_corr = clt_parameters.getGpuCorrLoGSigma(center_CLT.isMonochrome());
		final double radius =       (nrefine == 0)? clt_parameters.imp.cuas_rng_radius0: clt_parameters.imp.cuas_rng_radius;

		final int rng_combine = clt_parameters.imp.cuas_rng_combine; // // combine multiple scenes before intrascene correlation

		final int frame_center =     cuasMotion.getFrameCenter(nseq);
		final int half_accum_range = cuasMotion.getSeqLength()/2;
		final int start_scene = frame_center - half_accum_range;
		final int num_scenes = 2*half_accum_range + 1;
		boolean   clean_up = true;  // remove unneded data
		
		boolean   erase_source = true; // check everything works
		
		
		final double [][][] img_um_seq = getUMSequence(
				start_scene,  // final int       start_scene,
				num_scenes,   // final int       num_scenes,
				clean_up,     // final boolean   clean_up, // remove unneded data
				erase_source, // final boolean   erase_source, // erase source images (will not clean up this.img_um)
				debugLevel);  // final int debugLevel)

		double [][][] pXpYDs = cuasMotion.targetPxPyD(
				targets[nseq]); // final double [][] targets)

		///   	boolean batch_run =clt_parameters.batch_run; // may be modified for debug
		///		boolean cuas_debug =  clt_parameters.imp.cuas_debug;  // save debug images (and show them if not in batch mode)
		final int tilesX =     center_CLT.getTilesX();
		final int tilesY =     center_CLT.getTilesY();
		final int sensor_mask = -1; // all sensors
		final int num_sens =   center_CLT.getNumSensors();

		ImageDttParameters imgdtt_params = null;
		try {
			imgdtt_params = clt_parameters.img_dtt.clone();
		} catch (CloneNotSupportedException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		imgdtt_params.updateFromCuas(clt_parameters.imp, num_sens);



		boolean show_um = show_debug_single; // (show_corr & 4) != 0; // true; //  false; // true;
		if (show_um) {
			
			String [] titles_sensors = new String[img_um_seq[0].length];
			String [] titles_scenes =  new String[num_scenes];
			for (int i = 0; i < titles_sensors.length; i++) {
				titles_sensors[i] = "SENS-"+i;
			}
			for (int i = 0; i < num_scenes; i++) {
				titles_scenes[i] = cuasMotion.getSceneTitles()[i+start_scene];
			}

			ShowDoubleFloatArrays.showArraysHyperstack(
					img_um_seq,       // double[][][] pixels, 
					center_CLT.getWidth(),         // int          width, 
					center_CLT.getImageName()+"-img_um_seq-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +um_sigma,          // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_sensors,           // String []    titles, // all slices*frames titles or just slice titles or null
					titles_scenes,            // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
					true); // show);          // boolean      show)
		}


		final int mcorr_sel = ImageDttParameters.corrSelEncode(imgdtt_params,num_sens);
		ErsCorrection ers_reference = center_CLT.getErsCorrection();
		int num_used_sens = 0;
		for (int i = 0; i < num_sens; i++) if (((sensor_mask >> i) & 1) != 0) num_used_sens++;
		int [] channels = new int [num_used_sens];
		int nch = 0;
		for (int i = 0; i < num_sens; i++) if (((sensor_mask >> i) & 1) != 0) channels[nch++] = i;

		// from OF13477: 	public static double[][] correlateInterscene(
		final int num_pairs = Correlation2d.getNumPairs(num_sens);
		final float  [][][][]   fcorr_td_acc  = new float [tilesY][tilesX][][];
		final float  [][][]     accum_weights = new float [tilesY][tilesX][num_pairs];
		final boolean show_2d_corr = false;
		boolean show_accumulated_correlations = show_2d_corr || debugLevel > -5;
		//		boolean show_reference_correlations =  show_2d_corr || debugLevel > -5;
		//		final float  [][][]       fclt_corr = ((accum_2d_corr != null) || show_accumulated_correlations || show_reference_correlations) ?
		final float  [][][]       fclt_corr = ((accum_2d_corr != null) || show_accumulated_correlations) ?
				(new float [tilesX * tilesY][][]) : null;
		final boolean        use_rms = true; // DISPARITY_STRENGTH_INDEX means LMA RMS (18/04/2023) (from OF11168:		boolean        use_rms = true; )
		final boolean        no_map = false;

		ImageDtt image_dtt = new ImageDtt(
				center_CLT.getNumSensors(),
				clt_parameters.transform_size,
				imgdtt_params, // clt_parameters.img_dtt,
				center_CLT.isAux(),
				center_CLT.isMonochrome(),
				center_CLT.isLwir(),
				clt_parameters.getScaleStrength(center_CLT.isAux()),
				center_CLT.getGPU());
		image_dtt.getCorrelation2d(); // initiate image_dtt.correlation2d, needed if disparity_map != null  
		final double[][] disparity_map = no_map ? null : new double [image_dtt.getDisparityTitles().length][];

		TpTask[][] tp_tasks = null; // will contain last tp_tasks
		TpTask[][] tp_tasks_ref = null;
		final double [] window_nnorm= cuasMotion.getSegmentWindow( // normalized, sum == 1.0
				clt_parameters.imp.cuas_rng_coswnd, // boolean smooth)
				false); // boolean normalize) // not normalized, values up to 1.0

		final int num_scene_grp = (num_scenes + rng_combine - 1) / rng_combine;
		final double [] window_grp = new double [num_scene_grp];
		double sw = 0;
		for (int nscene = 0; nscene < num_scenes; nscene++) {
			double w = window_nnorm[nscene];
			window_grp[nscene / rng_combine] += w;
			sw += w;
		}
		//		final double sum_window = sw;

		for (int dseq = -half_accum_range; dseq <= half_accum_range; dseq++) {
			final double [][] ref_pXpYD =    pXpYDs[dseq + half_accum_range]; // center target
			final int nscene = frame_center + dseq;
			final int iscene = dseq + half_accum_range; // 0, 1,..., 2* half_accum_range.
			final int iscene_grp = iscene / rng_combine;
			final boolean first_in_grp = (iscene % rng_combine) == 0;
			final boolean last_in_grp =  (((iscene + 1) % rng_combine) == 0) || (dseq == half_accum_range);

			final int sm = -1; // merge_all? -1: sensor_mask;
			final double [][] dxyzatr_dt = (mb_en ? new double[][] { // for all, including ref
				scenes[nscene].getErsCorrection().getErsXYZ_dt(),
				scenes[nscene].getErsCorrection().getErsATR_dt()} : null);
			final String ts = scenes[nscene].getImageName();
			final double [] scene_xyz = ers_reference.getSceneXYZ(ts);
			final double [] scene_atr = ers_reference.getSceneATR(ts);
			if ((scene_atr==null) || (scene_xyz == null)) {  // should not happen
				System.out.println("renderKeyFrame() BUG1");
				continue;
			}
			prepareAltImages(
					image_dtt,        //  final ImageDtt      image_dtt,
					ref_pXpYD,        // final double [][]   ref_pXpYD,
					img_um_seq,       // final double [][][] img_um_seq,
					scene_xyz,        // final double []     scene_xyz,       
					scene_atr,        // final double []     scene_atr,
					radius,           // final double        radius,
					-half_accum_range,// final int           first_offs,
					nseq,             // final int           nseq,
					dseq,             // final int           dseq,
					debugLevel);      // final int           debugLevel)


			tp_tasks = new TpTask[(mb_en && (dxyzatr_dt != null)) ? 2 : 1][];
			if (mb_en && (dxyzatr_dt != null)) {
				double [][] motion_blur = OpticalFlow.getMotionBlur(
						center_CLT, // quadCLTs[ref_index],   // QuadCLT        ref_scene,
						scenes[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)
				// TODO: Add target movement (recalculate velocities) 

				QuadCLT.preRenderGPUFromDSI(
						!first_in_grp,    // final boolean     accumulate,
						window_nnorm[iscene],// final double      global_scale, // <=1.0
						sm,               // final int         sensor_mask,
						false,            // final boolean     merge_channels,
						0,                // final int         discard_border,
						0,                // final double      max_fold,
						0,                // final int         min_in_row_col,   // Minimal number of defined tiles in a row/column
						null,             // final Rectangle   full_woi_in,      // show larger than sensor WOI in tiles (or null)
						clt_parameters,   // CLTParameters     clt_parameters,
						null,             // double []         disparity_ref,
						// only center tiles here!
						ref_pXpYD, // 5x5,     // 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,  // now [2][ntiles];

						scene_xyz,        // double []         scene_xyz, // camera center in world coordinates. If null - no shift, no ers
						scene_atr,        // double []         scene_atr, // camera orientation relative to world frame
						scenes[nscene],   //final QuadCLT     scene,
						center_CLT,       // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
						tp_tasks, 			  // TpTask [][]       tasks,
						false, // clt_parameters.imp.show_mono_nan,         // final boolean     show_nan,
						ImageDtt.THREADS_MAX, // threadsMax,       // int               threadsMax,
						debugLevel);      // final int         debugLevel)
			} else {
				QuadCLT.preRenderGPUFromDSI(
						!first_in_grp,    // final  boolean     accumulate,
						window_nnorm[iscene], // final double      global_scale, // <=1.0
						sm,               // final int         sensor_mask,
						false,            // final boolean     merge_channels,
						0,                // final int         discard_border,
						0,                // final double      max_fold,
						0,                // final int         min_in_row_col,   // Minimal number of defined tiles in a row/column
						null,             // final Rectangle   full_woi_in,      // show larger than sensor WOI in tiles (or null)
						clt_parameters,   // CLTParameters     clt_parameters,
						null,             // double []         disparity_ref,
						// only center tiles here!
						ref_pXpYD, // 5x5,     // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
						// motion blur compensation 
						0.0,              // double            mb_tau,      // 0.008; // time constant, sec
						0.0,              // 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];
						scene_xyz,        // double []         scene_xyz, // camera center in world coordinates. If null - no shift, no ers
						scene_atr,        // double []         scene_atr, // camera orientation relative to world frame
						scenes[nscene],   //final QuadCLT     scene,
						center_CLT,       // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
						tp_tasks, 		  // TpTask [][]       tasks,
						false, // clt_parameters.imp.show_mono_nan,         // final boolean     show_nan,
						ImageDtt.THREADS_MAX, // threadsMax,       // int               threadsMax,
						debugLevel);      // final int         debugLevel)
			}
			// copy center tasks to the reference
			if (dseq==0) {
				tp_tasks_ref = tp_tasks; // for the middle frame
				for (int i = 0; i < tp_tasks.length; i++) {
					tp_tasks_ref[i] = tp_tasks[i].clone();
				}
			}
			// correlate only after the last in group, others - just accumulate in TD
			if (last_in_grp) {
				if (ranging_verify != null) { // [2*half_accum_range +1 +1][num_sens+1][] // [num_pix])
					double [][] scene_render = center_CLT.renderDoubleFromTDMono ( // [sensor][pixel]
							-1,     // int                 sensor_mask,
							null,   // int []              wh, // may be null, or {width, height}
							false); // boolean             use_reference
					for (int nsens = 0; nsens < scene_render.length; nsens++) {
						ranging_verify[iscene_grp][nsens] = scene_render[nsens];
					}
				}

				float  [][][][]     fcorr_td =       new float[tilesY][tilesX][][];
				image_dtt.quadCorrTD(
						imgdtt_params, // clt_parameters.img_dtt,            // final ImageDttParameters imgdtt_params,    // Now just extra correlation parameters, later will include, most others
						tp_tasks[0], // *** will be updated inside from GPU-calculated geometry
						fcorr_td, // fcorrs_td[nscene],                 // [tilesY][tilesX][pair][4*64] transform domain representation of 6 corr pairs
						clt_parameters.gpu_sigma_r,        // 0.9, 1.1
						clt_parameters.gpu_sigma_b,        // 0.9, 1.1
						clt_parameters.gpu_sigma_g,        // 0.6, 0.7
						clt_parameters.gpu_sigma_m,        //  =       0.4; // 0.7;
						gpu_sigma_rb_corr,                 // final double              gpu_sigma_rb_corr, //  = 0.5; // apply LPF after accumulating R and B correlation before G, monochrome ? 1.0 : gpu_sigma_rb_corr;
						gpu_sigma_corr,                    //  =    0.9;gpu_sigma_corr_m
						gpu_sigma_log_corr,                // final double              gpu_sigma_log_corr,   // hpf to reduce dynamic range for correlations
						clt_parameters.corr_red,           // +used
						clt_parameters.corr_blue,          // +used
						mcorr_sel,                         // final int                 mcorr_sel,    // Which pairs to correlate // +1 - all, +2 - dia, +4 - sq, +8 - neibs, +16 - hor + 32 - vert
						ImageDtt.THREADS_MAX,       // maximal number of threads to launch
						debugLevel);
				if (image_dtt.getGPU().getGpu_debug_level() > -1) {
					System.out.println("==ooo=after image_dtt.quadCorrTD()");
				}
				// Verify tasks are now updated
				OpticalFlow.accumulateCorrelations(
						//						window_full[iscene] * (2 * half_accum_range + 1),        // final double         weight, // keep same sum as before
						window_grp[iscene_grp],        // final double         weight, // keep same sum as before
						accum_weights, // final int [][][]     num_acc,          // number of accumulated tiles [tilesY][tilesX][pair]
						fcorr_td,      // final float [][][][] fcorr_td,         // [tilesY][tilesX][pair][256] sparse transform domain representation of corr pairs 
						fcorr_td_acc); // final float [][][][] fcorr_td_acc      // [tilesY][tilesX][pair][256] sparse transform domain representation of corr pairs
				if (image_dtt.getGPU().getGpu_debug_level() > -1) {
					System.out.println("==ooo=accumulateCorrelations()");
				}
			}
		} // for (int dseq = half_accum_range; dseq <= half_accum_range; dseq++) {
		//boolean show_um=true;
		if (show_um) {
			// show again - is it different?
			/*
			{
				String [] titles_sensors = new String[img_um_seq[0].length];
				String [] titles_scenes =  new String[num_scenes];
				for (int i = 0; i < titles_sensors.length; i++) {
					titles_sensors[i] = "SENS-"+i;
				}
				for (int i = 0; i < num_scenes; i++) {
					titles_scenes[i] = cuasMotion.getSceneTitles()[i+start_scene];
				}

				ShowDoubleFloatArrays.showArraysHyperstack(
						img_um_seq,       // double[][][] pixels, 
						center_CLT.getWidth(),         // int          width, 
						center_CLT.getImageName()+"-img_um_seq-AFTER-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +um_sigma,          // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
						titles_sensors,           // String []    titles, // all slices*frames titles or just slice titles or null
						titles_scenes,            // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
						true); // show);          // boolean      show)

			}
			 */

			double [][][] alt_images =     new double [2*half_accum_range+1][num_sens][];
			double [][][] orig_images =     new double [2*half_accum_range+1][num_sens][];
			for (int i = 0; i < alt_images.length; i++) {
				double [][][] alt_col_img = scenes[i + frame_center-half_accum_range].getAltImageData();
				double [][][] orig_col_img = scenes[i + frame_center-half_accum_range].getOrigImageData();
				for (int j = 0; j < num_sens; j++) {
					alt_images[i][j] = alt_col_img[j][0];
					orig_images[i][j] = orig_col_img[j][0];
				}
			}

			String [] titles_alt_scenes =  new String[alt_images.length];

			String [] titles_sensors = new String[img_um_seq[0].length];
			for (int i = 0; i < titles_sensors.length; i++) {
				titles_sensors[i] = "SENS-"+i;
			}
			for (int i = 0; i < titles_alt_scenes.length; i++) {
				titles_alt_scenes[i] = cuasMotion.getSceneTitles()[i];
			}

			ShowDoubleFloatArrays.showArraysHyperstack(
					alt_images,       // double[][][] pixels, 
					center_CLT.getWidth(),         // int          width, 
					center_CLT.getImageName()+"-alt_images-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +um_sigma+"-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
					"-B"+clt_parameters.imp.cuas_rng_blur,          // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_sensors,                // String []    titles, // all slices*frames titles or just slice titles or null
					titles_alt_scenes,             // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
					true); // show);               // boolean      show)
			ShowDoubleFloatArrays.showArraysHyperstack(
					orig_images,       // double[][][] pixels, 
					center_CLT.getWidth(),         // int          width, 
					center_CLT.getImageName()+"-orig_images-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +um_sigma+"-R"+clt_parameters.imp.cuas_rng_radius0+":"+clt_parameters.imp.cuas_rng_radius+
					"-B"+clt_parameters.imp.cuas_rng_blur,          // String       title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
					titles_sensors,                // String []    titles, // all slices*frames titles or just slice titles or null
					titles_alt_scenes,             // CuasMotionLMA.LMA_TITLES,  // String []    frame_titles, // frame titles or null
					true); // show);               // boolean      show)
		}

		// Normalize accumulated correlations
		OpticalFlow.accumulateCorrelationsAcOnly( 
				accum_weights,       // final float [][][]     num_acc,          // number of accumulated tiles [tilesY][tilesX][pair]
				fcorr_td_acc); // final float [][][][] fcorr_td_acc      // [tilesY][tilesX][pair][256] sparse transform domain representation of corr pairs 
		double [][][]     dcorr_tiles = (fclt_corr != null)? (new double [tp_tasks_ref[0].length][][]):null;
		image_dtt.clt_process_tl_correlations( // convert to pixel domain and process correlations already prepared in fcorr_td and/or fcorr_combo_td
				imgdtt_params, // clt_parameters.img_dtt,		   // final ImageDttParameters  imgdtt_params,   // Now just extra correlation parameters, later will include, most others
				fcorr_td_acc,		 	       // final float  [][][][]     fcorr_td,        // [tilesY][tilesX][pair][4*64] transform domain representation of all selected corr pairs
				accum_weights,                 // float [][][]                num_acc,         // number of accumulated tiles [tilesY][tilesX][pair] (or null)       
				null, // dcorr_weight,         // double []                 dcorr_weight,    // alternative to num_acc, compatible with CPU processing (only one non-zero enough)
				clt_parameters.gpu_corr_scale, //  final double              gpu_corr_scale,  //  0.75; // reduce GPU-generated correlation values
				//				clt_parameters.getGpuFatZero(center_CLT.isMonochrome()),   // final double     gpu_fat_zero,    // clt_parameters.getGpuFatZero(is_mono);absolute == 30.0
				rng_fz,                        // final double     gpu_fat_zero,    // clt_parameters.getGpuFatZero(is_mono);absolute == 30.0
				image_dtt.transform_size - 1,  // final int                 gpu_corr_rad,    // = transform_size - 1 ?
				// The tp_tasks data should be decoded from GPU to get coordinates
				tp_tasks_ref[0],               // final TpTask []           tp_tasks,        // data from the reference frame - will be applied to LMW for the integrated correlations
				null,                          // final double [][]         far_fgbg,        // null, or [nTile]{disp(fg)-disp(bg), str(fg)-str(bg)} hints for LMA FG/BG split 
				center_CLT.getErsCorrection().getRXY(false), // final double [][]         rXY,             // from geometryCorrection
				// next both can be nulls
				null,                          // final double [][][][]     clt_corr_out,   // sparse (by the first index) [type][tilesY][tilesX][(2*transform_size-1)*(2*transform_size-1)] or null
				// combo will be added as extra pair if mcorr_comb_width > 0 and clt_corr_out has a slot for it
				// to be converted to float
				dcorr_tiles,                   // final double  [][][]      dcorr_tiles,     // [tile][pair][(2*transform_size-1)*(2*transform_size-1)] // if null - will not calculate
				// When clt_mismatch is non-zero, no far objects extraction will be attempted
				use_rms,                       // final boolean             use_rms, // DISPARITY_STRENGTH_INDEX means LMA RMS (18/04/2023)
				//optional, may be null
				disparity_map,                 // final double [][]         disparity_map,   // [8][tilesY][tilesX], only [6][] is needed on input or null - do not calculate
				null,                          // final double [][]         ddnd,            // data for LY. SHould be either null or [num_sensors][]
				clt_parameters.correlate_lma,  // final boolean             run_lma,         // calculate LMA, false - CM only
				// define combining of all 2D correlation pairs for CM (LMA does not use them)
				imgdtt_params.mcorr_comb_width, //final int                 mcorr_comb_width,  // combined correlation tile width (set <=0 to skip combined correlations)
				imgdtt_params.mcorr_comb_height,//final int                 mcorr_comb_height, // combined correlation tile full height
				imgdtt_params.mcorr_comb_offset,//final int                 mcorr_comb_offset, // combined correlation tile height offset: 0 - centered (-height/2 to height/2), height/2 - only positive (0 to height)
				imgdtt_params.mcorr_comb_disp,	 //final double              mcorr_comb_disp,   // Combined tile per-pixel disparity for baseline == side of a square
				clt_parameters.clt_window,     // final int                 window_type,     // GPU: will not be used
				clt_parameters.tileX,          // final int                 debug_tileX,
				clt_parameters.tileY,          // final int                 debug_tileY,
				ImageDtt.THREADS_MAX,                    // final int                 threadsMax,      // maximal number of threads to launch
				"accumulated",                 // final String              debug_suffix,
				debugLevel + 0); // 2); // -1 );     // final int                 globalDebugLevel)
		ImageDtt.convertFcltCorr(
				dcorr_tiles, // double [][][] dcorr_tiles,// [tile][sparse, correlation pair][(2*transform_size-1)*(2*transform_size-1)] // if null - will not calculate
				fclt_corr);  // float  [][][] fclt_corr) //  new float [tilesX * tilesY][][] or null

		//		if (show_accumulated_correlations || (accum_2d_corr != null)){ // -1
		if (show_corr || (accum_2d_corr != null)){ // -1
			float [][] accum_2d_img = ImageDtt.corr_partial_dbg( // not used in lwir
					fclt_corr, // final float  [][][]     fcorr_data,       // [tile][pair][(2*transform_size-1)*(2*transform_size-1)] // if null - will not calculate
					tp_tasks_ref[0], // final TpTask []         tp_tasks,        //
					tilesX,    //final int               tilesX,
					tilesY,    //final int               tilesX,
					2*image_dtt.transform_size - 1,	// final int               corr_size,
					1000, // will be limited by available layersfinal int               layers0,
					clt_parameters.corr_border_contrast, // final double            border_contrast,
					ImageDtt.THREADS_MAX, // final int               threadsMax,     // maximal number of threads to launch
					debugLevel); // final int               globalDebugLevel)
			String [] titles = new String [accum_2d_img.length]; // dcorr_tiles[0].length];
			int ind_length = image_dtt.getCorrelation2d().getCorrTitles().length;
			System.arraycopy(image_dtt.getCorrelation2d().getCorrTitles(), 0, titles, 0, ind_length);
			String [] combos= {"combo-all", "combo-diameters"};
			for (int i = ind_length; i < titles.length; i++) {
				titles[i] = combos[i-ind_length]; //  "combo-"+(i - ind_length);
			}
			if ((accum_2d_corr != null)) {
				accum_2d_corr[0] = accum_2d_img;
			}
			if (show_accumulated_correlations) {
				ShowDoubleFloatArrays.showArrays( // out of boundary 15
						accum_2d_img,
						tilesX*(2*image_dtt.transform_size),
						tilesY*(2*image_dtt.transform_size),
						true,
						center_CLT.getImageName()+"-FZ"+clt_parameters.imp.cuas_rng_fz+
						"-CB"+clt_parameters.imp.cuas_rng_combine+(clt_parameters.imp.cuas_rng_coswnd?"-WC":"-WR")+"-UM"+(clt_parameters.imp.cuas_rng_um2?"D":"") +clt_parameters.imp.cuas_rng_um_sigma +  
						"-CORR-ACCUM_R"+clt_parameters.imp.cuas_rng_radius+
						"-F"+frame_center+":"+nrefine+"-M"+mcorr_sel+"-"+clt_parameters.imp.cuas_mcorr_sel_lma,
						titles);      // image_dtt.getCorrelation2d().getCorrTitles()); //CORR_TITLES);
			}
		}
		return disparity_map; // disparity_map
	}
	
	public void prepareAltImages( // Sets single scene defined by  nseq and dseq
			final ImageDtt      image_dtt,
			final double [][]   ref_pXpYD,
			final double [][][] img_um_seq,
			final double []     scene_xyz,       
			final double []     scene_atr,
			final double        radius,
			final int           first_offs, // not used
			final int           nseq, // absolute sequence number (to calculate center frame
			final int           dseq,
			final int           debugLevel) {
		final double radius_blur =  clt_parameters.imp.cuas_rng_blur;
		final double alt_scale =    clt_parameters.imp.cuas_rng_scale; // scale alt data to approximately match correlation values
		
		final int tilesX =     center_CLT.getTilesX();
		final int tileSize =   center_CLT.getTileSize();
		final int width = tilesX * tileSize;
		final int frame_center =      cuasMotion.getFrameCenter(nseq);
//		final int half_accum_range =  cuasMotion.getSeqLength()/2;
		ErsCorrection ers_reference = center_CLT.getErsCorrection();
		final int num_sens =          center_CLT.getNumSensors();

		final int nscene = frame_center + dseq;
		final int iscene = dseq - first_offs; // starts from 0
		
		final String ts = scenes[nscene].getImageName();
		final double []   scene_ers_xyz_dt = ers_reference.getSceneErsXYZ_dt(ts);
		final double []   scene_ers_atr_dt = ers_reference.getSceneErsATR_dt(ts);
		scenes[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 radius >0, apply mask (around getSensorsXY()). If 0 - bypass
		double [][][] img_sens_col = new double [num_sens][1][];
		if (radius > 0) {
			double [][][] sensors_xy = getSensorsXY(  //[tile][sensor]{x,y}
					clt_parameters,    // CLTParameters     clt_parameters,
					image_dtt,          // ImageDtt          image_dtt,
					ref_pXpYD,          // double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
					scene_xyz,          // double []         scene_xyz, // camera center in world coordinates. If null - no shift, no ers
					scene_atr,          // double []         scene_atr, // camera orientation relative to world frame
					scenes[nscene],     // final QuadCLT     scene,
					center_CLT,         // final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
					debugLevel);        // final int         debugLevel);
			double [][] masked = copyMaskedTargets ( // [sensor][pixel]
					radius,             // final double        radius,
					radius_blur,        // final double        radius_blur,
					alt_scale,          // final double        alt_scale,
					img_um_seq[iscene], // final double [][]   channel_img,
					sensors_xy,         // final double [][][] sensors_xy, // [tile][sensor]{x,y}
					width,              // final int           width,
					debugLevel);        // final int           debugLevel);
			for (int nsens = 0; nsens < num_sens; nsens++) {
				img_sens_col[nsens][0] = masked[nsens]; // .clone();
			}
		} else {
			for (int nsens = 0; nsens < num_sens; nsens++) {
				img_sens_col[nsens][0] = img_um_seq[iscene][nsens]; //.clone(); // Index 80 out of bounds for length 71
			}
		}
		scenes[nscene].setImageDataAlt(img_sens_col);
		return; //  true;	
	}
	
	
	
	public static double [][][] getSensorsXY( // will return PxPyD
			CLTParameters     clt_parameters,
			ImageDtt          image_dtt,
			double [][]       ref_pXpYD,    // alternative to disparity_ref when reference is not uniform
  			// motion blur compensation 
			double []         scene_xyz, // camera center in world coordinates. If null - no shift, no ers
			double []         scene_atr, // camera orientation relative to world frame
			final QuadCLT     scene,
			final QuadCLT     ref_scene, // now - may be null - for testing if scene is rotated ref
			final int         debugLevel){
		double [][] pXpYD;
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		// window in pixels!
			pXpYD=OpticalFlow.transformToScenePxPyD(
					ref_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
					ref_scene, // final QuadCLT     reference_QuadClt)
					scene);    // final QuadCLT     scene_QuadClt)
			scene.getErsCorrection().setupERS(); // NEW - did not help
//			ref_scene.getErsCorrection().setupERS(); // just in case - setup using instance parameters - inside
		
		int rendered_width = scene.getErsCorrection().getSensorWH()[0];
		final TpTask[] tp_tasks =  GpuQuad.setInterTasks( // "true" reference, with stereo actual reference will be offset
                scene.getNumSensors(),
                rendered_width,           // should match output size, pXpYD.length
                !scene.hasGPU(),          // final boolean             calcPortsCoordinatesAndDerivatives, // GPU can calculate them centreXY
                pXpYD,                    // final double [][]         pXpYD, // per-tile array of pX,pY,disparity triplets (or nulls)
                null,                     // final boolean []          selection, // may be null, if not null do not  process unselected tiles
                scene.getErsCorrection(), // final GeometryCorrection  geometryCorrection,
//              clt_parameters.imp.disparity_corr, // 04/07/2023 // 0.0,                      // final double              disparity_corr,
				clt_parameters.imp.disparity_corr + scene.getDispInfinityRef(), // 12/11/2025 - added ref_scene.getDispInfinityRef()
                -1, // 0, // margin,      // final int                 margin,      // do not use tiles if their centers are closer to the edges
                null,                     // final boolean []          valid_tiles,            
                ImageDtt.THREADS_MAX);              // final int                 threadsMax)  // maximal number of threads to launch
//	    scene.saveQuadClt(); // to re-load new set of Bayer images to the GPU (do nothing for CPU) and Geometry
	    /*
	     * ImageDtt image_dtt = new ImageDtt(
	    		scene.getNumSensors(),
	    		clt_parameters.transform_size,
	    		clt_parameters.img_dtt,
	    		scene.isAux(),
	    		scene.isMonochrome(),
	    		scene.isLwir(),
	    		clt_parameters.getScaleStrength(scene.isAux()),
	    		scene.getGPU());
	    		*/
//	    boolean use_reference = false;
	    image_dtt.preSetReferenceTD( // do not run execConvertDirect, exit after updating tasks
	    		clt_parameters.img_dtt,  // ,     // final ImageDttParameters  imgdtt_params,    // Now just extra correlation parameters, later will include, most others
				tp_tasks,                // final TpTask[]            tp_tasks,
				debugLevel);             // final int                 globalDebugLevel)
	    final double [][][] sensors_xy = new double [tp_tasks.length][][];
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nTile= ai.getAndIncrement(); nTile < tp_tasks.length; nTile = ai.getAndIncrement()) if ((tp_tasks[nTile] != null) && (tp_tasks[nTile].getTask() != 0)){
						sensors_xy[nTile] = tp_tasks[nTile].getDoubleXY();
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
	    return sensors_xy;
	}
	
	public static double [][]  copyMaskedTargets ( // [sensor][pixel]
			final double        radius,
			final double        radius_blur,
			final double        alt_scale,
			final double [][]   channel_img,
			final double [][][] sensors_xy, // [tile][sensor]{x,y}
			final int           width,
			final int           debugLevel){
		final int num_sens =  channel_img.length;
		final int num_pix =   channel_img[0].length;
		final int num_tiles = sensors_xy.length;
		final int height = channel_img[0].length / width;
		final double [][] masked_img = new double [num_sens][num_pix];
		final Thread[] threads = ImageDtt.newThreadArray();
		final AtomicInteger ai = new AtomicInteger(0);
		final double inner_radius = radius * (1.0 - radius_blur);
		final double outer_radius = radius * (1.0 + radius_blur);
		final double inner_radius2 = inner_radius * inner_radius;
		final double outer_radius2 = outer_radius * outer_radius;
		for (int ithread = 0; ithread < threads.length; ithread++) {
			threads[ithread] = new Thread() {
				public void run() {
					for (int nTile= ai.getAndIncrement(); nTile < num_tiles; nTile = ai.getAndIncrement()) if (sensors_xy[nTile] != null) {
						double [][] sens_xy = sensors_xy[nTile];
						for (int nsens = 0; nsens < num_sens; nsens++) {
							double xc = sens_xy[nsens][0];
							double yc = sens_xy[nsens][1];
							int py0 = Math.max((int) Math.floor(yc - outer_radius), 0);
							int py1 = Math.min((int) Math.ceil (yc + outer_radius), height-1);
							int px0 = Math.max((int) Math.floor(xc - outer_radius), 0);
							int px1 = Math.min((int) Math.ceil (xc + outer_radius), width-1);
							for (int py = py0; py <= py1; py++) {
								double dy = py - yc;
								double dy2 = dy*dy;
								for (int px = px0; px <= px1; px++) {
									int npix = py*width + px;
									double dx = px - xc;
									double r2 = dy2+dx*dx;
									if (r2 <= inner_radius2) {
										masked_img[nsens][npix] = channel_img[nsens][npix];
									} else if (r2 <= outer_radius2) {
										double r = Math.sqrt(r2);
										double k =  0.5 * (1.0 + Math.cos(Math.PI*(r - inner_radius) / (outer_radius - inner_radius)));
										masked_img[nsens][npix] = k* channel_img[nsens][npix];
									}
								}
							}
						}
					}
				}
			};
		}		      
		ImageDtt.startAndJoin(threads);
		return masked_img;
	}
	
	public void addUasData(
			final double [][][] targets_single) {
		int num_seq = targets_single.length;
		UasLogReader  uasLogReader = cuasMotion.getUasLogReader();
		String [] slice_titles = 	cuasMotion.getSliceTitles(); // timestamps
		int tilesX =     cuasMotion.getTilesX();
		int tilesY =    targets_single[0].length / tilesX;
		for (int nseq = 0; nseq < num_seq; nseq++) {
			String timestamp = slice_titles[nseq];
			double [] uas_pXpYDRange = uasLogReader.getUasPxPyDRange(timestamp); // px, py, d-  cuas_infinity (true disparity), range
			if (uas_pXpYDRange != null) {
				double px =        uas_pXpYDRange[0];
				double py =        uas_pXpYDRange[1];
				double disparity = uas_pXpYDRange[2]; // true disparity, w/o infinity offset
				double range =     uas_pXpYDRange[3];
				int tileX = (int) (px/GPUTileProcessor.DTT_SIZE);
				int tileY = (int) (py/GPUTileProcessor.DTT_SIZE);
				if ((tileX>=0) && (tileY>=0) && (tileX < tilesX)  && (tileY < tilesY)) {
					int ntile = tileX + tileY * tilesX;
					if (targets_single[nseq][ntile] == null) {
						targets_single[nseq][ntile] = CuasMotionLMA.getEmpty(); // all fields NaN
						targets_single[nseq][ntile][CuasMotionLMA.RSLT_FAIL] = CuasMotionLMA.FAIL_FL_ONLY; // should be enough not to treat as a target
					}
					double [] tile = targets_single[nseq][ntile];
					tile[CuasMotionLMA.RSLT_FL_PX] =    px;
					tile[CuasMotionLMA.RSLT_FL_PY] =    py;
					tile[CuasMotionLMA.RSLT_FL_DISP] =  disparity;
					tile[CuasMotionLMA.RSLT_FL_RANGE] = range;
				}
			}
		}		
	}
	
	// relies on calcMatchingTargetsLengths(.., true,...) called from recalcOmegas() to set [RSLT_GLOBAL] 
	public void saveTargetStats(
			final double [][][] targets_single) {
		final int tilesX = cuasMotion.getTilesX();
		final GeometryCorrection gc = center_CLT.getGeometryCorrection();
		final int tileSize = GPUTileProcessor.DTT_SIZE;
		double  cuas_infinity =      clt_parameters.imp.cuas_infinity ; // 0.63;  // disparity at infinity for targets
		int      flw_levels = clt_parameters.imp.cuas_flw_levels; // 1 - all, 2 - all and two halves, 3 - all, two halves and 4 quaters
		
		int num_seq = targets_single.length;
		int [][][] linked_targets = CuasMotion.getLinkedTargets(targets_single);
		int num_targets = linked_targets.length;
		StringBuffer sb = new StringBuffer();
		sb.append(getRangingLogParameters());
		String []  fields_uas = {"pX", "pY"," azimuth", "elevation","disparity","range"};
		String []  fields =     {"pX", "pY"," azimuth", "elevation","disparity","range","glength", "gdisp","gdisp-diff","gstr","grange"};
//		String []  global_fields = {"glength", "gdisp","gdisp-diff","gstr","grange"};
		// First header line
		sb.append("scene\t\tFlight log"); for (int i = 0; i < fields_uas.length; i++) sb.append("\t");
		for (int ntarg = 0; ntarg < num_targets; ntarg++) {
			sb.append("Target "+(ntarg+1));for (int i = 0; i < fields.length; i++) sb.append("\t"); // there will be 1 extra blank column
		}
		sb.append("\n"); // there will be 1 extra blank column
		// Second header line
		sb.append("index\ttimestamp\t");
		for (int ntarg = -1; ntarg < num_targets; ntarg++) {
			String suffix = (ntarg < 0) ? "-fl_log":("-"+(ntarg+1));
			if (ntarg < 0) {
				for (int i = 0; i < fields_uas.length; i++) {
					sb.append(fields[i]+suffix+"\t");
				}
			} else {
				for (int i = 0; i < fields.length; i++) {
					sb.append(fields[i]+suffix+"\t");
				}
			}
		}
		sb.append("\n"); // there will be 1 extra blank column
		String [] slice_titles = 	cuasMotion.getSliceTitles(); // timestamps
		UasLogReader  uasLogReader = cuasMotion.getUasLogReader();
		ErsCorrection ersCorrection = center_CLT.getErsCorrection();
		for (int nseq = 0; nseq < num_seq; nseq++) {
			String timestamp = slice_titles[nseq];
			sb.append(nseq+"\t"+timestamp+"\t");
			// get azimuth, elevation, target disparity from the log plus infinity, log range
			double [] uas_pXpYDRange = uasLogReader.getUasPxPyDRange(timestamp); // px, py, d-  cuas_infinity (true disparity), range
			double [][] az_el_oaz_oel= 	CuasMotion.getPixToAzElev(
					ersCorrection,   // ErsCorrection ersCorrection,
					uasLogReader,    // UasLogReader  uasLogReader,			
					Double.NaN,      // double        fps, // if NaN, will use default 60Hz. Used only for omegas
					uas_pXpYDRange[0],  // double         target_x,
					uas_pXpYDRange[1],  // double         target_y,
					0,                  // double         target_vx,
					0);                 // double         target_vy);
			
			sb.append(uas_pXpYDRange[0]+"\t"+uas_pXpYDRange[1]+"\t"+az_el_oaz_oel[0][0]+"\t"+az_el_oaz_oel[0][1]+"\t"+(uas_pXpYDRange[2]+cuas_infinity)+"\t"+uas_pXpYDRange[3]+"\t");

			for (int ntarg = 0; ntarg < num_targets; ntarg++) { 
				if (linked_targets[ntarg][nseq] != null) {
					int ntile = linked_targets[ntarg][nseq][0]; // [1] - ntarg not used here, it is always 0
					double [] target = targets_single[nseq][ntile];
					if (target == null) {
						System.out.println("saveTargetStats() BUG, (target == null). ntarg = "+ntarg+", nseq="+nseq+
								", linked_targets["+ntarg+"]["+nseq+"][0]="+linked_targets[ntarg][nseq][0]);
						System.out.println();
						for (int i = 0; i < fields.length; i++) {
							sb.append("\t");
						}
					} else {
						int tileX = ntile % tilesX;
						int tileY = ntile / tilesX;
					    double xc = tileSize * tileX + tileSize/2 + target[CuasMotionLMA.RSLT_X];  
					    double yc = tileSize * tileY + tileSize/2 + target[CuasMotionLMA.RSLT_Y];
						double vx = target[CuasMotionLMA.RSLT_VX];
						double vy = target[CuasMotionLMA.RSLT_VY];
						// calculate and output target azimuth, elevation, disparity (full) and range
						az_el_oaz_oel= 	CuasMotion.getPixToAzElev(
								ersCorrection, // ErsCorrection ersCorrection,
								uasLogReader,  // UasLogReader  uasLogReader,			
								Double.NaN,    // double        fps, // if NaN, will use default 60Hz. Used only for omegas
								xc,             // double         px, // null
								yc,             // double         py,
								vx,             // double         vx,
								vy);            // double         vy);
						
						sb.append(xc+"\t"+yc+"\t"+az_el_oaz_oel[0][0]+"\t"+az_el_oaz_oel[0][1]+"\t"+
								target[CuasMotionLMA.RSLT_DISPARITY]+"\t"+target[CuasMotionLMA.RSLT_RANGE]+"\t");
						if (target[CuasMotionLMA.RSLT_GLENGTH] > 0) {
							sb.append(target[CuasMotionLMA.RSLT_GLENGTH]+"\t"+
									target[CuasMotionLMA.RSLT_GDISPARITY]+"\t"+
									target[CuasMotionLMA.RSLT_GDISP_DIFF]+"\t"+
									target[CuasMotionLMA.RSLT_GDISP_STR]+"\t"+
									target[CuasMotionLMA.RSLT_GRANGE]+"\t");
						} else {
							for (int i = 0; i < (fields.length - fields_uas.length); i++) {
								sb.append("\t");
							}
						}
					}
				} else {
					for (int i = 0; i < fields.length; i++) {
						sb.append("\t");
					}
				}
			}
			sb.append("\n"); // there will be 1 extra blank column
		}
		String suffix = CuasMotion.getParametersSuffixRslt(clt_parameters, TARGET_STATS_SUFFIX)+".csv";
		center_CLT.saveStringInModelDirectory(sb.toString(), suffix, false);
		return;
	}


}
